diff --git a/examples/python/cohesive/custom_material.py b/examples/python/cohesive/custom_material.py new file mode 100755 index 000000000..509c7d10f --- /dev/null +++ b/examples/python/cohesive/custom_material.py @@ -0,0 +1,181 @@ +#!/usr/bin/env python3 +# pylint: disable=missing-module-docstring +# pylint: disable=missing-function-docstring + +import numpy as np +import akantu as aka + +spatial_dimension = 2 + +# ------------------------------------------------------------------------------ +class LinearCohesive(aka.MaterialCohesive): + """Material Cohesive Linear.""" + + def __init__(self, model, _id): + super().__init__(model, _id) + super().registerParamReal( + "G_c", aka._pat_readable | aka._pat_parsable, "Fracture energy" + ) + super().registerParamReal("beta", aka._pat_readable | aka._pat_parsable, "beta") + self.registerInternalReal("delta_max", 1) + self.beta = 0 + self.sigma_c = 0 + self.delta_c = 0 + + def initMaterial(self): + super().initMaterial() + self.sigma_c = self.getReal("sigma_c") + self.G_c = self.getReal("G_c") + self.beta = self.getReal("beta") + self.delta_c = 2 * self.G_c / self.sigma_c + + def checkInsertion(self, check_only): + model = self.getModel() + facets = self.getFacetFilter() + inserter = model.getElementInserter() + + for type_facet in facets.elementTypes(dim=(spatial_dimension - 1)): + facet_filter = facets(type_facet) + nb_facet = facet_filter.shape[0] + if nb_facet == 0: + continue + + fe_engine = model.getFEEngine("FacetsFEEngine") + + facets_check = inserter.getCheckFacets(type_facet) + insertion = inserter.getInsertionFacets(type_facet) + + nb_quad_facet = fe_engine.getNbIntegrationPoints(type_facet) + normals = fe_engine.getNormalsOnIntegrationPoints(type_facet) + facets_stresses = model.getStressOnFacets(type_facet).reshape( + normals.shape[0] // nb_quad_facet, + nb_quad_facet, + 2, + spatial_dimension, + spatial_dimension, + ) + + tangents = model.getTangents(type_facet) + + for facet, facet_stresses in zip(facet_filter, facets_stresses): + if not facets_check[facet]: + continue + + ref_stress = 0 + + for q, quad_stresses in enumerate(facet_stresses): + current_quad = facet * nb_quad_facet + q + normal = normals[current_quad, :].ravel() + tangent = tangents[current_quad, :].ravel() + + stress_1 = quad_stresses.T[0] + stress_2 = quad_stresses.T[1] + + avg_stress = stress_1 + stress_2 / 2.0 + traction = avg_stress.dot(normal) + + n = traction.dot(normal) + n = max(0, n) + t = traction.dot(tangent) + + ref_stress = max( + ref_stress, np.sqrt(n * n + t * t / (self.beta**2)) + ) + + if ref_stress > self.sigma_c: + insertion[facet] = True + + # constitutive law + def computeTraction(self, normals, el_type, ghost_type): + openings = self.getOpening(el_type, ghost_type) + tractions = self.getTraction(el_type, ghost_type) + + delta_max = self.getInternalReal("delta_max")(el_type) + + for el in range(normals.shape[0]): + normal = normals[el].ravel() + opening = openings[el].ravel() + + delta_n = opening.dot(normal) * normal + delta_s = opening - delta_n + + delta = ( + self.beta * np.linalg.norm(delta_s) ** 2 + np.linalg.norm(delta_n) ** 2 + ) + + delta_max[el] = max(delta_max[el], delta) + + tractions[el, :] = ( + (delta * delta_s + delta_n) + * self.sigma_c + / delta + * (1 - delta / self.delta_c) + ) + + +# register material to the MaterialFactory +def allocator(_dim, unused, model, _id): + return LinearCohesive(model, _id) + + +mat_factory = aka.MaterialFactory.getInstance() +mat_factory.registerAllocator("local_cohesive", allocator) + +# ------------------------------------------------------------------------- +# Initialization +# ------------------------------------------------------------------------- +aka.parseInput("local_material.dat") +mesh = aka.Mesh(spatial_dimension) +mesh.read("plate.msh") + +model = aka.SolidMechanicsModelCohesive(mesh) +model.initFull(_analysis_method=aka._static, _is_extrinsic=True) +model.initNewSolver(aka._explicit_lumped_mass) + +model.setBaseName("plate") +model.addDumpFieldVector("displacement") +model.addDumpFieldVector("external_force") +model.addDumpField("strain") +model.addDumpField("stress") +model.addDumpField("blocked_dofs") + +model.setBaseNameToDumper("cohesive elements", "cohesive") +model.addDumpFieldVectorToDumper("cohesive elements", "displacement") +model.addDumpFieldToDumper("cohesive elements", "damage") +model.addDumpFieldVectorToDumper("cohesive elements", "tractions") +model.addDumpFieldVectorToDumper("cohesive elements", "opening") + +# ------------------------------------------------------------------------- +# Boundary conditions +# ------------------------------------------------------------------------- +model.applyBC(aka.FixedValue(0.0, aka._x), "XBlocked") +model.applyBC(aka.FixedValue(0.0, aka._y), "YBlocked") + +trac = np.zeros(spatial_dimension) +traction = 0.095 +trac[int(aka._y)] = traction + +model.getExternalForce()[:] = 0 +model.applyBC(aka.FromTraction(trac), "Traction") + +print("Solve for traction ", traction) +solver = model.getNonLinearSolver("static") +solver.set("max_iterations", 100) +solver.set("threshold", 1e-10) +solver.set("convergence_type", aka.SolveConvergenceCriteria.residual) + +model.solveStep("static") +model.dump() +model.dump("cohesive elements") + +model.setTimeStep(model.getStableTimeStep() * 0.1) + +maxsteps = 100 + +for i in range(0, maxsteps): + print("{0}/{1}".format(i, maxsteps)) + model.checkCohesiveStress() + model.solveStep("explicit_lumped") + if i % 10 == 0: + model.dump() + model.dump("cohesive elements") diff --git a/examples/python/cohesive/local_material.dat b/examples/python/cohesive/local_material.dat new file mode 100644 index 000000000..b0ee626b6 --- /dev/null +++ b/examples/python/cohesive/local_material.dat @@ -0,0 +1,20 @@ +model solid_mechanics_model_cohesive [ + cohesive_inserter [ + bounding_box = [[0,10],[-10, 10]] + ] + + material elastic [ + name = virtual + rho = 1 # density + E = 1 # young's modulus + nu = 0.3 # poisson's ratio + finite_deformation = true + ] + + material local_cohesive [ + name = cohesive + sigma_c = 0.1 + G_c = 1e-2 + beta = 1. + ] +] diff --git a/examples/python/cohesive/plate.py b/examples/python/cohesive/plate.py index fccde62da..5f223c289 100644 --- a/examples/python/cohesive/plate.py +++ b/examples/python/cohesive/plate.py @@ -1,96 +1,97 @@ #!/usr/bin/env python3 """plate.py: Python example: plate with a hole breaking with cohesive elements""" __author__ = "Guillaume Anciaux" __credits__ = [ "Guillaume Anciaux ", ] -__copyright__ = "Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale" \ - " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ - " en Mécanique des Solides)" +__copyright__ = ( + "Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale" + " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" + " en Mécanique des Solides)" +) __license__ = "LGPLv3" import akantu as aka import numpy as np def solve(material_file, mesh_file, traction): aka.parseInput(material_file) spatial_dimension = 2 # ------------------------------------------------------------------------- # Initialization # ------------------------------------------------------------------------- mesh = aka.Mesh(spatial_dimension) mesh.read(mesh_file) model = aka.SolidMechanicsModelCohesive(mesh) - model.initFull(_analysis_method=aka._static, - _is_extrinsic=True) + model.initFull(_analysis_method=aka._static, _is_extrinsic=True) model.initNewSolver(aka._explicit_lumped_mass) - model.setBaseName('plate') - model.addDumpFieldVector('displacement') - model.addDumpFieldVector('external_force') - model.addDumpField('strain') - model.addDumpField('stress') - model.addDumpField('blocked_dofs') + model.setBaseName("plate") + model.addDumpFieldVector("displacement") + model.addDumpFieldVector("external_force") + model.addDumpField("strain") + model.addDumpField("stress") + model.addDumpField("blocked_dofs") - model.setBaseNameToDumper('cohesive elements', 'cohesive') - model.addDumpFieldVectorToDumper('cohesive elements', 'displacement') - model.addDumpFieldToDumper('cohesive elements', 'damage') - model.addDumpFieldVectorToDumper('cohesive elements', 'tractions') - model.addDumpFieldVectorToDumper('cohesive elements', 'opening') + model.setBaseNameToDumper("cohesive elements", "cohesive") + model.addDumpFieldVectorToDumper("cohesive elements", "displacement") + model.addDumpFieldToDumper("cohesive elements", "damage") + model.addDumpFieldVectorToDumper("cohesive elements", "tractions") + model.addDumpFieldVectorToDumper("cohesive elements", "opening") # ------------------------------------------------------------------------- # Boundary conditions # ------------------------------------------------------------------------- - model.applyBC(aka.FixedValue(0.0, aka._x), 'XBlocked') - model.applyBC(aka.FixedValue(0.0, aka._y), 'YBlocked') + model.applyBC(aka.FixedValue(0.0, aka._x), "XBlocked") + model.applyBC(aka.FixedValue(0.0, aka._y), "YBlocked") trac = np.zeros(spatial_dimension) trac[int(aka._y)] = traction - print('Solve for traction ', traction) + print("Solve for traction ", traction) model.getExternalForce()[:] = 0 - model.applyBC(aka.FromTraction(trac), 'Traction') + model.applyBC(aka.FromTraction(trac), "Traction") - solver = model.getNonLinearSolver('static') - solver.set('max_iterations', 100) - solver.set('threshold', 1e-10) + solver = model.getNonLinearSolver("static") + solver.set("max_iterations", 100) + solver.set("threshold", 1e-10) solver.set("convergence_type", aka.SolveConvergenceCriteria.residual) - model.solveStep('static') + model.solveStep("static") model.dump() - model.dump('cohesive elements') + model.dump("cohesive elements") - model.setTimeStep(model.getStableTimeStep()*0.1) + model.setTimeStep(model.getStableTimeStep() * 0.1) maxsteps = 100 for i in range(0, maxsteps): - print('{0}/{1}'.format(i, maxsteps)) + print("{0}/{1}".format(i, maxsteps)) model.checkCohesiveStress() - model.solveStep('explicit_lumped') + model.solveStep("explicit_lumped") if i % 10 == 0: model.dump() - model.dump('cohesive elements') + model.dump("cohesive elements") # ----------------------------------------------------------------------------- # main # ----------------------------------------------------------------------------- def main(): - mesh_file = 'plate.msh' - material_file = 'material.dat' + mesh_file = "plate.msh" + material_file = "material.dat" - traction = .095 + traction = 0.095 solve(material_file, mesh_file, traction) # ----------------------------------------------------------------------------- -if __name__ == '__main__': +if __name__ == "__main__": main() diff --git a/examples/python/custom-material/custom-material.py b/examples/python/custom-material/custom-material.py index a10a120e0..fae859398 100644 --- a/examples/python/custom-material/custom-material.py +++ b/examples/python/custom-material/custom-material.py @@ -1,189 +1,188 @@ #!/usr/bin/env python3 """ custom-material.py: Custom material example""" __author__ = "Guillaume Anciaux" __credits__ = [ "Guillaume Anciaux ", ] -__copyright__ = "Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale" \ - " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ - " en Mécanique des Solides)" +__copyright__ = ( + "Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale" + " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" + " en Mécanique des Solides)" +) __license__ = "LGPLv3" import numpy as np import akantu as aka # ------------------------------------------------------------------------------ class LocalElastic(aka.Material): - def __init__(self, model, _id): super().__init__(model, _id) - super().registerParamReal('E', - aka._pat_readable | aka._pat_parsable, - 'Youngs modulus') - super().registerParamReal('nu', - aka._pat_readable | aka._pat_parsable, - 'Poisson ratio') + super().registerParamReal( + "E", aka._pat_readable | aka._pat_parsable, "Youngs modulus" + ) + super().registerParamReal( + "nu", aka._pat_readable | aka._pat_parsable, "Poisson ratio" + ) def initMaterial(self): - nu = self.getReal('nu') - E = self.getReal('E') + nu = self.getReal("nu") + E = self.getReal("E") self.mu = E / (2 * (1 + nu)) - self.lame_lambda = nu * E / ( - (1. + nu) * (1. - 2. * nu)) + self.lame_lambda = nu * E / ((1.0 + nu) * (1.0 - 2.0 * nu)) # Second Lame coefficient (shear modulus) - self.lame_mu = E / (2. * (1. + nu)) + self.lame_mu = E / (2.0 * (1.0 + nu)) super().initMaterial() # declares all the parameters that are needed def getPushWaveSpeed(self, element): - rho = self.getReal('rho') + rho = self.getReal("rho") return np.sqrt((self.lame_lambda + 2 * self.lame_mu) / rho) # compute small deformation tensor @staticmethod def computeEpsilon(grad_u): - return 0.5 * (grad_u + np.einsum('aij->aji', grad_u)) + return 0.5 * (grad_u + np.einsum("aij->aji", grad_u)) # constitutive law def computeStress(self, el_type, ghost_type): grad_u = self.getGradU(el_type, ghost_type) sigma = self.getStress(el_type, ghost_type) n_quads = grad_u.shape[0] grad_u = grad_u.reshape((n_quads, 2, 2)) epsilon = self.computeEpsilon(grad_u) sigma = sigma.reshape((n_quads, 2, 2)) - trace = np.einsum('aii->a', grad_u) + 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) + np.einsum("a,ij->aij", trace, self.lame_lambda * np.eye(2)) + + 2.0 * self.lame_mu * epsilon + ) # constitutive law tangent modulii def computeTangentModuli(self, el_type, tangent_matrix, ghost_type): n_quads = tangent_matrix.shape[0] tangent = tangent_matrix.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 computePotentialEnergy(self, el_type): sigma = self.getStress(el_type) grad_u = self.getGradU(el_type) nquads = sigma.shape[0] stress = sigma.reshape(nquads, 2, 2) grad_u = grad_u.reshape((nquads, 2, 2)) epsilon = self.computeEpsilon(grad_u) energy_density = self.getPotentialEnergy(el_type) - energy_density[:, 0] = 0.5 * np.einsum('aij,aij->a', stress, epsilon) + energy_density[:, 0] = 0.5 * np.einsum("aij,aij->a", stress, epsilon) # register material to the MaterialFactory def allocator(_dim, unused, model, _id): return LocalElastic(model, _id) mat_factory = aka.MaterialFactory.getInstance() mat_factory.registerAllocator("local_elastic", allocator) # ------------------------------------------------------------------------------ # main # ------------------------------------------------------------------------------ spatial_dimension = 2 -aka.parseInput('material.dat') +aka.parseInput("material.dat") -mesh_file = 'bar.msh' +mesh_file = "bar.msh" max_steps = 250 time_step = 1e-3 # ------------------------------------------------------------------------------ # Initialization # ------------------------------------------------------------------------------ mesh = aka.Mesh(spatial_dimension) mesh.read(mesh_file) # parse input file -aka.parseInput('material.dat') +aka.parseInput("material.dat") model = aka.SolidMechanicsModel(mesh) model.initFull(_analysis_method=aka._explicit_lumped_mass) 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.applyBC(aka.FixedValue(0, aka._x), "XBlocked") model.applyBC(aka.FixedValue(0, aka._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. + x = position[i, 0] - 5.0 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)) + 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 # ------------------------------------------------------------------------------ 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') + 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))) + print( + "{0},{1},{2},{3},{4}".format(step, step * time_step, epot, ekin, (epot + ekin)) + ) diff --git a/python/py_fe_engine.cc b/python/py_fe_engine.cc index baa7e5225..11af53f0a 100644 --- a/python/py_fe_engine.cc +++ b/python/py_fe_engine.cc @@ -1,155 +1,160 @@ /** * @file py_fe_engine.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Wed Nov 27 2019 * @date last modification: Sat Dec 12 2020 * * @brief pybind11 interface to FEEngine * * * @section LICENSE * * Copyright (©) 2018-2021 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 "py_aka_array.hh" #include "py_aka_common.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { void register_fe_engine(py::module & mod) { py::class_(mod, "Element") .def(py::init([](ElementType type, UInt id) { return new Element{type, id, _not_ghost}; })) .def(py::init([](ElementType type, UInt id, GhostType ghost_type) { return new Element{type, id, ghost_type}; })) .def("__lt__", [](Element & self, const Element & other) { return (self < other); }) .def("__repr__", [](Element & self) { return std::to_string(self); }); mod.attr("ElementNull") = ElementNull; py::class_(mod, "FEEngine") .def( "getNbIntegrationPoints", [](FEEngine & fem, const ElementType & type, const GhostType & ghost_type) { return fem.getNbIntegrationPoints(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "gradientOnIntegrationPoints", [](FEEngine & fem, const Array & u, Array & nablauq, UInt nb_degree_of_freedom, ElementType type, GhostType ghost_type, const Array * filter_elements) { if (filter_elements == nullptr) { // This is due to the ArrayProxy that looses the // empty_filter information filter_elements = &empty_filter; } fem.gradientOnIntegrationPoints(u, nablauq, nb_degree_of_freedom, type, ghost_type, *filter_elements); }, py::arg("u"), py::arg("nablauq"), py::arg("nb_degree_of_freedom"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::arg("filter_elements") = nullptr) .def( "interpolateOnIntegrationPoints", [](FEEngine & self, const Array & u, Array & uq, UInt nb_degree_of_freedom, ElementType type, GhostType ghost_type, const Array * filter_elements) { if (filter_elements == nullptr) { // This is due to the ArrayProxy that looses the // empty_filter information filter_elements = &empty_filter; } self.interpolateOnIntegrationPoints(u, uq, nb_degree_of_freedom, type, ghost_type, *filter_elements); }, py::arg("u"), py::arg("uq"), py::arg("nb_degree_of_freedom"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::arg("filter_elements") = nullptr) .def( "interpolateOnIntegrationPoints", [](FEEngine & self, const Array & u, ElementTypeMapArray & uq, const ElementTypeMapArray * filter_elements) { self.interpolateOnIntegrationPoints(u, uq, filter_elements); }, py::arg("u"), py::arg("uq"), py::arg("filter_elements") = nullptr) .def( "computeIntegrationPointsCoordinates", [](FEEngine & self, ElementTypeMapArray & coordinates, const ElementTypeMapArray * filter_elements) -> decltype(auto) { return self.computeIntegrationPointsCoordinates(coordinates, filter_elements); }, py::arg("coordinates"), py::arg("filter_elements") = nullptr) .def( "assembleFieldLumped", [](FEEngine & fem, const std::function &, const Element &)> & field_funct, const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager, ElementType type, GhostType ghost_type) { fem.assembleFieldLumped(field_funct, matrix_id, dof_id, dof_manager, type, ghost_type); }, py::arg("field_funct"), py::arg("matrix_id"), py::arg("dof_id"), py::arg("dof_manager"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "assembleFieldMatrix", [](FEEngine & fem, const std::function &, const Element &)> & field_funct, const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager, ElementType type, GhostType ghost_type = _not_ghost) { fem.assembleFieldMatrix(field_funct, matrix_id, dof_id, dof_manager, type, ghost_type); }, py::arg("field_funct"), py::arg("matrix_id"), py::arg("dof_id"), py::arg("dof_manager"), py::arg("type"), py::arg("ghost_type") = _not_ghost) - .def("getElementInradius", [](FEEngine & self, const Element & element) { - return self.getElementInradius(element); - }); + .def("getElementInradius", + [](FEEngine & self, const Element & element) { + return self.getElementInradius(element); + }) + .def("getNormalsOnIntegrationPoints", + &FEEngine::getNormalsOnIntegrationPoints, py::arg("type"), + py::arg("ghost_type") = _not_ghost, + py::return_value_policy::reference); py::class_(mod, "IntegrationPoint"); } } // namespace akantu diff --git a/python/py_material.cc b/python/py_material.cc index 3164298dd..90a92c22b 100644 --- a/python/py_material.cc +++ b/python/py_material.cc @@ -1,261 +1,378 @@ /** * @file py_material.cc * * @author Guillaume Anciaux * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Thu Jun 20 2019 * @date last modification: Fri Apr 09 2021 * * @brief pybind11 interface to Material * * * @section LICENSE * * Copyright (©) 2018-2021 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 "py_aka_array.hh" #include "py_akantu_pybind11_compatibility.hh" /* -------------------------------------------------------------------------- */ #include #include #if defined(AKANTU_COHESIVE_ELEMENT) +#include #include #endif #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { namespace { template class PyMaterial : public _Material { public: /* Inherit the constructors */ using _Material::_Material; ~PyMaterial() override = default; void initMaterial() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, initMaterial, ); }; void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, _Material, computeStress, el_type, ghost_type); } void computeTangentModuli(ElementType el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, computeTangentModuli, el_type, tangent_matrix, ghost_type); } void computePotentialEnergy(ElementType el_type) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, computePotentialEnergy, el_type); } Real getPushWaveSpeed(const Element & element) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(Real, _Material, getPushWaveSpeed, element); } Real getShearWaveSpeed(const Element & element) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(Real, _Material, getShearWaveSpeed, element); } template void registerInternal(const std::string & name, UInt nb_component) { auto && internal = std::make_shared>(name, *this); AKANTU_DEBUG_INFO("alloc internal " << name << " " << &this->internals[name]); internal->initialize(nb_component); this->internals[name] = internal; } protected: std::map> internals; }; + /* ------------------------------------------------------------------------ */ + template + void register_material_classes(py::module & mod, const std::string & name) { + py::class_<_Material, Material, Parsable, PyMaterial<_Material>>( + mod, name.c_str(), py::multiple_inheritance()) + .def(py::init()); + } + +#if defined(AKANTU_COHESIVE_ELEMENT) + template class PyMaterialCohesive : public _Material { + public: + /* Inherit the constructors */ + using _Material::_Material; + + ~PyMaterialCohesive() override = default; + void initMaterial() override { + // NOLINTNEXTLINE + PYBIND11_OVERRIDE(void, _Material, initMaterial, ); + }; + + void computeStress(ElementType el_type, + GhostType ghost_type = _not_ghost) final {} + + void checkInsertion(bool check_only) override { + // NOLINTNEXTLINE + PYBIND11_OVERRIDE(void, _Material, checkInsertion, check_only); + } + + void computeTraction(const Array & normal, ElementType el_type, + GhostType ghost_type = _not_ghost) override { + // NOLINTNEXTLINE + PYBIND11_OVERRIDE_PURE(void, _Material, computeTraction, normal, el_type, + ghost_type); + } + + void computeTangentModuli(ElementType el_type, Array & tangent_matrix, + GhostType ghost_type = _not_ghost) final {} + + void computeTangentTraction(ElementType el_type, + Array & tangent_matrix, + const Array & normal, + GhostType ghost_type = _not_ghost) override { + // NOLINTNEXTLINE + PYBIND11_OVERRIDE(void, _Material, computeTangentTraction, el_type, + tangent_matrix, normal, ghost_type); + } + + template + void registerInternal(const std::string & name, UInt nb_component) { + auto && internal = + std::make_shared>(name, *this); + AKANTU_DEBUG_INFO("alloc internal " << name << " " + << &this->internals[name]); + + internal->initialize(nb_component); + this->internals[name] = internal; + } + + protected: + std::map> internals; + }; + + template + void register_material_cohesive_classes(py::module & mod, + const std::string & name) { + py::class_<_Material, Material, Parsable, PyMaterialCohesive<_Material>>( + mod, name.c_str(), py::multiple_inheritance()) + .def(py::init()); + } + +#endif + /* ------------------------------------------------------------------------ */ template void register_internal_field(py::module & mod, const std::string & name) { py::class_, ElementTypeMapArray, std::shared_ptr>>( mod, ("InternalField" + name).c_str()); } - /* ------------------------------------------------------------------------ */ - template - void register_material_classes(py::module & mod, const std::string & name) { - py::class_<_Material, Material, Parsable, PyMaterial<_Material>>( - mod, name.c_str(), py::multiple_inheritance()) - .def(py::init()); - } } // namespace /* -------------------------------------------------------------------------- */ void register_material(py::module & mod) { py::class_(mod, "MaterialFactory") .def_static( "getInstance", []() -> MaterialFactory & { return Material::getFactory(); }, py::return_value_policy::reference) .def("registerAllocator", [](MaterialFactory & self, const std::string id, py::function func) { self.registerAllocator( id, [func, id](UInt dim, const ID & /*unused*/, SolidMechanicsModel & model, const ID & option) -> std::unique_ptr { py::object obj = func(dim, id, model, option); auto & ptr = py::cast(obj); obj.release(); return std::unique_ptr(&ptr); }); }) .def("getPossibleAllocators", &MaterialFactory::getPossibleAllocators); register_internal_field(mod, "Real"); register_internal_field(mod, "UInt"); py::class_>( mod, "Material", py::multiple_inheritance()) .def(py::init()) .def( "getGradU", [](Material & self, ElementType el_type, GhostType ghost_type = _not_ghost) -> decltype(auto) { return self.getGradU(el_type, ghost_type); }, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "getStress", [](Material & self, ElementType el_type, GhostType ghost_type = _not_ghost) -> decltype(auto) { return self.getStress(el_type, ghost_type); }, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "getPotentialEnergy", [](Material & self, ElementType el_type) -> decltype(auto) { return self.getPotentialEnergy(el_type); }, py::arg("el_type"), py::return_value_policy::reference) .def( "getPotentialEnergy", [](Material & self, ElementType el_type, UInt index) -> Real { return self.getPotentialEnergy(el_type, index); }, py::arg("el_type"), py::arg("index")) .def("getPotentialEnergy", [](Material & self) -> Real { return self.getPotentialEnergy(); }) .def("initMaterial", &Material::initMaterial) .def("getModel", &Material::getModel) .def("registerInternalReal", [](Material & self, const std::string & name, UInt nb_component) { return dynamic_cast &>(self) .registerInternal(name, nb_component); }) .def("registerInternalUInt", [](Material & self, const std::string & name, UInt nb_component) { return dynamic_cast &>(self) .registerInternal(name, nb_component); }) .def( "getInternalReal", [](Material & self, const ID & id) -> decltype(auto) { return self.getInternal(id); }, py::arg("id"), py::return_value_policy::reference) .def( "getInternalUInt", [](Material & self, const ID & id) -> decltype(auto) { return self.getInternal(id); }, py::arg("id"), py::return_value_policy::reference) .def( "getElementFilter", [](Material & self) -> decltype(auto) { return self.getElementFilter(); }, py::return_value_policy::reference) /* * These functions override the `Parsable` interface. * This ensure that the `updateInternalParameters()` function is called. */ .def( "setReal", [](Material & self, const ID & name, const Real value) -> void { self.setParam(name, value); return; }, py::arg("name"), py::arg("value")) .def( "setBool", [](Material & self, const ID & name, const bool value) -> void { self.setParam(name, value); return; }, py::arg("name"), py::arg("value")) .def( "setString", [](Material & self, const ID & name, const std::string & value) -> void { self.setParam(name, value); return; }, py::arg("name"), py::arg("value")) .def( "setInt", [](Material & self, const ID & name, const int value) -> void { self.setParam(name, value); return; }, py::arg("name"), py::arg("value")) .def("getPushWaveSpeed", &Material::getPushWaveSpeed) .def("getShearWaveSpeed", &Material::getShearWaveSpeed) .def("__repr__", [](Material & self) { std::stringstream sstr; sstr << self; return sstr.str(); }); register_material_classes>(mod, "MaterialElastic2D"); register_material_classes>(mod, "MaterialElastic3D"); + +#if defined(AKANTU_COHESIVE_ELEMENT) + /* ------------------------------------------------------------------------ */ + py::class_>(mod, "MaterialCohesive", + py::multiple_inheritance()) + .def(py::init()) + .def("registerInternalReal", + [](MaterialCohesive & self, const std::string & name, + UInt nb_component) { + return dynamic_cast &>(self) + .registerInternal(name, nb_component); + }) + .def("registerInternalUInt", + [](MaterialCohesive & self, const std::string & name, + UInt nb_component) { + return dynamic_cast &>(self) + .registerInternal(name, nb_component); + }) + .def( + "getFacetFilter", + [](MaterialCohesive & self) -> decltype(auto) { + return self.getFacetFilter(); + }, + py::return_value_policy::reference) + .def( + "getFacetFilter", + [](MaterialCohesive & self, ElementType type, + GhostType ghost_type) -> decltype(auto) { + return self.getFacetFilter(type, ghost_type); + }, + py::arg("type"), py::arg("ghost_type") = _not_ghost, + py::return_value_policy::reference) + .def( + "getOpening", + [](MaterialCohesive & self, ElementType type, GhostType ghost_type) + -> decltype(auto) { return self.getOpening(type, ghost_type); }, + py::arg("type"), py::arg("ghost_type") = _not_ghost, + py::return_value_policy::reference) + .def( + "getTraction", + [](MaterialCohesive & self, ElementType type, GhostType ghost_type) + -> decltype(auto) { return self.getTraction(type, ghost_type); }, + py::arg("type"), py::arg("ghost_type") = _not_ghost, + py::return_value_policy::reference); + + register_material_cohesive_classes>( + mod, "MaterialCohesiveLinear2D"); + register_material_cohesive_classes>( + mod, "MaterialCohesiveLinear3D"); +#endif } } // namespace akantu diff --git a/python/py_solid_mechanics_model_cohesive.cc b/python/py_solid_mechanics_model_cohesive.cc index 2ce094bf2..81dce61dc 100644 --- a/python/py_solid_mechanics_model_cohesive.cc +++ b/python/py_solid_mechanics_model_cohesive.cc @@ -1,112 +1,133 @@ /** * @file py_solid_mechanics_model_cohesive.cc * * @author Nicolas Richart * * @date creation: Tue Jul 21 2020 * @date last modification: Tue Sep 29 2020 * * @brief pybind11 interface to SolidMechanicsModelCohesive * * * @section LICENSE * * Copyright (©) 2018-2021 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 "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ #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, \ [](SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }) void register_solid_mechanics_model_cohesive(py::module & mod) { py::class_(mod, "CohesiveElementInserter") .def("setLimit", &CohesiveElementInserter::setLimit) .def( "getCheckFacets", [](CohesiveElementInserter & self) -> decltype(auto) { return self.getCheckFacets(); }, py::return_value_policy::reference) .def( "getCheckFacets", [](CohesiveElementInserter & self, ElementType type, GhostType ghost_type) -> decltype(auto) { return self.getCheckFacets(type, ghost_type); }, + py::arg("type"), py::arg("ghost_type") = _not_ghost, + py::return_value_policy::reference) + .def( + "getInsertionFacets", + [](CohesiveElementInserter & self) -> decltype(auto) { + return self.getInsertionFacetsByElement(); + }, + py::return_value_policy::reference) + .def( + "getInsertionFacets", + [](CohesiveElementInserter & self, ElementType type, + GhostType ghost_type) -> decltype(auto) { + return self.getInsertionFacets(type, ghost_type); + }, + py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def("addPhysicalSurface", &CohesiveElementInserter::addPhysicalSurface) .def("addPhysicalVolume", &CohesiveElementInserter::addPhysicalVolume); py::class_( mod, "SolidMechanicsModelCohesiveOptions") .def(py::init(), py::arg("analysis_method") = _explicit_lumped_mass, py::arg("is_extrinsic") = false); py::class_( mod, "SolidMechanicsModelCohesive") .def(py::init(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "solid_mechanics_model") .def( "initFull", [](SolidMechanicsModel & self, const AnalysisMethod & analysis_method, bool is_extrinsic) { self.initFull(_analysis_method = analysis_method, _is_extrinsic = is_extrinsic); }, py::arg("_analysis_method") = _explicit_lumped_mass, py::arg("_is_extrinsic") = false) .def("checkCohesiveStress", &SolidMechanicsModelCohesive::checkCohesiveStress) .def("getElementInserter", &SolidMechanicsModelCohesive::getElementInserter, py::return_value_policy::reference) + .def("getStressOnFacets", &SolidMechanicsModelCohesive::getStressOnFacets, + py::arg("type"), py::arg("ghost_type") = _not_ghost, + py::return_value_policy::reference) + .def("getTangents", &SolidMechanicsModelCohesive::getTangents, + py::arg("type"), py::arg("ghost_type") = _not_ghost, + py::return_value_policy::reference) .def("updateAutomaticInsertion", &SolidMechanicsModelCohesive::updateAutomaticInsertion); } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh index 3a3dd4630..c3155d2d4 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh @@ -1,236 +1,236 @@ /** * @file material_cohesive.hh * * @author Nicolas Richart * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Jan 14 2021 * * @brief Specialization of the material class for cohesive elements * * * @section LICENSE * * Copyright (©) 2010-2021 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 "material.hh" /* -------------------------------------------------------------------------- */ #include "cohesive_internal_field.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_COHESIVE_HH_ #define AKANTU_MATERIAL_COHESIVE_HH_ /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModelCohesive; } namespace akantu { class MaterialCohesive : public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using MyFEEngineCohesiveType = FEEngineTemplate; public: MaterialCohesive(SolidMechanicsModel & model, const ID & id = ""); ~MaterialCohesive() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material computed parameter void initMaterial() override; /// compute tractions (including normals and openings) void computeTraction(GhostType ghost_type = _not_ghost); /// assemble residual void assembleInternalForces(GhostType ghost_type = _not_ghost) override; /// check stress for cohesive elements' insertion, by default it /// also updates the cohesive elements' data virtual void checkInsertion(bool /*check_only*/ = false) { AKANTU_TO_IMPLEMENT(); } /// interpolate stress on given positions for each element (empty /// implemantation to avoid the generic call to be done on cohesive elements) virtual void interpolateStress(const ElementType /*type*/, Array & /*result*/) {} /// compute the stresses - void computeAllStresses(GhostType /*ghost_type*/ = _not_ghost) override{}; + void computeAllStresses(GhostType /*ghost_type*/ = _not_ghost) final{}; // add the facet to be handled by the material UInt addFacet(const Element & element); protected: virtual void computeTangentTraction(ElementType /*el_type*/, Array & /*tangent_matrix*/, const Array & /*normal*/, GhostType /*ghost_type*/ = _not_ghost) { AKANTU_TO_IMPLEMENT(); } /// compute the normal void computeNormal(const Array & position, Array & normal, ElementType type, GhostType ghost_type); /// compute the opening void computeOpening(const Array & displacement, Array & opening, ElementType type, GhostType ghost_type); template void computeNormal(const Array & position, Array & normal, GhostType ghost_type); /// assemble stiffness void assembleStiffnessMatrix(GhostType ghost_type) override; /// constitutive law virtual void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// parallelism functions inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; protected: void updateEnergies(ElementType el_type) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the opening AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Opening, opening, Real); /// get the traction AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Traction, tractions, Real); /// get damage AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); /// get facet filter AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO(FacetFilter, facet_filter, const ElementTypeMapArray &); // AKANTU_GET_MACRO(ElementFilter, element_filter, const // ElementTypeMapArray &); /// compute reversible energy Real getReversibleEnergy(); /// compute dissipated energy Real getDissipatedEnergy(); /// compute contact energy Real getContactEnergy(); /// get energy Real getEnergy(const std::string & type) override; /// return the energy (identified by id) for the provided element Real getEnergy(const std::string & energy_id, ElementType type, UInt index) override { return Material::getEnergy(energy_id, type, index); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// list of facets assigned to this material ElementTypeMapArray facet_filter; /// Link to the cohesive fem object in the model FEEngine & fem_cohesive; private: /// reversible energy by quadrature point CohesiveInternalField reversible_energy; /// total energy by quadrature point CohesiveInternalField total_energy; protected: /// opening in all elements and quadrature points CohesiveInternalField opening; /// traction in all elements and quadrature points CohesiveInternalField tractions; /// traction due to contact CohesiveInternalField contact_tractions; /// normal openings for contact tractions CohesiveInternalField contact_opening; /// maximum displacement CohesiveInternalField delta_max; /// tell if the previous delta_max state is needed (in iterative schemes) bool use_previous_delta_max; /// tell if the previous opening state is needed (in iterative schemes) bool use_previous_opening; /// damage CohesiveInternalField damage; /// pointer to the solid mechanics model for cohesive elements SolidMechanicsModelCohesive * model; /// critical stress RandomInternalField sigma_c; /// critical displacement Real delta_c; /// array to temporarily store the normals Array normal; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "cohesive_internal_field_tmpl.hh" #include "material_cohesive_inline_impl.hh" #endif /* AKANTU_MATERIAL_COHESIVE_HH_ */