diff --git a/examples/python/custom-material/bi-material.py b/examples/python/custom-material/bi-material.py index 017acc202..37bf307b4 100644 --- a/examples/python/custom-material/bi-material.py +++ b/examples/python/custom-material/bi-material.py @@ -1,192 +1,184 @@ from __future__ import print_function # ------------------------------------------------------------- # -import akantu +import akantu as aka import subprocess import numpy as np import time # ------------------------------------------------------------- # class LocalElastic: # declares all the internals def initMaterial(self, internals, params): self.E = params['E'] self.nu = params['nu'] self.rho = params['rho'] # print(self.__dict__) # First Lame coefficient self.lame_lambda = self.nu * self.E / ( (1. + self.nu) * (1. - 2. * self.nu)) # Second Lame coefficient (shear modulus) self.lame_mu = self.E / (2. * (1. + self.nu)) all_factor = internals['factor'] all_quad_coords = internals['quad_coordinates'] for elem_type in all_factor.keys(): factor = all_factor[elem_type] quad_coords = all_quad_coords[elem_type] factor[:] = 1. factor[quad_coords[:, 1] < 0.5] = .5 # declares all the internals @staticmethod def registerInternals(): return ['potential', 'factor'] # declares all the internals @staticmethod def registerInternalSizes(): return [1, 1] # declares all the parameters that could be parsed @staticmethod def registerParam(): return ['E', 'nu'] # declares all the parameters that are needed def getPushWaveSpeed(self, params): return np.sqrt((self.lame_lambda + 2 * self.lame_mu) / self.rho) # compute small deformation tensor @staticmethod def computeEpsilon(grad_u): return 0.5 * (grad_u + np.einsum('aij->aji', grad_u)) # constitutive law def computeStress(self, grad_u, sigma, internals, params): n_quads = grad_u.shape[0] grad_u = grad_u.reshape((n_quads, 2, 2)) factor = internals['factor'].reshape(n_quads) epsilon = self.computeEpsilon(grad_u) sigma = sigma.reshape((n_quads, 2, 2)) trace = np.einsum('aii->a', grad_u) sigma[:, :, :] = ( np.einsum('a,ij->aij', trace, self.lame_lambda * np.eye(2)) + 2. * self.lame_mu * epsilon) # print(sigma.reshape((n_quads, 4))) # print(grad_u.reshape((n_quads, 4))) sigma[:, :, :] = np.einsum('aij, a->aij', sigma, factor) # constitutive law tangent modulii def computeTangentModuli(self, grad_u, tangent, internals, params): n_quads = tangent.shape[0] tangent = tangent.reshape(n_quads, 3, 3) factor = internals['factor'].reshape(n_quads) Miiii = self.lame_lambda + 2 * self.lame_mu Miijj = self.lame_lambda Mijij = self.lame_mu tangent[:, 0, 0] = Miiii tangent[:, 1, 1] = Miiii tangent[:, 0, 1] = Miijj tangent[:, 1, 0] = Miijj tangent[:, 2, 2] = Mijij tangent[:, :, :] = np.einsum('aij, a->aij', tangent, factor) # computes the energy density def getEnergyDensity(self, energy_type, energy_density, grad_u, stress, internals, params): nquads = stress.shape[0] stress = stress.reshape(nquads, 2, 2) grad_u = grad_u.reshape((nquads, 2, 2)) if energy_type != 'potential': raise RuntimeError('not known energy') epsilon = self.computeEpsilon(grad_u) energy_density[:, 0] = ( 0.5 * np.einsum('aij,aij->a', stress, epsilon)) # applies manually the boundary conditions def applyBC(model): nbNodes = model.getMesh().getNbNodes() position = model.getMesh().getNodes() displacement = model.getDisplacement() blocked_dofs = model.getBlockedDOFs() width = 1. height = 1. epsilon = 1e-8 for node in range(0, nbNodes): if((np.abs(position[node, 0]) < epsilon) or # left side (np.abs(position[node, 0] - width) < epsilon)): # right side blocked_dofs[node, 0] = True displacement[node, 0] = 0 * position[node, 0] + 0. if(np.abs(position[node, 1]) < epsilon): # lower side blocked_dofs[node, 1] = True displacement[node, 1] = - 1. if(np.abs(position[node, 1] - height) < epsilon): # upper side blocked_dofs[node, 1] = True displacement[node, 1] = 1. - -################################################################ -# main -################################################################ - -# main paraeters +# main parameters spatial_dimension = 2 mesh_file = 'square.msh' # call gmsh to generate the mesh ret = subprocess.call(['gmsh', '-2', 'square.geo', '-optimize', 'square.msh']) if ret != 0: raise Exception( 'execution of GMSH failed: do you have it installed ?') time.sleep(1) # read mesh -mesh = akantu.Mesh(spatial_dimension) +mesh = aka.Mesh(spatial_dimension) mesh.read(mesh_file) # create the custom material mat = LocalElastic() -akantu.registerNewPythonMaterial(mat, "local_elastic") +aka.registerNewPythonMaterial(mat, "local_elastic") -# init akantu -akantu.initialize('material.dat') +# parse input file +aka.parseInput('material.dat') # init the SolidMechanicsModel -model = akantu.SolidMechanicsModel(mesh) -model.initFull(_analysis_method=akantu._static) +model = aka.SolidMechanicsModel(mesh) +model.initFull(_analysis_method=aka._static) # configure the solver solver = model.getNonLinearSolver() solver.set("max_iterations", 2) solver.set("threshold", 1e-3) -solver.set("convergence_type", akantu._scc_solution) +solver.set("convergence_type", aka._scc_solution) # prepare the dumper model.setBaseName("bimaterial") model.addDumpFieldVector("displacement") model.addDumpFieldVector("internal_force") model.addDumpFieldVector("external_force") model.addDumpField("strain") model.addDumpField("stress") model.addDumpField("factor") model.addDumpField("blocked_dofs") # Boundary conditions applyBC(model) # solve the problem model.solveStep() # dump paraview files model.dump() - -# shutdown akantu -akantu.finalize() diff --git a/examples/python/custom-material/custom-material.py b/examples/python/custom-material/custom-material.py index 199d51dbe..b93cd14e5 100644 --- a/examples/python/custom-material/custom-material.py +++ b/examples/python/custom-material/custom-material.py @@ -1,210 +1,207 @@ #!/usr/bin/env python3 from __future__ import print_function ################################################################ import os import subprocess import numpy as np import akantu ################################################################ class FixedValue: def __init__(self, value, axis): self.value = value self.axis = axis def operator(self, node, flags, disp, coord): # sets the displacement to the desired value in the desired axis disp[self.axis] = self.value # sets the blocked dofs vector to true in the desired axis flags[self.axis] = True ################################################################ class LocalElastic: # declares all the internals def initMaterial(self, internals, params): self.E = params['E'] self.nu = params['nu'] self.rho = params['rho'] # First Lame coefficient self.lame_lambda = self.nu * self.E / ( (1 + self.nu) * (1 - 2 * self.nu)) # Second Lame coefficient (shear modulus) self.lame_mu = self.E / (2 * (1 + self.nu)) # declares all the internals @staticmethod def registerInternals(): return ['potential'] # declares all the internals @staticmethod def registerInternalSizes(): return [1] # declares all the parameters that could be parsed @staticmethod def registerParam(): return ['E', 'nu'] # declares all the parameters that are needed def getPushWaveSpeed(self, params): return np.sqrt((self.lame_lambda + 2 * self.lame_mu) / self.rho) # compute small deformation tensor @staticmethod def computeEpsilon(grad_u): return 0.5 * (grad_u + np.einsum('aij->aji', grad_u)) # constitutive law def computeStress(self, grad_u, sigma, internals, params): nquads = grad_u.shape[0] grad_u = grad_u.reshape((nquads, 2, 2)) epsilon = self.computeEpsilon(grad_u) sigma = sigma.reshape((nquads, 2, 2)) trace = np.trace(grad_u, axis1=1, axis2=2) sigma[:, :, :] = ( np.einsum('a,ij->aij', trace, self.lame_lambda * np.eye(2)) + 2.*self.lame_mu * epsilon) # constitutive law tangent modulii def computeTangentModuli(self, grad_u, tangent, internals, params): n_quads = tangent.shape[0] tangent = tangent.reshape(n_quads, 3, 3) Miiii = self.lame_lambda + 2 * self.lame_mu Miijj = self.lame_lambda Mijij = self.lame_mu tangent[:, 0, 0] = Miiii tangent[:, 1, 1] = Miiii tangent[:, 0, 1] = Miijj tangent[:, 1, 0] = Miijj tangent[:, 2, 2] = Mijij # computes the energy density def getEnergyDensity(self, energy_type, energy_density, grad_u, stress, internals, params): nquads = stress.shape[0] stress = stress.reshape(nquads, 2, 2) grad_u = grad_u.reshape((nquads, 2, 2)) if energy_type != 'potential': raise RuntimeError('not known energy') epsilon = self.computeEpsilon(grad_u) energy_density[:, 0] = ( 0.5 * np.einsum('aij,aij->a', stress, epsilon)) ################################################################ # main ################################################################ spatial_dimension = 2 -akantu.initialize('material.dat') +akantu.parseInput('material.dat') mesh_file = 'bar.msh' max_steps = 250 time_step = 1e-3 # if mesh was not created the calls gmsh to generate it if not os.path.isfile(mesh_file): ret = subprocess.call('gmsh -2 bar.geo bar.msh', shell=True) if ret != 0: raise Exception( 'execution of GMSH failed: do you have it installed ?') ################################################################ # Initialization ################################################################ mesh = akantu.Mesh(spatial_dimension) mesh.read(mesh_file) mat = LocalElastic() akantu.registerNewPythonMaterial(mat, "local_elastic") model = akantu.SolidMechanicsModel(mesh) model.initFull(_analysis_method=akantu._explicit_lumped_mass) # model.initFull(_analysis_method=akantu._implicit_dynamic) model.setBaseName("waves") model.addDumpFieldVector("displacement") model.addDumpFieldVector("acceleration") model.addDumpFieldVector("velocity") model.addDumpFieldVector("internal_force") model.addDumpFieldVector("external_force") model.addDumpField("strain") model.addDumpField("stress") model.addDumpField("blocked_dofs") ################################################################ # boundary conditions ################################################################ model.applyDirichletBC(FixedValue(0, akantu._x), "XBlocked") model.applyDirichletBC(FixedValue(0, akantu._y), "YBlocked") ################################################################ # initial conditions ################################################################ displacement = model.getDisplacement() nb_nodes = mesh.getNbNodes() position = mesh.getNodes() pulse_width = 1 A = 0.01 for i in range(0, nb_nodes): # Sinus * Gaussian x = position[i, 0] - 5. L = pulse_width k = 0.1 * 2 * np.pi * 3 / L displacement[i, 0] = A * \ np.sin(k * x) * np.exp(-(k * x) * (k * x) / (L * L)) ################################################################ # timestep value computation ################################################################ time_factor = 0.8 stable_time_step = model.getStableTimeStep() * time_factor print("Stable Time Step = {0}".format(stable_time_step)) print("Required Time Step = {0}".format(time_step)) time_step = stable_time_step * time_factor model.setTimeStep(time_step) ################################################################ # loop for evolution of motion dynamics ################################################################ model.assembleInternalForces() print("step,step * time_step,epot,ekin,epot + ekin") for step in range(0, max_steps + 1): model.solveStep() if step % 10 == 0: model.dump() epot = model.getEnergy('potential') ekin = model.getEnergy('kinetic') # output energy calculation to screen print("{0},{1},{2},{3},{4}".format(step, step * time_step, epot, ekin, (epot + ekin))) - -akantu.finalize() - diff --git a/examples/python/dynamics/dynamics.py b/examples/python/dynamics/dynamics.py index 23d34898e..c8be57fec 100644 --- a/examples/python/dynamics/dynamics.py +++ b/examples/python/dynamics/dynamics.py @@ -1,130 +1,129 @@ #!/usr/bin/env python3 from __future__ import print_function ################################################################ import os import subprocess import numpy as np import akantu ################################################################ class FixedValue: def __init__(self, value, axis): self.value = value self.axis = axis def operator(self, node, flags, disp, coord): # sets the displacement to the desired value in the desired axis disp[self.axis] = self.value # sets the blocked dofs vector to true in the desired axis flags[self.axis] = True ################################################################ def main(): spatial_dimension = 2 - akantu.initialize('material.dat') + akantu.parseInput('material.dat') mesh_file = 'bar.msh' max_steps = 250 time_step = 1e-3 # if mesh was not created the calls gmsh to generate it if not os.path.isfile(mesh_file): ret = subprocess.call('gmsh -2 bar.geo bar.msh', shell=True) if ret != 0: raise Exception( 'execution of GMSH failed: do you have it installed ?') ################################################################ # Initialization ################################################################ mesh = akantu.Mesh(spatial_dimension) mesh.read(mesh_file) model = akantu.SolidMechanicsModel(mesh) model.initFull(_analysis_method=akantu._explicit_lumped_mass) # model.initFull(_analysis_method=akantu._implicit_dynamic) model.setBaseName("waves") model.addDumpFieldVector("displacement") model.addDumpFieldVector("acceleration") model.addDumpFieldVector("velocity") model.addDumpFieldVector("internal_force") model.addDumpFieldVector("external_force") model.addDumpField("strain") model.addDumpField("stress") model.addDumpField("blocked_dofs") ################################################################ # boundary conditions ################################################################ model.applyDirichletBC(FixedValue(0, akantu._x), "XBlocked") model.applyDirichletBC(FixedValue(0, akantu._y), "YBlocked") ################################################################ # initial conditions ################################################################ displacement = model.getDisplacement() nb_nodes = mesh.getNbNodes() position = mesh.getNodes() pulse_width = 1 A = 0.01 for i in range(0, nb_nodes): # Sinus * Gaussian x = position[i, 0] - 5. L = pulse_width k = 0.1 * 2 * np.pi * 3 / L displacement[i, 0] = A * \ np.sin(k * x) * np.exp(-(k * x) * (k * x) / (L * L)) ################################################################ # timestep value computation ################################################################ time_factor = 0.8 stable_time_step = model.getStableTimeStep() * time_factor print("Stable Time Step = {0}".format(stable_time_step)) print("Required Time Step = {0}".format(time_step)) time_step = stable_time_step * time_factor model.setTimeStep(time_step) ################################################################ # loop for evolution of motion dynamics ################################################################ model.assembleInternalForces() print("step,step * time_step,epot,ekin,epot + ekin") for step in range(0, max_steps + 1): model.solveStep() if step % 10 == 0: model.dump() epot = model.getEnergy('potential') ekin = model.getEnergy('kinetic') # output energy calculation to screen print("{0},{1},{2},{3},{4}".format(step, step * time_step, epot, ekin, (epot + ekin))) - akantu.finalize() return ################################################################ if __name__ == "__main__": main() diff --git a/examples/python/plate-hole/plate.py b/examples/python/plate-hole/plate.py index bc821c0ea..bafded612 100644 --- a/examples/python/plate-hole/plate.py +++ b/examples/python/plate-hole/plate.py @@ -1,113 +1,112 @@ #!/usr/bin/env python3 from __future__ import print_function import akantu import numpy as np ################################################################ # Dirichlet Boudary condition functor: fix the displacement ################################################################ class FixedValue: def __init__(self, value, axis): self.value = value self.axis = axis def operator(self, node, flags, disp, coord): # sets the displacement to the desired value in the desired axis disp[self.axis] = self.value # sets the blocked dofs vector to true in the desired axis flags[self.axis] = True ################################################################ # Neumann Boudary condition functor: impose a traction ################################################################ class FromTraction: def __init__(self, traction): self.traction = traction def operator(self, quad_point, force, coord, normals): # sets the force to the desired value in the desired axis force[:] = self.traction ################################################################ def solve(material_file, mesh_file, traction): - akantu.initialize(material_file) + akantu.parseInput(material_file) spatial_dimension = 2 ################################################################ # Initialization ################################################################ mesh = akantu.Mesh(spatial_dimension) mesh.read(mesh_file) model = akantu.SolidMechanicsModel(mesh) model.initFull(akantu.SolidMechanicsModelOptions(akantu._static)) model.setBaseName("plate") model.addDumpFieldVector("displacement") model.addDumpFieldVector("external_force") model.addDumpField("strain") model.addDumpField("stress") model.addDumpField("blocked_dofs") ################################################################ # Boundary conditions ################################################################ model.applyDirichletBC(FixedValue(0.0, akantu._x), "XBlocked") model.applyDirichletBC(FixedValue(0.0, akantu._y), "YBlocked") trac = np.zeros(spatial_dimension) trac[1] = traction print("Solve for traction ", traction) model.getForce()[:] = 0 model.applyNeumannBC(FromTraction(trac), "Traction") solver = model.getNonLinearSolver() solver.set("max_iterations", 2) solver.set("threshold", 1e-10) solver.set("convergence_type", akantu._scc_residual) model.solveStep() model.dump() - akantu.finalize() ################################################################ # main ################################################################ def main(): import os mesh_file = 'plate.msh' # if mesh was not created the calls gmsh to generate it if not os.path.isfile(mesh_file): import subprocess ret = subprocess.call( 'gmsh -2 plate.geo {0}'.format(mesh_file), shell=True) if not ret == 0: raise Exception( 'execution of GMSH failed: do you have it installed ?') material_file = 'material.dat' traction = 1. solve(material_file, mesh_file, traction) ################################################################ if __name__ == "__main__": main() diff --git a/python/swig/aka_common.i b/python/swig/aka_common.i index 7a6f7c384..f25044211 100644 --- a/python/swig/aka_common.i +++ b/python/swig/aka_common.i @@ -1,143 +1,157 @@ /** * @file aka_common.i * * @author Guillaume Anciaux * @author Fabian Barras * @author Nicolas Richart * * @date creation: Fri Dec 12 2014 * @date last modification: Wed Jan 13 2016 * * @brief wrapper to aka_common.hh * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ %{ //#include "aka_config.hh" #include "aka_common.hh" #include "aka_csr.hh" #include "element.hh" #include "python_functor.hh" #include "parser.hh" %} -%include "stl.i" -%include "parser.i" - namespace akantu { %ignore getUserParser; %ignore ghost_types; %ignore initialize(int & argc, char ** & argv); %ignore initialize(const std::string & input_file, int & argc, char ** & argv); + %ignore finalize; extern const Array empty_filter; } +%include "stl.i" +%include "parser.i" + + %typemap(in) (int argc, char *argv[]) { int i = 0; if (!PyList_Check($input)) { PyErr_SetString(PyExc_ValueError, "Expecting a list"); return NULL; } $1 = PyList_Size($input); $2 = new char *[$1+1]; for (i = 0; i < $1; i++) { PyObject *s = PyList_GetItem($input,i); if (!PyString_Check(s)) { free($2); PyErr_SetString(PyExc_ValueError, "List items must be strings"); return NULL; } $2[i] = PyString_AsString(s); } $2[i] = 0; } %typemap(freearg) (int argc, char *argv[]) { %#if defined(__INTEL_COMPILER) //#pragma warning ( disable : 383 ) %#elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu %#elif (defined(__GNUC__) || defined(__GNUG__)) %# if __cplusplus > 199711L %# pragma GCC diagnostic ignored "-Wmaybe-uninitialized" %# endif %#endif delete [] $2; %#if defined(__INTEL_COMPILER) //#pragma warning ( disable : 383 ) %#elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu %#elif (defined(__GNUC__) || defined(__GNUG__)) %# if __cplusplus > 199711L %# pragma GCC diagnostic pop %# endif %#endif } %inline %{ namespace akantu { #if defined(AKANTU_USE_MPI) const int MPI=1; #else const int MPI=0; #endif void _initializeWithoutArgv(const std::string & input_file) { int nb_args = 0; char ** null = nullptr; initialize(input_file, nb_args, null); } + void __finalize() { + finalize(); + } + #define AKANTU_PP_STR_TO_TYPE2(s, data, elem) \ ({ BOOST_PP_STRINGIZE(elem) , elem}) PyObject * getElementTypes(){ - - std::map element_types{ - BOOST_PP_SEQ_FOR_EACH_I( + + std::map element_types{ + BOOST_PP_SEQ_FOR_EACH_I( AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(AKANTU_ek_regular_ELEMENT_TYPE), - BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_STR_TO_TYPE2, akantu, AKANTU_ek_regular_ELEMENT_TYPE))}; - + BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_STR_TO_TYPE2, akantu, AKANTU_ek_regular_ELEMENT_TYPE))}; + return akantu::PythonFunctor::convertToPython(element_types); } } - %} %pythoncode %{ import sys as _aka_sys - def initialize(input_file="", argv=_aka_sys.argv): + def __initialize(input_file="", argv=_aka_sys.argv): if _aka_sys.modules[__name__].MPI == 1: try: from mpi4py import MPI except ImportError: pass _initializeWithoutArgv(input_file) + def initialize(*args, **kwargs): + raise RuntimeError("No need to call initialize," + " use parseInput to read an input file") + + def finalize(*args, **kwargs): + raise RuntimeError("No need to call finalize") + + def parseInput(input_file): + getStaticParser().parse(input_file) %} %include "aka_config.hh" %include "aka_common.hh" %include "aka_element_classes_info.hh" %include "element.hh" %include "aka_error.hh" diff --git a/python/swig/akantu.i b/python/swig/akantu.i index aecdc6d80..abfd3456c 100644 --- a/python/swig/akantu.i +++ b/python/swig/akantu.i @@ -1,80 +1,84 @@ /** * @file akantu.i * * @author Guillaume Anciaux * @author Fabian Barras * @author Nicolas Richart * * @date creation: Fri Dec 12 2014 * @date last modification: Mon Nov 23 2015 * * @brief Main swig file for akantu' python interface * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ %module akantu %exception { try { $action } catch (akantu::debug::Exception e) { PyErr_SetString(PyExc_IndexError,e.what()); return NULL; } } #define __attribute__(x) %ignore akantu::operator <<; %include "aka_common.i" %include "aka_csr.i" %include "aka_array.i" %define print_self(MY_CLASS) %extend akantu::MY_CLASS { std::string __str__() { std::stringstream sstr; sstr << *($self); return sstr.str(); } } %enddef %include "mesh.i" %include "mesh_utils.i" %include "model.i" %include "communicator.i" %include "solid_mechanics_model.i" #if defined(AKANTU_COHESIVE_ELEMENT) %include "solid_mechanics_model_cohesive.i" #endif #if defined(AKANTU_HEAT_TRANSFER) %include "heat_transfer_model.i" #endif #if defined(AKANTU_STRUCTURAL_MECHANICS) %include "load_functions.i" %include "structural_mechanics_model.i" #endif + +%pythoncode %{ + __initialize() +%} diff --git a/python/swig/model.i b/python/swig/model.i index ca5f4f497..8f0353180 100644 --- a/python/swig/model.i +++ b/python/swig/model.i @@ -1,135 +1,132 @@ /** * @file model.i * * @author Guillaume Anciaux * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Dec 12 2014 * @date last modification: Wed Nov 11 2015 * * @brief model wrapper * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ %{ #include "boundary_condition_python_functor.hh" #include "model_solver.hh" #include "non_linear_solver.hh" %} namespace akantu { %ignore Model::createSynchronizerRegistry; %ignore Model::getSynchronizerRegistry; %ignore Model::createParallelSynch; %ignore Model::getDOFSynchronizer; %ignore Model::registerFEEngineObject; %ignore Model::unregisterFEEngineObject; %ignore Model::getFEEngineBoundary; // %ignore Model::getFEEngine; %ignore Model::getFEEngineClass; %ignore Model::getFEEngineClassBoundary; %ignore Model::setParser; %ignore Model::updateDataForNonLocalCriterion; %ignore IntegrationPoint::operator=; %ignore FEEngine::getNbIntegrationPoints; %ignore FEEngine::getShapes; %ignore FEEngine::getShapesDerivatives; %ignore FEEngine::getIntegrationPoints; %ignore FEEngine::getIGFEMElementTypes; %ignore FEEngine::interpolateOnIntegrationPoints(const Array &,ElementTypeMapArray &,const ElementTypeMapArray *) const; %ignore FEEngine::interpolateOnIntegrationPoints(const Array &,ElementTypeMapArray &) const; %ignore FEEngine::interpolateOnIntegrationPoints(const Array &,Array &,UInt,const ElementType&,const GhostType &,const Array< UInt > &) const; %ignore FEEngine::interpolateOnIntegrationPoints(const Array &,Array &,UInt,const ElementType&,const GhostType &) const; %ignore FEEngine::onNodesAdded; %ignore FEEngine::onNodesRemoved; %ignore FEEngine::onElementsAdded; %ignore FEEngine::onElementsChanged; %ignore FEEngine::onElementsRemoved; %ignore FEEngine::elementTypes; } %include "sparse_matrix.i" %include "fe_engine.hh" %rename(FreeBoundaryDirichlet) akantu::BC::Dirichlet::FreeBoundary; %rename(FreeBoundaryNeumann) akantu::BC::Neumann::FreeBoundary; %rename(PythonBoundary) akantu::BC::Dirichlet::PythonFunctor; %include "boundary_condition_functor.hh" %include "boundary_condition.hh" %include "boundary_condition_python_functor.hh" %include "communication_buffer.hh" %include "data_accessor.hh" //%include "synchronizer.hh" //%include "synchronizer_registry.hh" %include "model.hh" %include "non_linear_solver.hh" %include "model_options.hh" %extend akantu::Model { void initFullImpl( const akantu::ModelOptions & options = akantu::ModelOptions()){ $self->initFull(options); }; %insert("python") %{ def initFull(self, *args, **kwargs): if len(args) == 0: import importlib as __aka_importlib _module = __aka_importlib.import_module(self.__module__) _cls = getattr(_module, "{0}Options".format(self.__class__.__name__)) options = _cls() - dir(options) if len(kwargs) > 0: for key, val in kwargs.items(): if key[0] == '_': key = key[1:] - print(key, val) setattr(options, key, val) else: options = args[0] - print(options) self.initFullImpl(options) %} void solveStep(const akantu::ID & id=""){ $self->solveStep(id); } akantu::NonLinearSolver & getNonLinearSolver(const akantu::ID & id=""){ return $self->getNonLinearSolver(id); } } %extend akantu::NonLinearSolver { void set(const std::string & id, akantu::Real val){ if (id == "max_iterations") $self->set(id, int(val)); else if (id == "convergence_type") $self->set(id, akantu::SolveConvergenceCriteria(UInt(val))); else $self->set(id, val); } } diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index b208f073d..0edab13c7 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,561 +1,568 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Model of Solid Mechanics * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" #include "non_local_manager.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition, public NonLocalManagerCallback, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; using MyFEEngineType = FEEngineTemplate; protected: using EventManager = EventHandlerManager; public: SolidMechanicsModel( Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0, const ModelType model_type = ModelType::_solid_mechanics_model); ~SolidMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl( const ModelOptions & options = SolidMechanicsModelOptions()) override; /// initialize all internal arrays for materials virtual void initMaterials(); /// initialize the model void initModel() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this adds f_{ext} or f_{int} to the residual void assembleResidual(const ID & residual_part) override; bool canSplitResidual() override { return true; } /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep() override; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } protected: /// update the current position vector void updateCurrentPosition(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register an empty material of a given type Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param); /// reassigns materials depending on the material selector virtual void reassignMaterial(); /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials virtual void assignMaterialToElements(const ElementTypeMapArray * filter = nullptr); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(const ElementType & type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* NonLocalManager inherited members */ /* ------------------------------------------------------------------------ */ protected: void initializeNonLocal() override; void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; void computeNonLocalStresses(const GhostType & ghost_type) override; void insertIntegrationPointsInNeighborhoods(const GhostType & ghost_type) override; /// update the values of the non local internal void updateLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /// copy the results of the averaging in the materials void updateNonLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: void splitElementByMaterial(const Array & elements, std::vector> & elements_per_mat) const; template void splitByMaterial(const Array & elements, Operation && op) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array & nodes_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array &, const Array &, const ElementTypeMapArray &, const ChangedElementsEvent &) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); //! decide wether a field is a material internal or not bool isInternal(const std::string & field_name, const ElementKind & element_kind); #ifndef SWIG //! give the amount of data per element virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, const ElementKind & kind); //! flatten a given material internal field ElementTypeMapArray & flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(const ElementKind & kind); #endif dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the SolidMechanicsModel::previous_displacement vector AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition const Array & getCurrentPosition(); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *displacement_increment, Array &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, Array &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); + /// get the SolidMechanicsModel::external_force vector (external forces) + AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); + /// get the SolidMechanicsModel::force vector (external forces) - AKANTU_GET_MACRO(Force, *external_force, Array &); + Array & getForce() { + AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, " + "use getExternalForce instead"); + return *external_force; + } /// get the SolidMechanicsModel::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the SolidMechanicsModel::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get a particular material (by material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); } /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for energy Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); - AKANTU_GET_MACRO(MaterialLocalNumbering, - material_local_numbering, const ElementTypeMapArray &); + AKANTU_GET_MACRO(MaterialLocalNumbering, material_local_numbering, + const ElementTypeMapArray &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, material_local_numbering, UInt); AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, *material_selector, MaterialSelector &); AKANTU_SET_MACRO(MaterialSelector, material_selector, std::shared_ptr); /// Access the non_local_manager interface AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); /// get the FEEngine object to integrate or interpolate on the boundary FEEngine & getFEEngineBoundary(const ID & name = "") override; protected: friend class Material; protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array * displacement; UInt displacement_release{0}; /// displacements array at the previous time step (used in finite deformation) Array * previous_displacement{nullptr}; /// increment of displacement Array * displacement_increment{nullptr}; /// lumped mass array Array * mass{nullptr}; /// Check if materials need to recompute the mass array bool need_to_reassemble_lumped_mass{true}; /// Check if materials need to recompute the mass matrix bool need_to_reassemble_mass{true}; /// velocities array Array * velocity{nullptr}; /// accelerations array Array * acceleration{nullptr}; /// accelerations array // Array * increment_acceleration; /// external forces array Array * external_force{nullptr}; /// internal forces array Array * internal_force{nullptr}; /// array specifing if a degree of freedom is blocked or not Array * blocked_dofs{nullptr}; /// array of current position used during update residual Array * current_position{nullptr}; UInt current_position_release{0}; /// Arrays containing the material index for each element ElementTypeMapArray material_index; /// Arrays containing the position in the element filter of the material /// (material's local numbering) ElementTypeMapArray material_local_numbering; /// list of used materials std::vector> materials; /// mapping between material name and material internal id std::map materials_names_to_id; /// class defining of to choose a material std::shared_ptr material_selector; /// flag defining if the increment must be computed or not bool increment_flag; /// tells if the material are instantiated bool are_materials_instantiated; using flatten_internal_map = std::map, ElementTypeMapArray *>; /// map a registered internals to be flattened for dump purposes flatten_internal_map registered_internals; /// non local manager std::unique_ptr non_local_manager; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { using FromStress = FromHigherDim; using FromTraction = FromSameDim; } // namespace Neumann } // namespace BC } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" #include "solid_mechanics_model_inline_impl.cc" #include "solid_mechanics_model_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/src/solver/sparse_solver_mumps.cc b/src/solver/sparse_solver_mumps.cc index bd797a5f8..ba9350fd8 100644 --- a/src/solver/sparse_solver_mumps.cc +++ b/src/solver/sparse_solver_mumps.cc @@ -1,444 +1,451 @@ /** * @file sparse_solver_mumps.cc * * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Tue Feb 20 2018 * * @brief implem of SparseSolverMumps class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * @subsection Ctrl_param Control parameters * * ICNTL(1), * ICNTL(2), * ICNTL(3) : output streams for error, diagnostics, and global messages * * ICNTL(4) : verbose level : 0 no message - 4 all messages * * ICNTL(5) : type of matrix, 0 assembled, 1 elementary * * ICNTL(6) : control the permutation and scaling(default 7) see mumps doc for * more information * * ICNTL(7) : determine the pivot order (default 7) see mumps doc for more * information * * ICNTL(8) : describe the scaling method used * * ICNTL(9) : 1 solve A x = b, 0 solve At x = b * * ICNTL(10) : number of iterative refinement when NRHS = 1 * * ICNTL(11) : > 0 return statistics * * ICNTL(12) : only used for SYM = 2, ordering strategy * * ICNTL(13) : * * ICNTL(14) : percentage of increase of the estimated working space * * ICNTL(15-17) : not used * * ICNTL(18) : only used if ICNTL(5) = 0, 0 matrix centralized, 1 structure on * host and mumps give the mapping, 2 structure on host and distributed matrix * for facto, 3 distributed matrix * * ICNTL(19) : > 0, Shur complement returned * * ICNTL(20) : 0 rhs dense, 1 rhs sparse * * ICNTL(21) : 0 solution in rhs, 1 solution distributed in ISOL_loc and SOL_loc * allocated by user * * ICNTL(22) : 0 in-core, 1 out-of-core * * ICNTL(23) : maximum memory allocatable by mumps pre proc * * ICNTL(24) : controls the detection of "null pivot rows" * * ICNTL(25) : * * ICNTL(26) : * * ICNTL(27) : * * ICNTL(28) : 0 automatic choice, 1 sequential analysis, 2 parallel analysis * * ICNTL(29) : 0 automatic choice, 1 PT-Scotch, 2 ParMetis */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "dof_manager_default.hh" #include "dof_synchronizer.hh" #include "sparse_matrix_aij.hh" #if defined(AKANTU_USE_MPI) #include "mpi_communicator_data.hh" #endif #include "sparse_solver_mumps.hh" /* -------------------------------------------------------------------------- */ // static std::ostream & operator <<(std::ostream & stream, const DMUMPS_STRUC_C // & _this) { // stream << "DMUMPS Data [" << std::endl; // stream << " + job : " << _this.job << std::endl; // stream << " + par : " << _this.par << std::endl; // stream << " + sym : " << _this.sym << std::endl; // stream << " + comm_fortran : " << _this.comm_fortran << std::endl; // stream << " + nz : " << _this.nz << std::endl; // stream << " + irn : " << _this.irn << std::endl; // stream << " + jcn : " << _this.jcn << std::endl; // stream << " + nz_loc : " << _this.nz_loc << std::endl; // stream << " + irn_loc : " << _this.irn_loc << std::endl; // stream << " + jcn_loc : " << _this.jcn_loc << std::endl; // stream << "]"; // return stream; // } namespace akantu { /* -------------------------------------------------------------------------- */ SparseSolverMumps::SparseSolverMumps(DOFManagerDefault & dof_manager, const ID & matrix_id, const ID & id, const MemoryID & memory_id) : SparseSolver(dof_manager, matrix_id, id, memory_id), dof_manager(dof_manager), master_rhs_solution(0, 1) { AKANTU_DEBUG_IN(); this->prank = communicator.whoAmI(); #ifdef AKANTU_USE_MPI this->parallel_method = _fully_distributed; #else // AKANTU_USE_MPI this->parallel_method = _not_parallel; #endif // AKANTU_USE_MPI AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseSolverMumps::~SparseSolverMumps() { AKANTU_DEBUG_IN(); mumpsDataDestroy(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::mumpsDataDestroy() { +#ifdef AKANTU_USE_MPI + int finalized = 0; + MPI_Finalized(&finalized); + if (finalized) // Da fuck !? + return; +#endif + if (this->is_initialized) { this->mumps_data.job = _smj_destroy; // destroy dmumps_c(&this->mumps_data); this->is_initialized = false; } } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::destroyInternalData() { mumpsDataDestroy(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::checkInitialized() { if (this->is_initialized) return; this->initialize(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::setOutputLevel() { // Output setup icntl(1) = 0; // error output icntl(2) = 0; // diagnostics output icntl(3) = 0; // information icntl(4) = 0; #if !defined(AKANTU_NDEBUG) DebugLevel dbg_lvl = debug::debugger.getDebugLevel(); if (AKANTU_DEBUG_TEST(dblDump)) { strcpy(this->mumps_data.write_problem, "mumps_matrix.mtx"); } // clang-format off icntl(1) = (dbg_lvl >= dblWarning) ? 6 : 0; icntl(3) = (dbg_lvl >= dblInfo) ? 6 : 0; icntl(2) = (dbg_lvl >= dblTrace) ? 6 : 0; icntl(4) = dbg_lvl >= dblDump ? 4 : dbg_lvl >= dblTrace ? 3 : dbg_lvl >= dblInfo ? 2 : dbg_lvl >= dblWarning ? 1 : 0; // clang-format on #endif } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::initMumpsData() { auto & A = dof_manager.getMatrix(matrix_id); // Default Scaling icntl(8) = 77; // Assembled matrix icntl(5) = 0; /// Default centralized dense second member icntl(20) = 0; icntl(21) = 0; // automatic choice for analysis icntl(28) = 0; UInt size = A.size(); if (prank == 0) { this->master_rhs_solution.resize(size); } this->mumps_data.nz_alloc = 0; this->mumps_data.n = size; switch (this->parallel_method) { case _fully_distributed: icntl(18) = 3; // fully distributed this->mumps_data.nz_loc = A.getNbNonZero(); this->mumps_data.irn_loc = A.getIRN().storage(); this->mumps_data.jcn_loc = A.getJCN().storage(); break; case _not_parallel: case _master_slave_distributed: icntl(18) = 0; // centralized if (prank == 0) { this->mumps_data.nz = A.getNbNonZero(); this->mumps_data.irn = A.getIRN().storage(); this->mumps_data.jcn = A.getJCN().storage(); } else { this->mumps_data.nz = 0; this->mumps_data.irn = nullptr; this->mumps_data.jcn = nullptr; } break; default: AKANTU_ERROR("This case should not happen!!"); } } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::initialize() { AKANTU_DEBUG_IN(); this->mumps_data.par = 1; // The host is part of computations switch (this->parallel_method) { case _not_parallel: break; case _master_slave_distributed: this->mumps_data.par = 0; // The host is not part of the computations /* FALLTHRU */ /* [[fallthrough]]; un-comment when compiler will get it */ case _fully_distributed: #ifdef AKANTU_USE_MPI const auto & mpi_data = dynamic_cast( communicator.getCommunicatorData()); MPI_Comm mpi_comm = mpi_data.getMPICommunicator(); this->mumps_data.comm_fortran = MPI_Comm_c2f(mpi_comm); #else AKANTU_ERROR( "You cannot use parallel method to solve without activating MPI"); #endif break; } const auto & A = dof_manager.getMatrix(matrix_id); this->mumps_data.sym = 2 * (A.getMatrixType() == _symmetric); this->prank = communicator.whoAmI(); this->setOutputLevel(); this->mumps_data.job = _smj_initialize; // initialize dmumps_c(&this->mumps_data); this->setOutputLevel(); this->is_initialized = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::analysis() { AKANTU_DEBUG_IN(); initMumpsData(); this->mumps_data.job = _smj_analyze; // analyze dmumps_c(&this->mumps_data); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::factorize() { AKANTU_DEBUG_IN(); auto & A = dof_manager.getMatrix(matrix_id); if (parallel_method == _fully_distributed) this->mumps_data.a_loc = A.getA().storage(); else { if (prank == 0) this->mumps_data.a = A.getA().storage(); } this->mumps_data.job = _smj_factorize; // factorize dmumps_c(&this->mumps_data); this->printError(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::solve(Array & x, const Array & b) { auto & synch = this->dof_manager.getSynchronizer(); if (this->prank == 0) { this->master_rhs_solution.resize(this->dof_manager.getSystemSize()); synch.gather(b, this->master_rhs_solution); } else { synch.gather(b); } this->solveInternal(); if (this->prank == 0) { synch.scatter(x, this->master_rhs_solution); } else { synch.scatter(x); } } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::solve() { this->master_rhs_solution.copy(this->dof_manager.getGlobalResidual()); this->solveInternal(); this->dof_manager.setGlobalSolution(this->master_rhs_solution); this->dof_manager.splitSolutionPerDOFs(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::solveInternal() { AKANTU_DEBUG_IN(); this->checkInitialized(); const auto & A = dof_manager.getMatrix(matrix_id); this->setOutputLevel(); if (this->last_profile_release != A.getProfileRelease()) { this->analysis(); this->last_profile_release = A.getProfileRelease(); } if (AKANTU_DEBUG_TEST(dblDump)) { A.saveMatrix("solver_mumps" + std::to_string(prank) + ".mtx"); } if (this->last_value_release != A.getValueRelease()) { this->factorize(); this->last_value_release = A.getValueRelease(); } if (prank == 0) { this->mumps_data.rhs = this->master_rhs_solution.storage(); } this->mumps_data.job = _smj_solve; // solve dmumps_c(&this->mumps_data); this->printError(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::printError() { Vector _info_v(2); _info_v[0] = info(1); // to get errors _info_v[1] = -info(1); // to get warnings dof_manager.getCommunicator().allReduce(_info_v, SynchronizerOperation::_min); _info_v[1] = -_info_v[1]; if (_info_v[0] < 0) { // < 0 is an error switch (_info_v[0]) { case -10: { AKANTU_CUSTOM_EXCEPTION( debug::SingularMatrixException(dof_manager.getMatrix(matrix_id))); break; } case -9: { icntl(14) += 10; if (icntl(14) != 90) { // std::cout << "Dynamic memory increase of 10%" << std::endl; AKANTU_DEBUG_WARNING("MUMPS dynamic memory is insufficient it will be " "increased allowed to use 10% more"); // change releases to force a recompute this->last_value_release--; this->last_profile_release--; this->solve(); } else { AKANTU_ERROR("The MUMPS workarray is too small INFO(2)=" << info(2) << "No further increase possible"); } break; } default: AKANTU_ERROR("Error in mumps during solve process, check mumps " "user guide INFO(1) = " << _info_v[1]); } } else if (_info_v[1] > 0) { AKANTU_DEBUG_WARNING("Warning in mumps during solve process, check mumps " "user guide INFO(1) = " << _info_v[1]); } } } // akantu diff --git a/test/test_model/patch_tests/patch_test_linear_anisotropic_explicit.py b/test/test_model/patch_tests/patch_test_linear_anisotropic_explicit.py index f8b455dca..7ba974329 100644 --- a/test/test_model/patch_tests/patch_test_linear_anisotropic_explicit.py +++ b/test/test_model/patch_tests/patch_test_linear_anisotropic_explicit.py @@ -1,86 +1,85 @@ #!/usr/bin/env python3 # ------------------------------------------------------------------------------ __author__ = "Guillaume Anciaux" __copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ " en Mécanique des Solides)" __credits__ = ["Guillaume Anciaux"] __license__ = "L-GPLv3" __maintainer__ = "Guillaume Anciaux" __email__ = "guillaume.anciaux@epfl.ch" # ------------------------------------------------------------------------------ from patch_test_linear_solid_mechanics_fixture import TestPatchTestSMMLinear import akantu import unittest import numpy as np # Stiffness tensor, rotated by hand C = np.array([[[[112.93753505, 1.85842452538e-10, -4.47654358027e-10], [1.85847317471e-10, 54.2334345331, -3.69840984824], [-4.4764768395e-10, -3.69840984824, 56.848605217]], [[1.85847781609e-10, 25.429294233, -3.69840984816], [25.429294233, 3.31613847493e-10, -8.38797920011e-11], [-3.69840984816, -8.38804581349e-11, -1.97875715813e-10]], [[-4.47654358027e-10, -3.69840984816, 28.044464917], [-3.69840984816, 2.09374961813e-10, 9.4857455224e-12], [28.044464917, 9.48308098714e-12, -2.1367885239e-10]]], [[[1.85847781609e-10, 25.429294233, -3.69840984816], [25.429294233, 3.31613847493e-10, -8.38793479119e-11], [-3.69840984816, -8.38795699565e-11, -1.97876381947e-10]], [[54.2334345331, 3.31617400207e-10, 2.09372075233e-10], [3.3161562385e-10, 115.552705733, -3.15093728886e-10], [2.09372075233e-10, -3.15090176173e-10, 54.2334345333]], [[-3.69840984824, -8.38795699565e-11, 9.48219280872e-12], [-8.38795699565e-11, -3.1509195253e-10, 25.4292942335], [9.48441325477e-12, 25.4292942335, 3.69840984851]]], [[[-4.47653469848e-10, -3.69840984816, 28.044464917], [-3.69840984816, 2.09374073634e-10, 9.48752187924e-12], [28.044464917, 9.48552347779e-12, -2.1367885239e-10]], [[-3.69840984824, -8.3884899027e-11, 9.48219280872e-12], [-8.3884899027e-11, -3.150972816e-10, 25.4292942335], [9.48041645188e-12, 25.4292942335, 3.69840984851]], [[56.848605217, -1.97875493768e-10, -2.13681516925e-10], [-1.97877270125e-10, 54.2334345333, 3.69840984851], [-2.13683293282e-10, 3.69840984851, 112.93753505]]]]) def foo(self): self.initModel( akantu.SolidMechanicsModelOptions(akantu._explicit_lumped_mass), "material_anisotropic.dat") coordinates = self.mesh.getNodes() displacement = self.model.getDisplacement() # set the position of all nodes to the static solution self.setLinearDOF(displacement, coordinates) for s in range(0, 100): self.model.solveStep() ekin = self.model.getEnergy("kinetic") self.assertAlmostEqual(0, ekin, delta=1e-16) self.checkDisplacements() self.checkStrains() def foo(pstrain): strain = (pstrain + pstrain.transpose()) / 2. stress = np.zeros((self.dim, self.dim)) for i in range(0, self.dim): for j in range(0, self.dim): stress[i, j] = 0 for k in range(0, self.dim): for l in range(0, self.dim): stress[i, j] += C[i][j][k][l] * strain(k, l) return stress self.checkStresses(foo) -akantu.initialize() suite = TestPatchTestSMMLinear.TYPED_TEST(foo, "AnisotropicExplicit") unittest.TextTestRunner(verbosity=1).run(suite) diff --git a/test/test_model/patch_tests/patch_test_linear_elastic_explicit.py b/test/test_model/patch_tests/patch_test_linear_elastic_explicit.py index e777b6132..16e9470d6 100644 --- a/test/test_model/patch_tests/patch_test_linear_elastic_explicit.py +++ b/test/test_model/patch_tests/patch_test_linear_elastic_explicit.py @@ -1,41 +1,41 @@ #!/usr/bin/env python3 # ------------------------------------------------------------------------------ __author__ = "Guillaume Anciaux" __copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ " en Mécanique des Solides)" __credits__ = ["Guillaume Anciaux"] __license__ = "L-GPLv3" __maintainer__ = "Guillaume Anciaux" __email__ = "guillaume.anciaux@epfl.ch" # ------------------------------------------------------------------------------ from patch_test_linear_solid_mechanics_fixture import TestPatchTestSMMLinear import akantu def foo(self): filename = "material_check_stress_plane_stress.dat" if self.plane_strain: filename = "material_check_stress_plane_strain.dat" self.initModel( akantu.SolidMechanicsModelOptions(akantu._explicit_lumped_mass), filename) coordinates = self.mesh.getNodes() displacement = self.model.getDisplacement() # set the position of all nodes to the static solution self.setLinearDOF(displacement, coordinates) for s in range(0, 100): self.model.solveStep() ekin = self.model.getEnergy("kinetic") self.assertAlmostEqual(0, ekin, delta=1e-16) self.checkAll() -akantu.initialize() + TestPatchTestSMMLinear.TYPED_TEST(foo, "Explicit") diff --git a/test/test_model/patch_tests/patch_test_linear_elastic_static.py b/test/test_model/patch_tests/patch_test_linear_elastic_static.py index 23e2a7bfd..a3ccc8620 100644 --- a/test/test_model/patch_tests/patch_test_linear_elastic_static.py +++ b/test/test_model/patch_tests/patch_test_linear_elastic_static.py @@ -1,36 +1,36 @@ #!/usr/bin/env python3 # ------------------------------------------------------------------------------ __author__ = "Guillaume Anciaux" __copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ " en Mécanique des Solides)" __credits__ = ["Guillaume Anciaux"] __license__ = "L-GPLv3" __maintainer__ = "Guillaume Anciaux" __email__ = "guillaume.anciaux@epfl.ch" # ------------------------------------------------------------------------------ from patch_test_linear_solid_mechanics_fixture import TestPatchTestSMMLinear import akantu def foo(self): filename = "material_check_stress_plane_stress.dat" if self.plane_strain: filename = "material_check_stress_plane_strain.dat" self.initModel(akantu.SolidMechanicsModelOptions(akantu._static), filename) solver = self.model.getNonLinearSolver() solver.set("max_iterations", 2) solver.set("threshold", 2e-4) solver.set("convergence_type", akantu._scc_residual) self.model.solveStep() self.checkAll() -akantu.initialize() + TestPatchTestSMMLinear.TYPED_TEST(foo, "Static") diff --git a/test/test_model/patch_tests/patch_test_linear_fixture.py b/test/test_model/patch_tests/patch_test_linear_fixture.py index 1d0bd0d8e..d7195686d 100644 --- a/test/test_model/patch_tests/patch_test_linear_fixture.py +++ b/test/test_model/patch_tests/patch_test_linear_fixture.py @@ -1,141 +1,141 @@ #!/usr/bin/env python3 # ------------------------------------------------------------------------------ __author__ = "Guillaume Anciaux" __copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ " en Mécanique des Solides)" __credits__ = ["Guillaume Anciaux"] __license__ = "L-GPLv3" __maintainer__ = "Guillaume Anciaux" __email__ = "guillaume.anciaux@epfl.ch" # ------------------------------------------------------------------------------ import akantu import unittest import numpy as np import sys class TestPatchTestLinear(unittest.TestCase): alpha = np.array([[0.01, 0.02, 0.03, 0.04], [0.05, 0.06, 0.07, 0.08], [0.09, 0.10, 0.11, 0.12]]) gradient_tolerance = 1e-13 result_tolerance = 1e-13 dofs_tolerance = 1e-15 def __init__(self, test_name, elem_type_str, functor=None): self.test_name = test_name self.functor = functor self.elem_type = akantu.getElementTypes()[elem_type_str] self.elem_type_str = elem_type_str self.dim = akantu.Mesh.getSpatialDimension(self.elem_type) super().__init__(test_name) def _internal_call(self): self.functor(self) def __getattr__(self, key): if key == self.test_name: return self._internal_call def setUp(self): self.mesh = akantu.Mesh(self.dim, self.elem_type_str) self.mesh.read(str(self.elem_type_str) + ".msh") akantu.MeshUtils.buildFacets(self.mesh) self.mesh.createBoundaryGroupFromGeometry() self.model = self.model_type(self.mesh) def tearDown(self): del self.model del self.mesh def initModel(self, method, material_file): - akantu.getStaticParser().parse(material_file) + akantu.parseInput(material_file) akantu.setDebugLevel(akantu.dblError) self.model.initFull(method) self.applyBC() if method != akantu._static: self.model.setTimeStep(0.8 * self.model.getStableTimeStep()) def applyBC(self): boundary = self.model.getBlockedDOFs() for name, eg in self.mesh.getElementGroups().items(): nodes = eg['node_group'] boundary[nodes, :] = True def applyBConDOFs(self, dofs): coordinates = self.mesh.getNodes() for name, eg in self.mesh.getElementGroups().items(): nodes = eg['node_group'].flatten() dofs[nodes] = self.setLinearDOF(dofs[nodes], coordinates[nodes]) def prescribed_gradient(self, dof): gradient = self.alpha[:dof.shape[1], 1:self.dim + 1] return gradient def checkGradient(self, gradient, dofs): pgrad = self.prescribed_gradient(dofs).T gradient = gradient.reshape(gradient.shape[0], *pgrad.shape) diff = gradient[:] - pgrad norm = np.abs(pgrad).max() gradient_error = np.abs(diff).max() / norm self.assertAlmostEqual(0, gradient_error, delta=self.gradient_tolerance) def checkResults(self, presult_func, results, dofs): presult = presult_func(self.prescribed_gradient(dofs)).flatten() remaining_size = np.prod(np.array(results.shape[1:])) results = results.reshape((results.shape[0], remaining_size)) for result in results: diff = result - presult norm = np.abs(result).max() if norm == 0: result_error = np.abs(diff).max() else: result_error = np.abs(diff).max() / norm self.assertAlmostEqual(0., result_error, delta=self.result_tolerance) def setLinearDOF(self, dof, coord): nb_dofs = dof.shape[1] dof[:] = np.einsum('ik,ak->ai', self.alpha[:nb_dofs, 1:self.dim + 1], coord) for i in range(0, nb_dofs): dof[:, i] += self.alpha[i, 0] return dof def checkDOFs(self, dofs): coordinates = self.mesh.getNodes() ref_dofs = np.zeros_like(dofs) self.setLinearDOF(ref_dofs, coordinates) diff = dofs - ref_dofs dofs_error = np.abs(diff).max() self.assertAlmostEqual(0., dofs_error, delta=self.dofs_tolerance) @classmethod def TYPED_TEST(cls, functor, label): suite = unittest.TestSuite() for type_name, _type in akantu.getElementTypes().items(): if type_name == "_point_1": continue method_name = type_name + '_' + label test_case = cls(method_name, type_name, functor) suite.addTest(test_case) result = unittest.TextTestRunner(verbosity=1).run(suite) ret = not result.wasSuccessful() sys.exit(ret) diff --git a/test/test_model/patch_tests/patch_test_linear_heat_transfer_explicit.py b/test/test_model/patch_tests/patch_test_linear_heat_transfer_explicit.py index c827fe6d8..999fc03d7 100644 --- a/test/test_model/patch_tests/patch_test_linear_heat_transfer_explicit.py +++ b/test/test_model/patch_tests/patch_test_linear_heat_transfer_explicit.py @@ -1,36 +1,35 @@ #!/usr/bin/env python3 # ------------------------------------------------------------------------------ __author__ = "Guillaume Anciaux" __copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ " en Mécanique des Solides)" __credits__ = ["Guillaume Anciaux"] __license__ = "L-GPLv3" __maintainer__ = "Guillaume Anciaux" __email__ = "guillaume.anciaux@epfl.ch" # ------------------------------------------------------------------------------ from patch_test_linear_heat_transfer_fixture import TestPatchTestHTMLinear import akantu def foo(self): self.initModel( akantu.HeatTransferModelOptions(akantu._explicit_lumped_mass), "heat_transfer_input.dat") coordinates = self.mesh.getNodes() temperature = self.model.getTemperature() # set the position of all nodes to the static solution self.setLinearDOF(temperature, coordinates) for s in range(0, 100): self.model.solveStep() self.checkAll() -akantu.initialize() TestPatchTestHTMLinear.TYPED_TEST(foo, "Explicit") diff --git a/test/test_model/patch_tests/patch_test_linear_heat_transfer_implicit.py b/test/test_model/patch_tests/patch_test_linear_heat_transfer_implicit.py index 25674e2a3..24f070737 100644 --- a/test/test_model/patch_tests/patch_test_linear_heat_transfer_implicit.py +++ b/test/test_model/patch_tests/patch_test_linear_heat_transfer_implicit.py @@ -1,34 +1,33 @@ #!/usr/bin/env python3 # ------------------------------------------------------------------------------ __author__ = "Guillaume Anciaux" __copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ " en Mécanique des Solides)" __credits__ = ["Guillaume Anciaux"] __license__ = "L-GPLv3" __maintainer__ = "Guillaume Anciaux" __email__ = "guillaume.anciaux@epfl.ch" # ------------------------------------------------------------------------------ from patch_test_linear_heat_transfer_fixture import TestPatchTestHTMLinear import akantu def foo(self): self.initModel(akantu.HeatTransferModelOptions(akantu._implicit_dynamic), "heat_transfer_input.dat") coordinates = self.mesh.getNodes() temperature = self.model.getTemperature() # set the position of all nodes to the static solution self.setLinearDOF(temperature, coordinates) for s in range(0, 100): self.model.solveStep() self.checkAll() -akantu.initialize() TestPatchTestHTMLinear.TYPED_TEST(foo, "Explicit") diff --git a/test/test_model/patch_tests/patch_test_linear_heat_transfer_static.py b/test/test_model/patch_tests/patch_test_linear_heat_transfer_static.py index ede63a7b8..58e009d54 100644 --- a/test/test_model/patch_tests/patch_test_linear_heat_transfer_static.py +++ b/test/test_model/patch_tests/patch_test_linear_heat_transfer_static.py @@ -1,33 +1,32 @@ #!/usr/bin/env python3 # ------------------------------------------------------------------------------ __author__ = "Guillaume Anciaux" __copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ " en Mécanique des Solides)" __credits__ = ["Guillaume Anciaux"] __license__ = "L-GPLv3" __maintainer__ = "Guillaume Anciaux" __email__ = "guillaume.anciaux@epfl.ch" # ------------------------------------------------------------------------------ from patch_test_linear_heat_transfer_fixture import TestPatchTestHTMLinear import akantu def foo(self): self.initModel(akantu.HeatTransferModelOptions(akantu._static), "heat_transfer_input.dat") solver = self.model.getNonLinearSolver() solver.set("max_iterations", 2) solver.set("threshold", 2e-4) solver.set("convergence_type", akantu._scc_residual) self.model.solveStep() self.checkAll() -akantu.initialize() TestPatchTestHTMLinear.TYPED_TEST(foo, "Static") diff --git a/test/test_python_interface/test_boundary_condition_functors.py b/test/test_python_interface/test_boundary_condition_functors.py index 6190931b8..dc2d6816d 100644 --- a/test/test_python_interface/test_boundary_condition_functors.py +++ b/test/test_python_interface/test_boundary_condition_functors.py @@ -1,83 +1,82 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # ------------------------------------------------------------------------------ __author__ = "Lucas Frérot" __copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ " en Mécanique des Solides)" __credits__ = ["Lucas Frérot"] __license__ = "L-GPLv3" __maintainer__ = "Lucas Frérot" __email__ = "lucas.frerot@epfl.ch" # ------------------------------------------------------------------------------ import sys import os import numpy as np import akantu as aka ###################################################################### # Boundary conditions founctors ###################################################################### class FixedValue: """Functor for Dirichlet boundary conditions""" def __init__(self, value, axis): self.value = value axis_dict = {'x':0, 'y':1, 'z':2} self.axis = axis_dict[axis] def operator(self, node, flags, primal, coord): primal[self.axis] = self.value flags[self.axis] = True #--------------------------------------------------------------------- class FromStress: """Functor for Neumann boundary conditions""" def __init__(self, stress): self.stress = stress def operator(self, quad_point, dual, coord, normals): dual[:] = np.dot(self.stress,normals) ###################################################################### def main(): - aka.initialize("input_test.dat") + aka.parseInput("input_test.dat") mesh = aka.Mesh(2) mesh.read('mesh_dcb_2d.msh') model = aka.SolidMechanicsModel(mesh, 2) model.initFull() model.applyDirichletBC(FixedValue(0.0, 'x'), "edge") stress = np.array([[1, 0], [0, 0]]) blocked_nodes = mesh.getElementGroup("edge").getNodes().flatten() boundary = model.getBlockedDOFs() # Testing that nodes are correctly blocked for n in blocked_nodes: if not boundary[n, 0]: return -1 boundary.fill(False) model.applyNeumannBC(FromStress(stress), "edge") force = model.getForce() # Checking that nodes have a force in the correct direction for n in blocked_nodes: if not force[n, 0] > 0: return -1 - aka.finalize() return 0 if __name__ == "__main__": sys.exit(main()) diff --git a/test/test_python_interface/test_mesh_interface.py b/test/test_python_interface/test_mesh_interface.py index d48fc08a2..2ec7e8358 100644 --- a/test/test_python_interface/test_mesh_interface.py +++ b/test/test_python_interface/test_mesh_interface.py @@ -1,37 +1,34 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # ------------------------------------------------------------------------------ __author__ = "Lucas Frérot" __copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ " en Mécanique des Solides)" __credits__ = ["Lucas Frérot"] __license__ = "L-GPLv3" __maintainer__ = "Lucas Frérot" __email__ = "lucas.frerot@epfl.ch" # ------------------------------------------------------------------------------ import sys import os import akantu as aka def main(): - aka.initialize() - mesh = aka.Mesh(2) mesh.read('mesh_dcb_2d.msh') # Tests the getNbElement() function if mesh.getNbElement(aka._quadrangle_8) \ != mesh.getNbElementByDimension(2): raise Exception("Number of elements wrong, should be"\ " {}".format(mesh.getNbElementByDimension(2))) return -1 # TODO test the other functions in Mesh - aka.finalize() return 0 if __name__ == "__main__": sys.exit(main()) diff --git a/test/test_python_interface/test_multiple_init.py b/test/test_python_interface/test_multiple_init.py index d70e99934..8876d0e96 100755 --- a/test/test_python_interface/test_multiple_init.py +++ b/test/test_python_interface/test_multiple_init.py @@ -1,37 +1,37 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # ------------------------------------------------------------------------------ __author__ = "Fabian Barras" __copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ " en Mécanique des Solides)" __credits__ = ["Fabian Barras"] __license__ = "L-GPLv3" __maintainer__ = "Fabian Barras" __email__ = "fabian.barras@epfl.ch" # ------------------------------------------------------------------------------ import sys import os import akantu as aka -aka.initialize('input_test.dat') +aka.parseInput('input_test.dat') print('First initialisation') mesh = aka.Mesh(2) mesh.read('mesh_dcb_2d.msh') model = aka.SolidMechanicsModel(mesh) model.initFull(aka.SolidMechanicsModelOptions(aka._static)) del model del mesh print('Second initialisation') mesh = aka.Mesh(2) mesh.read('mesh_dcb_2d.msh') model = aka.SolidMechanicsModel(mesh) model.initFull(aka.SolidMechanicsModelOptions(aka._static)) del model del mesh print('All right')