diff --git a/examples/python/custom-material/bi-material.py b/examples/python/custom-material/bi-material.py index 096e020f3..8fa6956ee 100644 --- a/examples/python/custom-material/bi-material.py +++ b/examples/python/custom-material/bi-material.py @@ -1,188 +1,193 @@ from __future__ import print_function # ------------------------------------------------------------- # import akantu as aka import subprocess import numpy as np import time import os # ------------------------------------------------------------- # 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 parameters spatial_dimension = 2 mesh_file = 'square.msh' if not os.path.isfile(mesh_file): # call gmsh to generate the mesh ret = subprocess.call( 'gmsh -format msh2 -2 square.geo -optimize square.msh', shell=True) if ret != 0: raise Exception( 'execution of GMSH failed: do you have it installed ?') time.sleep(2) # read mesh mesh = aka.Mesh(spatial_dimension) mesh.read(mesh_file) # create the custom material mat = LocalElastic() -aka.registerNewPythonMaterial(mat, "local_elastic") + +mat_factory = aka.MaterialFactory.getInstance() + +mat_factory.registerAllocator( + "local_elastic", + lambda a, _id__, model, _id: LocalElastic(model, _id)) # parse input file aka.parseInput('material.dat') # init the SolidMechanicsModel 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", 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() diff --git a/python/pybind11/py_aka_solid_mechanics_model.cc b/python/pybind11/py_aka_solid_mechanics_model.cc index 6da0b2e88..b05cc926d 100644 --- a/python/pybind11/py_aka_solid_mechanics_model.cc +++ b/python/pybind11/py_aka_solid_mechanics_model.cc @@ -1,122 +1,163 @@ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "non_linear_solver.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace py = pybind11; + +/* -------------------------------------------------------------------------- */ + +namespace akantu { + +class MaterialPython : public Material { + +public: + MaterialPython(SolidMechanicsModel & model, PyObject * obj, + const ID & id = ""); + + ~MaterialPython() override = default; + +public: + void registerInternals(); + + void initMaterial() override; + + void computeStress(ElementType el_type, + GhostType ghost_type = _not_ghost) override; + + void computeTangentModuli(const ElementType & el_type, + Array & tangent_matrix, + GhostType ghost_type = _not_ghost) override; + + Real getPushWaveSpeed(const Element & element) const override; + + Real getEnergy(const std::string & type) override; + + Real getEnergyForType(const std::string & type, ElementType el_type); + +protected: + std::map local_params; + std::map>> internals; +}; + +} // namespace akantu + +/* -------------------------------------------------------------------------- */ namespace _aka = akantu; namespace { /* -------------------------------------------------------------------------- */ #define def_deprecated(func_name, mesg) \ def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); }) #define def_function_nocopy(func_name) \ def(#func_name, \ [](_aka::SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](_aka::SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }) /* -------------------------------------------------------------------------- */ py::module & register_solid_mechanics_models(py::module & mod) { py::class_<_aka::SolidMechanicsModelOptions>( mod, "_aka::SolidMechanicsModelOptions") .def(py::init<_aka::AnalysisMethod>(), py::arg("analysis_method") = _aka::_explicit_lumped_mass); py::class_<_aka::Model>(mod, "Model") .def("setBaseName", &_aka::Model::setBaseName) .def("addDumpFieldVector", &_aka::Model::addDumpFieldVector) .def("addDumpField", &_aka::Model::addDumpField) .def("dump", &_aka::Model::dump); py::class_<_aka::NonLinearSolver>(mod, "NonLinearSolver") .def("set", [](_aka::NonLinearSolver & self, const std::string & id, const _aka::Real & val) { if (id == "max_iterations") self.set(id, int(val)); else if (id == "convergence_type") self.set(id, akantu::SolveConvergenceCriteria(_aka::UInt(val))); else self.set(id, val); }) .def("set", &_aka::NonLinearSolver::set<_aka::SolveConvergenceCriteria>); py::class_<_aka::ModelSolver>(mod, "ModelSolver") .def("getNonLinearSolver", (_aka::NonLinearSolver & (_aka::ModelSolver::*)(const _aka::ID &)) & _aka::ModelSolver::getNonLinearSolver, py::arg("solver_id") = "", py::return_value_policy::reference) .def("solveStep", &_aka::ModelSolver::solveStep, py::arg("solver_id") = ""); py::class_<_aka::SolidMechanicsModel, _aka::Model, _aka::ModelSolver>( mod, "_aka::SolidMechanicsModel") .def(py::init<_aka::Mesh &, _aka::UInt, const _aka::ID &, const _aka::MemoryID &, const _aka::ModelType>(), py::arg("mesh"), py::arg("spatial_dimension") = _aka::_all_dimensions, py::arg("id") = "solid_mechanics_model", py::arg("memory_id") = 0, py::arg("model_type") = _aka::ModelType::_solid_mechanics_model) .def("initFull", [](_aka::SolidMechanicsModel & self, const _aka::SolidMechanicsModelOptions & options) { self.initFull(options); }, py::arg("_analysis_method") = _aka::SolidMechanicsModelOptions()) .def("initFull", [](_aka::SolidMechanicsModel & self, const _aka::AnalysisMethod & _analysis_method) { self.initFull(_aka::SolidMechanicsModelOptions(_analysis_method)); }, py::arg("_analysis_method")) .def_deprecated("applyDirichletBC", "Deprecated: use applyBC") .def("applyBC", [](_aka::SolidMechanicsModel & self, _aka::BC::DirichletFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("applyBC", [](_aka::SolidMechanicsModel & self, _aka::BC::NeumannFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("setTimeStep", &_aka::SolidMechanicsModel::setTimeStep, py::arg("time_step"), py::arg("solver_id") = "") .def("getEnergy", py::overload_cast( &_aka::SolidMechanicsModel::getEnergy), py::arg("energy_id")) .def_function(assembleStiffnessMatrix) .def_function(assembleInternalForces) .def_function(getStableTimeStep) .def_function_nocopy(getExternalForce) .def_function_nocopy(getDisplacement) .def_function_nocopy(getPreviousDisplacement) .def_function_nocopy(getIncrement) .def_function_nocopy(getMass) .def_function_nocopy(getVelocity) .def_function_nocopy(getAcceleration) .def_function_nocopy(getInternalForce) .def_function_nocopy(getBlockedDOFs) .def_function_nocopy(getIncrementFlag); + py::class_<_aka::MaterialFactory>(mod, "MaterialFactory"); + return mod; } // namespace } // namespace