diff --git a/.gitlab-ci.d/code-quality.yaml b/.gitlab-ci.d/code-quality.yaml index 73a000acd..0a89fbe0d 100644 --- a/.gitlab-ci.d/code-quality.yaml +++ b/.gitlab-ci.d/code-quality.yaml @@ -1,63 +1,63 @@ .code_quality_common: stage: code_quality allow_failure: true rules: - if: "$CODE_QUALITY_DISABLED" when: never - if: "$CI_COMMIT_TAG || $CI_COMMIT_BRANCH" .code_quality_gitlab_template: extends: - .code_quality_common - image: docker:19.03.12 + image: docker:20.10.18 allow_failure: true services: - - docker:19.03.12-dind + - docker:20.10.18-dind variables: DOCKER_DRIVER: overlay2 DOCKER_TLS_CERTDIR: "" - CODE_QUALITY_IMAGE: "registry.gitlab.com/gitlab-org/ci-cd/codequality:0.85.24" + CODE_QUALITY_IMAGE: "registry.gitlab.com/gitlab-org/ci-cd/codequality:0.87.0" needs: [] script: - export SOURCE_CODE=$PWD - | if ! docker info &>/dev/null; then if [ -z "$DOCKER_HOST" -a "$KUBERNETES_PORT" ]; then export DOCKER_HOST='tcp://localhost:2375' fi fi - | # this is required to avoid undesirable reset of Docker image ENV variables being set on build stage function propagate_env_vars() { CURRENT_ENV=$(printenv) for VAR_NAME; do echo $CURRENT_ENV | grep "${VAR_NAME}=" > /dev/null && echo "--env $VAR_NAME " done } - docker pull --quiet "$CODE_QUALITY_IMAGE" - | docker run --rm \ $(propagate_env_vars \ SOURCE_CODE \ TIMEOUT_SECONDS \ CODECLIMATE_DEBUG \ CODECLIMATE_DEV \ REPORT_STDOUT \ REPORT_FORMAT \ ENGINE_MEMORY_LIMIT_BYTES \ CODECLIMATE_PREFIX \ ) \ --volume "$PWD":/code \ --volume /var/run/docker.sock:/var/run/docker.sock \ "$CODE_QUALITY_IMAGE" /code .clang_tools: extends: - .code_quality_common - .debian_bullseye_clang before_script: - if [ 'x${CI_MERGE_REQUEST_ID}' != 'x' ]; then - git fetch origin $CI_MERGE_REQUEST_TARGET_BRANCH_NAME - git diff --name-only $CI_COMMIT_SHA $CI_MERGE_REQUEST_TARGET_BRANCH_NAME > file_list - FILE_LIST_ARG='-f file_list' - fi 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_mesh.cc b/python/py_mesh.cc index 4f2972024..3c4de188a 100644 --- a/python/py_mesh.cc +++ b/python/py_mesh.cc @@ -1,240 +1,254 @@ /** * @file py_mesh.cc * * @author Guillaume Anciaux * @author Philip Mueller * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Sun Jun 16 2019 * @date last modification: Mon Mar 15 2021 * * @brief pybind11 interface to Mesh * * * @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 "aka_config.hh" /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { namespace { /* ------------------------------------------------------------------------ */ template decltype(auto) register_element_type_map_array(py::module & mod, const std::string & name) { return py::class_, std::shared_ptr>>( mod, ("ElementTypeMapArray" + name).c_str()) .def(py::init(), py::arg("id") = "by_element_type_array", py::arg("parent_id") = "no_parent") .def( "__call__", [](ElementTypeMapArray & self, ElementType type, GhostType ghost_type) -> decltype(auto) { return self(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference, py::keep_alive<0, 1>()) .def( "elementTypes", [](ElementTypeMapArray & self, UInt _dim, GhostType _ghost_type, ElementKind _kind) -> std::vector { auto types = self.elementTypes(_dim, _ghost_type, _kind); std::vector _types; for (auto && t : types) { _types.push_back(t); } return _types; }, py::arg("dim") = _all_dimensions, py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_regular) .def( "initialize", [](ElementTypeMapArray & self, const Mesh & mesh, GhostType ghost_type = _casper, UInt nb_component = 1, UInt spatial_dimension = UInt(-2), ElementKind element_kind = _ek_not_defined, bool with_nb_element = false, bool with_nb_nodes_per_element = false, T default_value = T(), bool do_not_default = false) { self.initialize( mesh, _ghost_type = ghost_type, _nb_component = nb_component, _spatial_dimension = (spatial_dimension == UInt(-2) ? mesh.getSpatialDimension() : spatial_dimension), _element_kind = element_kind, _with_nb_element = with_nb_element, _with_nb_nodes_per_element = with_nb_nodes_per_element, _default_value = default_value, _do_not_default = do_not_default); }, py::arg("mesh"), py::arg("ghost_type") = _casper, py::arg("nb_component") = 1, py::arg("spatial_dimension") = UInt(-2), py::arg("element_kind") = _ek_not_defined, py::arg("with_nb_element") = false, py::arg("with_nb_nodes_per_element") = false, py::arg("default_value") = T(), py::arg("do_not_default") = false); } } // namespace /* -------------------------------------------------------------------------- */ void register_mesh(py::module & mod) { py::class_(mod, "PeriodicSlaves") .def( "__iter__", [](Mesh::PeriodicSlaves & _this) { py::make_iterator(_this.begin(), _this.end()); }, py::keep_alive<0, 1>()); py::class_(mod, "MeshData") .def( "getElementalDataUInt", [](MeshData & _this, const ID & name) -> decltype(auto) { return _this.getElementalData(name); }, py::return_value_policy::reference) .def( "getElementalDataReal", [](MeshData & _this, const ID & name) -> decltype(auto) { return _this.getElementalData(name); }, py::return_value_policy::reference); py::class_(mod, "Mesh", py::multiple_inheritance()) .def(py::init(), py::arg("spatial_dimension"), py::arg("id") = "mesh") .def("read", &Mesh::read, py::arg("filename"), py::arg("mesh_io_type") = _miot_auto, "read the mesh from a file") .def( "getNodes", [](Mesh & self) -> decltype(auto) { return self.getNodes(); }, py::return_value_policy::reference) .def("getNbNodes", &Mesh::getNbNodes) .def( "getConnectivity", [](Mesh & self, ElementType type) -> decltype(auto) { return self.getConnectivity(type); }, py::return_value_policy::reference) + .def( + "getConnectivities", + [](Mesh & self) -> decltype(auto) { + return self.getConnectivities(); + }, + py::return_value_policy::reference) .def( "addConnectivityType", [](Mesh & self, ElementType type, GhostType ghost_type) -> void { self.addConnectivityType(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost) .def("distribute", [](Mesh & self) { self.distribute(); }) - .def("isDistributed", [](const Mesh& self) { return self.isDistributed(); }) + .def("isDistributed", + [](const Mesh & self) { return self.isDistributed(); }) .def("fillNodesToElements", &Mesh::fillNodesToElements, py::arg("dimension") = _all_dimensions) .def("getAssociatedElements", [](Mesh & self, const UInt & node, py::list list) { Array elements; self.getAssociatedElements(node, elements); for (auto && element : elements) { list.append(element); } }) .def("makePeriodic", [](Mesh & self, const SpatialDirection & direction) { self.makePeriodic(direction); }) .def( "getNbElement", [](Mesh & self, const UInt spatial_dimension, GhostType ghost_type, ElementKind kind) { return self.getNbElement(spatial_dimension, ghost_type, kind); }, py::arg("spatial_dimension") = _all_dimensions, py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_not_defined) .def( "getNbElement", [](Mesh & self, ElementType type, GhostType ghost_type) { return self.getNbElement(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost) .def_static( "getSpatialDimension", [](ElementType & type) { return Mesh::getSpatialDimension(type); }) .def( "getDataReal", [](Mesh & _this, const ID & name, ElementType type, GhostType ghost_type) -> decltype(auto) { return _this.getData(name, type, ghost_type); }, py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "hasDataReal", [](Mesh & _this, const ID & name, ElementType type, GhostType ghost_type) -> bool { return _this.hasData(name, type, ghost_type); }, py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def("isPeriodic", [](const Mesh & _this) { return _this.isPeriodic(); }) .def("getPeriodicMaster", &Mesh::getPeriodicMaster) .def("getPeriodicSlaves", &Mesh::getPeriodicSlaves) .def("isPeriodicSlave", &Mesh::isPeriodicSlave) .def("isPeriodicMaster", &Mesh::isPeriodicMaster) + .def( + "getMeshFacets", + [](const Mesh & self) -> const Mesh & { + return self.getMeshFacets(); + }, + py::return_value_policy::reference) .def("initMeshFacets", &Mesh::initMeshFacets, py::arg("id") = "mesh_facets", py::return_value_policy::reference); /* ------------------------------------------------------------------------ */ py::class_(mod, "MeshUtils") .def_static("buildFacets", &MeshUtils::buildFacets); py::class_(mod, "MeshAccessor") .def(py::init(), py::arg("mesh")) .def( "resizeConnectivity", [](MeshAccessor & self, UInt new_size, ElementType type, GhostType gt) -> void { self.resizeConnectivity(new_size, type, gt); }, py::arg("new_size"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "resizeNodes", [](MeshAccessor & self, UInt new_size) -> void { self.resizeNodes(new_size); }, py::arg("new_size")) .def("makeReady", &MeshAccessor::makeReady); register_element_type_map_array(mod, "Real"); register_element_type_map_array(mod, "UInt"); + register_element_type_map_array(mod, "bool"); // register_element_type_map_array(mod, "String"); } } // namespace akantu diff --git a/python/py_solid_mechanics_model_cohesive.cc b/python/py_solid_mechanics_model_cohesive.cc index ba64b18c5..81dce61dc 100644 --- a/python/py_solid_mechanics_model_cohesive.cc +++ b/python/py_solid_mechanics_model_cohesive.cc @@ -1,96 +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("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/common/aka_common.cc b/src/common/aka_common.cc index 4dca7b03b..f15de9893 100644 --- a/src/common/aka_common.cc +++ b/src/common/aka_common.cc @@ -1,156 +1,160 @@ /** * @file aka_common.cc * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Wed Dec 09 2020 * * @brief Initialization of global variables * * * @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 "aka_common.hh" #include "aka_random_generator.hh" #include "communicator.hh" #include "cppargparse.hh" #include "parser.hh" #include "communication_tag.hh" /* -------------------------------------------------------------------------- */ #include +#include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void initialize(int & argc, char **& argv) { AKANTU_DEBUG_IN(); initialize("", argc, argv); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void initialize(const std::string & input_file, int & argc, char **& argv) { AKANTU_DEBUG_IN(); Communicator & comm = Communicator::getStaticCommunicator(argc, argv); Tag::setMaxTag(comm.getMaxTag()); debug::debugger.setParallelContext(comm.whoAmI(), comm.getNbProc()); debug::setDebugLevel(dblError); static_argparser.setParallelContext(comm.whoAmI(), comm.getNbProc()); static_argparser.setExternalExitFunction(debug::exit); static_argparser.addArgument("--aka_input_file", "Akantu's input file", 1, cppargparse::_string, std::string()); static_argparser.addArgument( "--aka_debug_level", std::string("Akantu's overall debug level") + std::string(" (0: error, 1: exceptions, 4: warnings, 5: info, ..., " "100: dump") + std::string(" more info on levels can be foind in aka_error.hh)"), 1, cppargparse::_integer, (long int)(dblWarning)); static_argparser.addArgument( "--aka_print_backtrace", "Should Akantu print a backtrace in case of error", 0, cppargparse::_boolean, false, true); static_argparser.addArgument("--aka_seed", "The seed to use on prank 0", 1, cppargparse::_integer); static_argparser.parse(argc, argv, cppargparse::_remove_parsed); std::string infile = static_argparser["aka_input_file"]; if (infile.empty()) { infile = input_file; } debug::debugger.printBacktrace(static_argparser["aka_print_backtrace"]); if (not infile.empty()) { readInputFile(infile); } long int seed; - if (static_argparser.has("aka_seed")) { + char * env_seed = std::getenv("AKA_SEED"); + if (env_seed != nullptr) { + seed = std::atol(env_seed); + } else if (static_argparser.has("aka_seed")) { seed = static_argparser["aka_seed"]; } else { seed = static_parser.getParameter("seed", time(nullptr), _ppsc_current_scope); } seed *= (comm.whoAmI() + 1); RandomGenerator::seed(seed); long int dbl_level = static_argparser["aka_debug_level"]; debug::setDebugLevel(DebugLevel(dbl_level)); AKANTU_DEBUG_INFO("Random seed set to " << seed); std::atexit(finalize); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void finalize() {} /* -------------------------------------------------------------------------- */ void readInputFile(const std::string & input_file) { static_parser.parse(input_file); } /* -------------------------------------------------------------------------- */ cppargparse::ArgumentParser & getStaticArgumentParser() { return static_argparser; } /* -------------------------------------------------------------------------- */ Parser & getStaticParser() { return static_parser; } /* -------------------------------------------------------------------------- */ const ParserSection & getUserParser() { return *(static_parser.getSubSections(ParserType::_user).first); } std::unique_ptr Communicator::static_communicator; std::ostream & operator<<(std::ostream & stream, NodeFlag flag) { using under = std::underlying_type_t; auto digits = static_cast( std::log(std::numeric_limits::max() + 1) / std::log(16)); std::ios_base::fmtflags ff; ff = stream.flags(); auto value = static_cast>(flag); stream << "0x" << std::hex << std::setw(digits) << std::setfill('0') << value; stream.flags(ff); return stream; } } // namespace akantu diff --git a/src/common/aka_grid_dynamic.hh b/src/common/aka_grid_dynamic.hh index c0386aa45..668da3734 100644 --- a/src/common/aka_grid_dynamic.hh +++ b/src/common/aka_grid_dynamic.hh @@ -1,532 +1,530 @@ /** * @file aka_grid_dynamic.hh * * @author Aurelia Isabel Cuba Ramos * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Thu Feb 21 2013 * @date last modification: Tue Feb 09 2021 * * @brief Grid that is auto balanced * * * @section LICENSE * * Copyright (©) 2014-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 "aka_array.hh" #include "aka_common.hh" #include "aka_types.hh" #include "mesh_accessor.hh" #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_AKA_GRID_DYNAMIC_HH_ #define AKANTU_AKA_GRID_DYNAMIC_HH_ namespace akantu { class Mesh; template class SpatialGrid { public: explicit SpatialGrid(UInt dimension) : dimension(dimension), spacing(dimension), center(dimension), lower(dimension), upper(dimension), empty_cell() {} SpatialGrid(UInt dimension, const Vector & spacing, const Vector & center) : dimension(dimension), spacing(spacing), center(center), lower(dimension), upper(dimension), empty_cell() { for (UInt i = 0; i < dimension; ++i) { lower(i) = std::numeric_limits::max(); upper(i) = -std::numeric_limits::max(); } } virtual ~SpatialGrid() = default; class neighbor_cells_iterator; class cells_iterator; class CellID { public: CellID() = default; explicit CellID(UInt dimention) : ids(dimention) {} void setID(UInt dir, Int id) { ids(dir) = id; } Int getID(UInt dir) const { return ids(dir); } bool operator<(const CellID & id) const { return std::lexicographical_compare( ids.storage(), ids.storage() + ids.size(), id.ids.storage(), id.ids.storage() + id.ids.size()); } bool operator==(const CellID & id) const { return std::equal(ids.storage(), ids.storage() + ids.size(), id.ids.storage()); } bool operator!=(const CellID & id) const { return !(operator==(id)); } - class neighbor_cells_iterator - : private std::iterator { + class neighbor_cells_iterator { public: neighbor_cells_iterator(const CellID & cell_id, bool end) : cell_id(cell_id), position(cell_id.ids.size(), end ? 1 : -1) { this->updateIt(); if (end) { this->it++; } } neighbor_cells_iterator & operator++() { UInt i = 0; for (; i < position.size() && position(i) == 1; ++i) { ; } if (i == position.size()) { ++it; return *this; } for (UInt j = 0; j < i; ++j) { position(j) = -1; } position(i)++; updateIt(); return *this; } neighbor_cells_iterator operator++(int) { neighbor_cells_iterator tmp(*this); operator++(); return tmp; }; bool operator==(const neighbor_cells_iterator & rhs) const { return cell_id == rhs.cell_id && it == rhs.it; }; bool operator!=(const neighbor_cells_iterator & rhs) const { return !operator==(rhs); }; CellID operator*() const { CellID cur_cell_id(cell_id); cur_cell_id.ids += position; return cur_cell_id; }; private: void updateIt() { it = 0; for (UInt i = 0; i < position.size(); ++i) { it = it * 3 + (position(i) + 1); } } private: /// central cell id const CellID & cell_id; // number representing the current neighbor in base 3; UInt it; // current cell shift Vector position; }; class Neighbors { public: explicit Neighbors(const CellID & cell_id) : cell_id(cell_id) {} decltype(auto) begin() { return neighbor_cells_iterator(cell_id, false); } decltype(auto) end() { return neighbor_cells_iterator(cell_id, true); } private: const CellID & cell_id; }; decltype(auto) neighbors() { return Neighbors(*this); } private: friend class cells_iterator; Vector ids; }; /* ------------------------------------------------------------------------ */ class Cell { public: using iterator = typename std::vector::iterator; using const_iterator = typename std::vector::const_iterator; Cell() : id(), data() {} explicit Cell(const CellID & cell_id) : id(cell_id), data() {} bool operator==(const Cell & cell) const { return id == cell.id; } bool operator!=(const Cell & cell) const { return id != cell.id; } Cell & add(const T & d) { data.push_back(d); return *this; } iterator begin() { return data.begin(); } const_iterator begin() const { return data.begin(); } iterator end() { return data.end(); } const_iterator end() const { return data.end(); } private: CellID id; std::vector data; }; private: using cells_container = std::map; public: const Cell & getCell(const CellID & cell_id) const { auto it = cells.find(cell_id); if (it != cells.end()) { return it->second; } return empty_cell; } decltype(auto) beginCell(const CellID & cell_id) { auto it = cells.find(cell_id); if (it != cells.end()) { return it->second.begin(); } return empty_cell.begin(); } decltype(auto) endCell(const CellID & cell_id) { auto it = cells.find(cell_id); if (it != cells.end()) { return it->second.end(); } return empty_cell.end(); } decltype(auto) beginCell(const CellID & cell_id) const { auto it = cells.find(cell_id); if (it != cells.end()) { return it->second.begin(); } return empty_cell.begin(); } decltype(auto) endCell(const CellID & cell_id) const { auto it = cells.find(cell_id); if (it != cells.end()) { return it->second.end(); } return empty_cell.end(); } /* ------------------------------------------------------------------------ */ - class cells_iterator - : private std::iterator { + class cells_iterator { public: explicit cells_iterator(typename std::map::const_iterator it) : it(it) {} cells_iterator & operator++() { this->it++; return *this; } cells_iterator operator++(int /*unused*/) { cells_iterator tmp(*this); operator++(); return tmp; }; bool operator==(const cells_iterator & rhs) const { return it == rhs.it; }; bool operator!=(const cells_iterator & rhs) const { return !operator==(rhs); }; CellID operator*() const { CellID cur_cell_id(this->it->first); return cur_cell_id; }; private: /// map iterator typename std::map::const_iterator it; }; public: template Cell & insert(const T & d, const vector_type & position) { auto && cell_id = getCellID(position); auto && it = cells.find(cell_id); if (it == cells.end()) { Cell cell(cell_id); auto & tmp = (cells[cell_id] = cell).add(d); for (UInt i = 0; i < dimension; ++i) { Real posl = center(i) + cell_id.getID(i) * spacing(i); Real posu = posl + spacing(i); if (posl <= lower(i)) { lower(i) = posl; } if (posu > upper(i)) { upper(i) = posu; } } return tmp; } return it->second.add(d); } /* ------------------------------------------------------------------------ */ inline decltype(auto) begin() const { auto begin = this->cells.begin(); return cells_iterator(begin); } inline decltype(auto) end() const { auto end = this->cells.end(); return cells_iterator(end); } template CellID getCellID(const vector_type & position) const { CellID cell_id(dimension); for (UInt i = 0; i < dimension; ++i) { cell_id.setID(i, getCellID(position(i), i)); } return cell_id; } void printself(std::ostream & stream, int indent = 0) const { std::string space(indent, AKANTU_INDENT); std::streamsize prec = stream.precision(); std::ios_base::fmtflags ff = stream.flags(); stream.setf(std::ios_base::showbase); stream.precision(5); stream << space << "SpatialGrid<" << debug::demangle(typeid(T).name()) << "> [" << std::endl; stream << space << " + dimension : " << this->dimension << std::endl; stream << space << " + lower bounds : {"; for (UInt i = 0; i < lower.size(); ++i) { if (i != 0) { stream << ", "; } stream << lower(i); }; stream << "}" << std::endl; stream << space << " + upper bounds : {"; for (UInt i = 0; i < upper.size(); ++i) { if (i != 0) { stream << ", "; } stream << upper(i); }; stream << "}" << std::endl; stream << space << " + spacing : {"; for (UInt i = 0; i < spacing.size(); ++i) { if (i != 0) { stream << ", "; } stream << spacing(i); }; stream << "}" << std::endl; stream << space << " + center : {"; for (UInt i = 0; i < center.size(); ++i) { if (i != 0) { stream << ", "; } stream << center(i); }; stream << "}" << std::endl; stream << space << " + nb_cells : " << this->cells.size() << "/"; Vector dist(this->dimension); dist = upper; dist -= lower; for (UInt i = 0; i < this->dimension; ++i) { dist(i) /= spacing(i); } UInt nb_cells = std::ceil(dist(0)); for (UInt i = 1; i < this->dimension; ++i) { nb_cells *= std::ceil(dist(i)); } stream << nb_cells << std::endl; stream << space << "]" << std::endl; stream.precision(prec); stream.flags(ff); } void saveAsMesh(Mesh & mesh) const; private: /* -------------------------------------------------------------------------- */ inline UInt getCellID(Real position, UInt direction) const { AKANTU_DEBUG_ASSERT(direction < center.size(), "The direction asked (" << direction << ") is out of range " << center.size()); Real dist_center = position - center(direction); Int id = std::floor(dist_center / spacing(direction)); // if(dist_center < 0) id--; return id; } friend class GridSynchronizer; public: AKANTU_GET_MACRO(LowerBounds, lower, const Vector &); AKANTU_GET_MACRO(UpperBounds, upper, const Vector &); AKANTU_GET_MACRO(Spacing, spacing, const Vector &); AKANTU_SET_MACRO(Spacing, spacing, Vector &); AKANTU_GET_MACRO(Center, center, const Vector &); AKANTU_SET_MACRO(Center, center, Vector &); protected: UInt dimension; cells_container cells; Vector spacing; Vector center; Vector lower; Vector upper; Cell empty_cell; }; /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const SpatialGrid & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "mesh.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template void SpatialGrid::saveAsMesh(Mesh & mesh) const { ElementType type = _not_defined; switch (dimension) { case 1: type = _segment_2; break; case 2: type = _quadrangle_4; break; case 3: type = _hexahedron_8; break; } MeshAccessor mesh_accessor(mesh); auto & connectivity = mesh_accessor.getConnectivity(type); auto & nodes = mesh_accessor.getNodes(); auto & uint_data = mesh.getDataPointer("tag_1", type); Vector pos(dimension); UInt global_id = 0; for (auto & cell_pair : cells) { UInt cur_node = nodes.size(); UInt cur_elem = connectivity.size(); const CellID & cell_id = cell_pair.first; for (UInt i = 0; i < dimension; ++i) { pos(i) = center(i) + cell_id.getID(i) * spacing(i); } nodes.push_back(pos); for (UInt i = 0; i < dimension; ++i) { pos(i) += spacing(i); } nodes.push_back(pos); connectivity.push_back(cur_node); switch (dimension) { case 1: connectivity(cur_elem, 1) = cur_node + 1; break; case 2: pos(0) -= spacing(0); nodes.push_back(pos); pos(0) += spacing(0); pos(1) -= spacing(1); nodes.push_back(pos); connectivity(cur_elem, 1) = cur_node + 3; connectivity(cur_elem, 2) = cur_node + 1; connectivity(cur_elem, 3) = cur_node + 2; break; case 3: pos(1) -= spacing(1); pos(2) -= spacing(2); nodes.push_back(pos); pos(1) += spacing(1); nodes.push_back(pos); pos(0) -= spacing(0); nodes.push_back(pos); pos(1) -= spacing(1); pos(2) += spacing(2); nodes.push_back(pos); pos(0) += spacing(0); nodes.push_back(pos); pos(0) -= spacing(0); pos(1) += spacing(1); nodes.push_back(pos); connectivity(cur_elem, 1) = cur_node + 2; connectivity(cur_elem, 2) = cur_node + 3; connectivity(cur_elem, 3) = cur_node + 4; connectivity(cur_elem, 4) = cur_node + 5; connectivity(cur_elem, 5) = cur_node + 6; connectivity(cur_elem, 6) = cur_node + 1; connectivity(cur_elem, 7) = cur_node + 7; break; } uint_data.push_back(global_id); ++global_id; } } } // namespace akantu #endif /* AKANTU_AKA_GRID_DYNAMIC_HH_ */ diff --git a/src/common/aka_math_tmpl.hh b/src/common/aka_math_tmpl.hh index a27337a82..d5b37876f 100644 --- a/src/common/aka_math_tmpl.hh +++ b/src/common/aka_math_tmpl.hh @@ -1,844 +1,843 @@ /** * @file aka_math_tmpl.hh * * @author Ramin Aghababaei * @author Guillaume Anciaux * @author Alejandro M. Aragón * @author Emil Gallyamov * @author David Simon Kammer * @author Daniel Pino Muñoz * @author Mohit Pundir * @author Mathilde Radiguet * @author Nicolas Richart * @author Leonardo Snozzi * @author Peter Spijker * @author Marco Vocialta * * @date creation: Wed Aug 04 2010 * @date last modification: Fri Dec 11 2020 * * @brief Implementation of the inline functions of the math toolkit * * * @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 "aka_blas_lapack.hh" #include "aka_math.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { namespace Math { /* ------------------------------------------------------------------------ */ inline void matrix_vector(UInt im, UInt in, Real * A, // NOLINT(readability-non-const-parameter) Real * x, // NOLINT(readability-non-const-parameter) Real * y, Real alpha) { #ifdef AKANTU_USE_BLAS /// y = alpha*op(A)*x + beta*y char tran_A = 'N'; int incx = 1; int incy = 1; double beta = 0.; int m = im; int n = in; aka_gemv(&tran_A, &m, &n, &alpha, A, &m, x, &incx, &beta, y, &incy); #else std::fill_n(y, im, 0.); for (UInt i = 0; i < im; ++i) { for (UInt j = 0; j < in; ++j) { y[i] += A[i + j * im] * x[j]; } y[i] *= alpha; } #endif } /* ------------------------------------------------------------------------ */ inline void matrixt_vector(UInt im, UInt in, Real * A, // NOLINT(readability-non-const-parameter) Real * x, // NOLINT(readability-non-const-parameter) Real * y, Real alpha) { #ifdef AKANTU_USE_BLAS /// y = alpha*op(A)*x + beta*y char tran_A = 'T'; int incx = 1; int incy = 1; double beta = 0.; int m = im; int n = in; aka_gemv(&tran_A, &m, &n, &alpha, A, &m, x, &incx, &beta, y, &incy); #else std::fill_n(y, in, 0.); for (UInt i = 0; i < im; ++i) { for (UInt j = 0; j < in; ++j) { - y[j] += A[j * im + i] * x[i]; + y[j] += alpha * A[j * im + i] * x[i]; } - y[i] *= alpha; } #endif } /* ------------------------------------------------------------------------ */ inline void matrix_matrix(UInt im, UInt in, UInt ik, Real * A, // NOLINT(readability-non-const-parameter) Real * B, // NOLINT(readability-non-const-parameter) Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C char trans_a = 'N'; char trans_b = 'N'; double beta = 0.; int m = im, n = in, k = ik; aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &m, B, &k, &beta, C, &m); #else std::fill_n(C, im * in, 0.); for (UInt j = 0; j < in; ++j) { UInt _jb = j * ik; UInt _jc = j * im; for (UInt i = 0; i < im; ++i) { for (UInt l = 0; l < ik; ++l) { UInt _la = l * im; C[i + _jc] += A[i + _la] * B[l + _jb]; } C[i + _jc] *= alpha; } } #endif } /* ------------------------------------------------------------------------ */ inline void matrixt_matrix(UInt im, UInt in, UInt ik, Real * A, // NOLINT(readability-non-const-parameter) Real * B, // NOLINT(readability-non-const-parameter) Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C char trans_a = 'T'; char trans_b = 'N'; double beta = 0.; int m = im, n = in, k = ik; aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &k, B, &k, &beta, C, &m); #else std::fill_n(C, im * in, 0.); for (UInt j = 0; j < in; ++j) { UInt _jc = j * im; UInt _jb = j * ik; for (UInt i = 0; i < im; ++i) { UInt _ia = i * ik; for (UInt l = 0; l < ik; ++l) { C[i + _jc] += A[l + _ia] * B[l + _jb]; } C[i + _jc] *= alpha; } } #endif } /* ------------------------------------------------------------------------ */ inline void matrix_matrixt(UInt im, UInt in, UInt ik, Real * A, // NOLINT(readability-non-const-parameter) Real * B, // NOLINT(readability-non-const-parameter) Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C char trans_a = 'N'; char trans_b = 'T'; double beta = 0.; int m = im, n = in, k = ik; aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &m, B, &n, &beta, C, &m); #else std::fill_n(C, im * in, 0.); for (UInt j = 0; j < in; ++j) { UInt _jc = j * im; for (UInt i = 0; i < im; ++i) { for (UInt l = 0; l < ik; ++l) { UInt _la = l * im; UInt _lb = l * in; C[i + _jc] += A[i + _la] * B[j + _lb]; } C[i + _jc] *= alpha; } } #endif } /* ------------------------------------------------------------------------ */ inline void matrixt_matrixt(UInt im, UInt in, UInt ik, Real * A, // NOLINT(readability-non-const-parameter) Real * B, // NOLINT(readability-non-const-parameter) Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C char trans_a = 'T'; char trans_b = 'T'; double beta = 0.; int m = im, n = in, k = ik; aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &k, B, &n, &beta, C, &m); #else std::fill_n(C, im * in, 0.); for (UInt j = 0; j < in; ++j) { UInt _jc = j * im; for (UInt i = 0; i < im; ++i) { UInt _ia = i * ik; for (UInt l = 0; l < ik; ++l) { UInt _lb = l * in; C[i + _jc] += A[l + _ia] * B[j + _lb]; } C[i + _jc] *= alpha; } } #endif } /* ------------------------------------------------------------------------ */ inline void aXplusY(UInt n, Real alpha, Real * x, // NOLINT(readability-non-const-parameter) Real * y) { #ifdef AKANTU_USE_BLAS /// y := alpha x + y int incx = 1, incy = 1; aka_axpy(&n, &alpha, x, &incx, y, &incy); #else for (UInt i = 0; i < n; ++i) { *(y++) += alpha * *(x++); } #endif } /* ------------------------------------------------------------------------ */ inline Real vectorDot(Real * v1, // NOLINT(readability-non-const-parameter) Real * v2, // NOLINT(readability-non-const-parameter) UInt in) { #ifdef AKANTU_USE_BLAS /// d := v1 . v2 int incx = 1, incy = 1, n = in; Real d = aka_dot(&n, v1, &incx, v2, &incy); #else Real d = 0; for (UInt i = 0; i < in; ++i) { d += v1[i] * v2[i]; } #endif return d; } /* ------------------------------------------------------------------------ */ template inline void matMul(UInt m, UInt n, UInt k, Real alpha, Real * A, // NOLINT(readability-non-const-parameter) Real * B, // NOLINT(readability-non-const-parameter) Real /*beta*/, Real * C) { if (tr_A) { if (tr_B) { matrixt_matrixt(m, n, k, A, B, C, alpha); } else { matrixt_matrix(m, n, k, A, B, C, alpha); } } else { if (tr_B) { matrix_matrixt(m, n, k, A, B, C, alpha); } else { matrix_matrix(m, n, k, A, B, C, alpha); } } } /* ------------------------------------------------------------------------ */ template inline void matVectMul(UInt m, UInt n, Real alpha, Real * A, // NOLINT(readability-non-const-parameter) Real * x, // NOLINT(readability-non-const-parameter) Real /*beta*/, Real * y) { if (tr_A) { matrixt_vector(m, n, A, x, y, alpha); } else { matrix_vector(m, n, A, x, y, alpha); } } /* ------------------------------------------------------------------------ */ template inline void matrixEig(UInt n, T * A, // NOLINT(readability-non-const-parameter) T * d, T * V) { // Matrix A is row major, so the lapack function in fortran will // process A^t. Asking for the left eigenvectors of A^t will give the // transposed right eigenvectors of A so in the C++ code the right // eigenvectors. char jobvr{'N'}; if (V != nullptr) { jobvr = 'V'; // compute left eigenvectors } char jobvl{'N'}; // compute right eigenvectors auto * di = new T[n]; // imaginary part of the eigenvalues int info; int N = n; T wkopt; int lwork = -1; // query and allocate the optimal workspace aka_geev(&jobvl, &jobvr, &N, A, &N, d, di, nullptr, &N, V, &N, &wkopt, &lwork, &info); lwork = int(wkopt); auto * work = new T[lwork]; // solve the eigenproblem aka_geev(&jobvl, &jobvr, &N, A, &N, d, di, nullptr, &N, V, &N, work, &lwork, &info); AKANTU_DEBUG_ASSERT( info == 0, "Problem computing eigenvalues/vectors. DGEEV exited with the value " << info); delete[] work; delete[] di; // I hope for you that there was no complex eigenvalues !!! } /* ------------------------------------------------------------------------ */ inline void matrix22_eigenvalues(Real * A, // NOLINT(readability-non-const-parameter) Real * Adiag) { /// d = determinant of Matrix A Real d = det2(A); /// b = trace of Matrix A Real b = A[0] + A[3]; Real c = std::sqrt(b * b - 4 * d); Adiag[0] = .5 * (b + c); Adiag[1] = .5 * (b - c); } /* ------------------------------------------------------------------------ */ inline void matrix33_eigenvalues(Real * A, // NOLINT(readability-non-const-parameter) Real * Adiag) { matrixEig(3, A, Adiag); } /* ------------------------------------------------------------------------ */ template inline void eigenvalues(Real * A, // NOLINT(readability-non-const-parameter) Real * d) { if (dim == 1) { d[0] = A[0]; } else if (dim == 2) { matrix22_eigenvalues(A, d); } // else if(dim == 3) { matrix33_eigenvalues(A, d); } else { matrixEig(dim, A, d); } } /* ------------------------------------------------------------------------ */ inline Real det2(const Real * mat) { return mat[0] * mat[3] - mat[1] * mat[2]; } /* ------------------------------------------------------------------------ */ inline Real det3(const Real * mat) { return mat[0] * (mat[4] * mat[8] - mat[7] * mat[5]) - mat[3] * (mat[1] * mat[8] - mat[7] * mat[2]) + mat[6] * (mat[1] * mat[5] - mat[4] * mat[2]); } /* ------------------------------------------------------------------------ */ template inline Real det(const Real * mat) { if (n == 1) { return *mat; } if (n == 2) { return det2(mat); } if (n == 3) { return det3(mat); } return det(n, mat); } /* ------------------------------------------------------------------------ */ template inline T det(UInt n, const T * A) { int N = n; int info; auto * ipiv = new int[N + 1]; auto * LU = new T[N * N]; std::copy(A, A + N * N, LU); // LU factorization of A aka_getrf(&N, &N, LU, &N, ipiv, &info); if (info > 0) { AKANTU_ERROR("Singular matrix - cannot factorize it (info: " << info << " )"); } // det(A) = det(L) * det(U) = 1 * det(U) = product_i U_{ii} T det = 1.; for (int i = 0; i < N; ++i) { det *= (2 * (ipiv[i] == i) - 1) * LU[i * n + i]; } delete[] ipiv; delete[] LU; return det; } /* ------------------------------------------------------------------------ */ inline void normal2(const Real * vec, Real * normal) { normal[0] = vec[1]; normal[1] = -vec[0]; normalize2(normal); } /* ------------------------------------------------------------------------ */ inline void normal3(const Real * vec1, const Real * vec2, Real * normal) { vectorProduct3(vec1, vec2, normal); normalize3(normal); } /* ------------------------------------------------------------------------ */ inline void normalize2(Real * vec) { Real norm = norm2(vec); vec[0] /= norm; vec[1] /= norm; } /* ------------------------------------------------------------------------ */ inline void normalize3(Real * vec) { Real norm = norm3(vec); vec[0] /= norm; vec[1] /= norm; vec[2] /= norm; } /* ------------------------------------------------------------------------ */ inline Real norm2(const Real * vec) { return sqrt(vec[0] * vec[0] + vec[1] * vec[1]); } /* ------------------------------------------------------------------------ */ inline Real norm3(const Real * vec) { return sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]); } /* ------------------------------------------------------------------------ */ inline Real norm(UInt n, const Real * vec) { Real norm = 0.; for (UInt i = 0; i < n; ++i) { norm += vec[i] * vec[i]; } return sqrt(norm); } /* ------------------------------------------------------------------------ */ inline void inv2(const Real * mat, Real * inv) { Real det_mat = det2(mat); inv[0] = mat[3] / det_mat; inv[1] = -mat[1] / det_mat; inv[2] = -mat[2] / det_mat; inv[3] = mat[0] / det_mat; } /* ------------------------------------------------------------------------ */ inline void inv3(const Real * mat, Real * inv) { Real det_mat = det3(mat); inv[0] = (mat[4] * mat[8] - mat[7] * mat[5]) / det_mat; inv[1] = (mat[2] * mat[7] - mat[8] * mat[1]) / det_mat; inv[2] = (mat[1] * mat[5] - mat[4] * mat[2]) / det_mat; inv[3] = (mat[5] * mat[6] - mat[8] * mat[3]) / det_mat; inv[4] = (mat[0] * mat[8] - mat[6] * mat[2]) / det_mat; inv[5] = (mat[2] * mat[3] - mat[5] * mat[0]) / det_mat; inv[6] = (mat[3] * mat[7] - mat[6] * mat[4]) / det_mat; inv[7] = (mat[1] * mat[6] - mat[7] * mat[0]) / det_mat; inv[8] = (mat[0] * mat[4] - mat[3] * mat[1]) / det_mat; } /* ------------------------------------------------------------------------ */ template inline void inv(const Real * A, Real * Ainv) { if (n == 1) { *Ainv = 1. / *A; } else if (n == 2) { inv2(A, Ainv); } else if (n == 3) { inv3(A, Ainv); } else { inv(n, A, Ainv); } } /* ------------------------------------------------------------------------ */ template inline void inv(UInt n, const T * A, T * invA) { int N = n; int info; auto * ipiv = new int[N + 1]; int lwork = N * N; auto * work = new T[lwork]; std::copy(A, A + n * n, invA); aka_getrf(&N, &N, invA, &N, ipiv, &info); if (info > 0) { AKANTU_ERROR("Singular matrix - cannot factorize it (info: " << info << " )"); } aka_getri(&N, invA, &N, ipiv, work, &lwork, &info); if (info != 0) { AKANTU_ERROR("Cannot invert the matrix (info: " << info << " )"); } delete[] ipiv; delete[] work; } /* ------------------------------------------------------------------------ */ template inline void solve(UInt n, const T * A, T * x, const T * b) { int N = n; int info; auto * ipiv = new int[N]; auto * lu_A = new T[N * N]; std::copy(A, A + N * N, lu_A); aka_getrf(&N, &N, lu_A, &N, ipiv, &info); if (info > 0) { AKANTU_ERROR("Singular matrix - cannot factorize it (info: " << info << " )"); } char trans = 'N'; int nrhs = 1; std::copy(b, b + N, x); aka_getrs(&trans, &N, &nrhs, lu_A, &N, ipiv, x, &N, &info); if (info != 0) { AKANTU_ERROR("Cannot solve the system (info: " << info << " )"); } delete[] ipiv; delete[] lu_A; } /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ inline Real matrixDoubleDot22(Real * A, // NOLINT(readability-non-const-parameter) Real * B // NOLINT(readability-non-const-parameter) ) { Real d; d = A[0] * B[0] + A[1] * B[1] + A[2] * B[2] + A[3] * B[3]; return d; } /* ------------------------------------------------------------------------ */ inline Real matrixDoubleDot33(Real * A, // NOLINT(readability-non-const-parameter) Real * B // NOLINT(readability-non-const-parameter) ) { Real d; d = A[0] * B[0] + A[1] * B[1] + A[2] * B[2] + A[3] * B[3] + A[4] * B[4] + A[5] * B[5] + A[6] * B[6] + A[7] * B[7] + A[8] * B[8]; return d; } /* ------------------------------------------------------------------------ */ inline Real matrixDoubleDot(UInt n, Real * A, // NOLINT(readability-non-const-parameter) Real * B // NOLINT(readability-non-const-parameter) ) { Real d = 0.; for (UInt i = 0; i < n; ++i) { for (UInt j = 0; j < n; ++j) { d += A[i * n + j] * B[i * n + j]; } } return d; } /* ------------------------------------------------------------------------ */ inline void vectorProduct3(const Real * v1, const Real * v2, Real * res) { res[0] = v1[1] * v2[2] - v1[2] * v2[1]; res[1] = v1[2] * v2[0] - v1[0] * v2[2]; res[2] = v1[0] * v2[1] - v1[1] * v2[0]; } /* ------------------------------------------------------------------------ */ inline Real vectorDot2(const Real * v1, const Real * v2) { return (v1[0] * v2[0] + v1[1] * v2[1]); } /* ------------------------------------------------------------------------ */ inline Real vectorDot3(const Real * v1, const Real * v2) { return (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]); } /* ------------------------------------------------------------------------ */ inline Real distance_2d(const Real * x, const Real * y) { return std::sqrt((y[0] - x[0]) * (y[0] - x[0]) + (y[1] - x[1]) * (y[1] - x[1])); } /* ------------------------------------------------------------------------ */ inline Real triangle_inradius(const Vector & coord1, const Vector & coord2, const Vector & coord3) { /** * @f{eqnarray*}{ * r &=& A / s \\ * A &=& 1/4 * \sqrt{(a + b + c) * (a - b + c) * (a + b - c) (-a + b + c)} * \\ s &=& \frac{a + b + c}{2} * @f} */ Real a, b, c; a = coord1.distance(coord2); b = coord2.distance(coord3); c = coord1.distance(coord3); Real s; s = (a + b + c) * 0.5; return std::sqrt((s - a) * (s - b) * (s - c) / s); } /* ------------------------------------------------------------------------ */ inline Real distance_3d(const Real * x, const Real * y) { return std::sqrt((y[0] - x[0]) * (y[0] - x[0]) + (y[1] - x[1]) * (y[1] - x[1]) + (y[2] - x[2]) * (y[2] - x[2])); } /* ------------------------------------------------------------------------ */ inline Real tetrahedron_volume(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4) { Real xx[9]; xx[0] = coord2[0]; xx[1] = coord2[1]; xx[2] = coord2[2]; xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; auto vol = det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol -= det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol += det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2]; xx[6] = coord3[0]; xx[7] = coord3[1]; xx[8] = coord3[2]; vol -= det3(xx); vol /= 6; return vol; } /* ------------------------------------------------------------------------ */ inline Real tetrahedron_inradius(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4) { auto l12 = distance_3d(coord1, coord2); auto l13 = distance_3d(coord1, coord3); auto l14 = distance_3d(coord1, coord4); auto l23 = distance_3d(coord2, coord3); auto l24 = distance_3d(coord2, coord4); auto l34 = distance_3d(coord3, coord4); auto s1 = (l12 + l23 + l13) * 0.5; s1 = std::sqrt(s1 * (s1 - l12) * (s1 - l23) * (s1 - l13)); auto s2 = (l12 + l24 + l14) * 0.5; s2 = std::sqrt(s2 * (s2 - l12) * (s2 - l24) * (s2 - l14)); auto s3 = (l23 + l34 + l24) * 0.5; s3 = std::sqrt(s3 * (s3 - l23) * (s3 - l34) * (s3 - l24)); auto s4 = (l13 + l34 + l14) * 0.5; s4 = std::sqrt(s4 * (s4 - l13) * (s4 - l34) * (s4 - l14)); auto volume = tetrahedron_volume(coord1, coord2, coord3, coord4); return 3 * volume / (s1 + s2 + s3 + s4); } /* ------------------------------------------------------------------------ */ inline void barycenter(const Real * coord, UInt nb_points, UInt spatial_dimension, Real * barycenter) { std::fill_n(barycenter, spatial_dimension, 0.); for (UInt n = 0; n < nb_points; ++n) { UInt offset = n * spatial_dimension; for (UInt i = 0; i < spatial_dimension; ++i) { barycenter[i] += coord[offset + i] / (Real)nb_points; } } } /* ------------------------------------------------------------------------ */ inline void vector_2d(const Real * x, const Real * y, Real * res) { res[0] = y[0] - x[0]; res[1] = y[1] - x[1]; } /* ------------------------------------------------------------------------ */ inline void vector_3d(const Real * x, const Real * y, Real * res) { res[0] = y[0] - x[0]; res[1] = y[1] - x[1]; res[2] = y[2] - x[2]; } /* ------------------------------------------------------------------------ */ /// Combined absolute and relative tolerance test proposed in /// Real-time collision detection by C. Ericson (2004) inline bool are_float_equal(const Real x, const Real y) { Real abs_max = std::max(std::abs(x), std::abs(y)); abs_max = std::max(abs_max, Real(1.)); return std::abs(x - y) <= (tolerance * abs_max); } /* ------------------------------------------------------------------------ */ inline bool isnan(Real x) { #if defined(__INTEL_COMPILER) #pragma warning(push) #pragma warning(disable : 1572) #endif // defined(__INTEL_COMPILER) // x = x return false means x = quiet_NaN return !(x == x); #if defined(__INTEL_COMPILER) #pragma warning(pop) #endif // defined(__INTEL_COMPILER) } /* ------------------------------------------------------------------------ */ inline bool are_vector_equal(UInt n, Real * x, Real * y) { bool test = true; for (UInt i = 0; i < n; ++i) { test &= are_float_equal(x[i], y[i]); } return test; } /* ------------------------------------------------------------------------ */ inline bool intersects(Real x_min, Real x_max, Real y_min, Real y_max) { return not((x_max < y_min) or (x_min > y_max)); } /* ------------------------------------------------------------------------ */ inline bool is_in_range(Real a, Real x_min, Real x_max) { return ((a >= x_min) and (a <= x_max)); } /* ------------------------------------------------------------------------ */ template inline T pow(T x) { return (pow

(x) * x); } template <> inline UInt pow<0, UInt>(__attribute__((unused)) UInt x) { return (1); } template <> inline Real pow<0, Real>(__attribute__((unused)) Real x) { return (1.); } /* ------------------------------------------------------------------------ */ template Real NewtonRaphson::solve(const Functor & funct, Real x_0) { Real x = x_0; Real f_x = funct.f(x); UInt iter = 0; while (std::abs(f_x) > this->tolerance && iter < this->max_iteration) { x -= f_x / funct.f_prime(x); f_x = funct.f(x); iter++; } AKANTU_DEBUG_ASSERT(iter < this->max_iteration, "Newton Raphson (" << funct.name << ") solve did not converge in " << this->max_iteration << " iterations (tolerance: " << this->tolerance << ")"); return x; } } // namespace Math } // namespace akantu diff --git a/src/mesh/element_type_map.hh b/src/mesh/element_type_map.hh index c5722a3c0..187f901c4 100644 --- a/src/mesh/element_type_map.hh +++ b/src/mesh/element_type_map.hh @@ -1,492 +1,491 @@ /** * @file element_type_map.hh * * @author Lucas Frerot * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Thu Mar 11 2021 * * @brief storage class by element type * * * @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 "aka_array.hh" #include "aka_named_argument.hh" #include "element.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_TYPE_MAP_HH_ #define AKANTU_ELEMENT_TYPE_MAP_HH_ namespace akantu { class FEEngine; } // namespace akantu namespace akantu { namespace { DECLARE_NAMED_ARGUMENT(all_ghost_types); DECLARE_NAMED_ARGUMENT(default_value); DECLARE_NAMED_ARGUMENT(element_kind); DECLARE_NAMED_ARGUMENT(ghost_type); DECLARE_NAMED_ARGUMENT(nb_component); DECLARE_NAMED_ARGUMENT(nb_component_functor); DECLARE_NAMED_ARGUMENT(with_nb_element); DECLARE_NAMED_ARGUMENT(with_nb_nodes_per_element); DECLARE_NAMED_ARGUMENT(spatial_dimension); DECLARE_NAMED_ARGUMENT(do_not_default); DECLARE_NAMED_ARGUMENT(element_filter); } // namespace template class ElementTypeMap; /* -------------------------------------------------------------------------- */ /* ElementTypeMapBase */ /* -------------------------------------------------------------------------- */ /// Common non templated base class for the ElementTypeMap class class ElementTypeMapBase { public: virtual ~ElementTypeMapBase() = default; }; /* -------------------------------------------------------------------------- */ /* ElementTypeMap */ /* -------------------------------------------------------------------------- */ template class ElementTypeMap : public ElementTypeMapBase { public: ElementTypeMap(); ~ElementTypeMap() override; inline static std::string printType(const SupportType & type, GhostType ghost_type); /*! Tests whether a type is present in the object * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return true if the type is present. */ inline bool exists(const SupportType & type, GhostType ghost_type = _not_ghost) const; /*! get the stored data corresponding to a type * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ inline const Stored & operator()(const SupportType & type, GhostType ghost_type = _not_ghost) const; /*! get the stored data corresponding to a type * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ inline Stored & operator()(const SupportType & type, GhostType ghost_type = _not_ghost); /*! insert data of a new type (not yet present) into the map. THIS METHOD IS * NOT ARRAY SAFE, when using ElementTypeMapArray, use setArray instead * @param data to insert * @param type type of data (if this type is already present in the map, * an exception is thrown). * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ template inline Stored & operator()(U && insertee, const SupportType & type, GhostType ghost_type = _not_ghost); public: /// print helper virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ /*! iterator allows to iterate over type-data pairs of the map. The interface * expects the SupportType to be ElementType. */ using DataMap = std::map; /// helper class to use in range for constructions - class type_iterator - : private std::iterator { + class type_iterator { public: using value_type = const SupportType; using pointer = const SupportType *; using reference = const SupportType &; protected: using DataMapIterator = typename ElementTypeMap::DataMap::const_iterator; public: type_iterator(DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek); type_iterator(const type_iterator & it); type_iterator() = default; inline reference operator*(); inline reference operator*() const; inline type_iterator & operator++(); type_iterator operator++(int); inline bool operator==(const type_iterator & other) const; inline bool operator!=(const type_iterator & other) const; type_iterator & operator=(const type_iterator & it); private: DataMapIterator list_begin; DataMapIterator list_end; UInt dim; ElementKind kind; }; /// helper class to use in range for constructions class ElementTypesIteratorHelper { public: using Container = ElementTypeMap; using iterator = typename Container::type_iterator; ElementTypesIteratorHelper(const Container & container, UInt dim, GhostType ghost_type, ElementKind kind) : container(std::cref(container)), dim(dim), ghost_type(ghost_type), kind(kind) {} template ElementTypesIteratorHelper(const Container & container, use_named_args_t /*unused*/, pack &&... _pack) : ElementTypesIteratorHelper( container, OPTIONAL_NAMED_ARG(spatial_dimension, _all_dimensions), OPTIONAL_NAMED_ARG(ghost_type, _not_ghost), OPTIONAL_NAMED_ARG(element_kind, _ek_not_defined)) {} ElementTypesIteratorHelper(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(ElementTypesIteratorHelper &&) noexcept = default; iterator begin(); iterator end(); private: std::reference_wrapper container; UInt dim; GhostType ghost_type; ElementKind kind; }; private: ElementTypesIteratorHelper elementTypesImpl(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; template ElementTypesIteratorHelper elementTypesImpl(const use_named_args_t & /*unused*/, pack &&... _pack) const; public: /*! * \param _pack * \parblock * represent optional parameters: * \li \c _spatial_dimension filter for elements of given spatial * dimension * \li \c _ghost_type filter for a certain ghost_type * \li \c _element_kind filter for elements of given kind * \endparblock */ template std::enable_if_t::value, ElementTypesIteratorHelper> elementTypes(pack &&... _pack) const { return elementTypesImpl(use_named_args, std::forward(_pack)...); } template std::enable_if_t::value, ElementTypesIteratorHelper> elementTypes(pack &&... _pack) const { return elementTypesImpl(std::forward(_pack)...); } /*! Get an iterator to the beginning of a subset datamap. This method expects * the SupportType to be ElementType. * @param dim optional: iterate over data of dimension dim (e.g. when * iterating over (surface) facets of a 3D mesh, dim would be 2). * by default, all dimensions are considered. * @param ghost_type optional: by default, the data map for non-ghost * elements is iterated over. * @param kind optional: the kind of element to search for (see * aka_common.hh), by default all kinds are considered * @return an iterator to the first stored data matching the filters * or an iterator to the end of the map if none match*/ [[deprecated("Use elementTypes instead")]] inline type_iterator firstType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; /*! Get an iterator to the end of a subset datamap. This method expects * the SupportType to be ElementType. * @param dim optional: iterate over data of dimension dim (e.g. when * iterating over (surface) facets of a 3D mesh, dim would be 2). * by default, all dimensions are considered. * @param ghost_type optional: by default, the data map for non-ghost * elements is iterated over. * @param kind optional: the kind of element to search for (see * aka_common.hh), by default all kinds are considered * @return an iterator to the last stored data matching the filters * or an iterator to the end of the map if none match */ [[deprecated("Use elementTypes instead")]] inline type_iterator lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; /*! Direct access to the underlying data map. for internal use by daughter * classes only * @param ghost_type whether to return the data map or the ghost_data map * @return the raw map */ inline DataMap & getData(GhostType ghost_type); /*! Direct access to the underlying data map. for internal use by daughter * classes only * @param ghost_type whether to return the data map or the ghost_data map * @return the raw map */ inline const DataMap & getData(GhostType ghost_type) const; /* ------------------------------------------------------------------------ */ protected: DataMap data; DataMap ghost_data; }; /* -------------------------------------------------------------------------- */ /* Some typedefs */ /* -------------------------------------------------------------------------- */ template class ElementTypeMapArray : public ElementTypeMap>, SupportType> { public: using value_type = T; using array_type = Array; protected: using parent = ElementTypeMap>, SupportType>; using DataMap = typename parent::DataMap; public: using type_iterator = typename parent::type_iterator; /// standard assigment (copy) operator auto operator=(const ElementTypeMapArray & other) -> ElementTypeMapArray &; ElementTypeMapArray(const ElementTypeMapArray & other); /// explicit copy void copy(const ElementTypeMapArray & other); /*! Constructor * @param id optional: identifier (string) * @param parent_id optional: parent identifier. for organizational purposes * only */ ElementTypeMapArray(const ID & id = "by_element_type_array", const ID & parent_id = "no_parent") : parent(), id(parent_id + ":" + id), name(id){}; /*! allocate memory for a new array * @param size number of tuples of the new array * @param nb_component tuple size * @param type the type under which the array is indexed in the map * @param ghost_type whether to add the field to the data map or the * ghost_data map * @param default_value the default value to use to fill the array * @return a reference to the allocated array */ inline Array & alloc(UInt size, UInt nb_component, const SupportType & type, GhostType ghost_type, const T & default_value = T()); /*! allocate memory for a new array in both the data and the ghost_data map * @param size number of tuples of the new array * @param nb_component tuple size * @param type the type under which the array is indexed in the map * @param default_value the default value to use to fill the array */ inline void alloc(UInt size, UInt nb_component, const SupportType & type, const T & default_value = T()); /* get a reference to the array of certain type * @param type data filed under type is returned * @param ghost_type optional: by default the non-ghost map is searched * @return a reference to the array */ inline const Array & operator()(const SupportType & type, GhostType ghost_type = _not_ghost) const; /// access the data of an element, this combine the map and array accessor inline const T & operator()(const Element & element, UInt component = 0) const; /// access the data of an element, this combine the map and array accessor inline T & operator()(const Element & element, UInt component = 0); /// access the data of an element, this combine the map and array accessor inline decltype(auto) get(const Element & element); inline decltype(auto) get(const Element & element) const; /* get a reference to the array of certain type * @param type data filed under type is returned * @param ghost_type optional: by default the non-ghost map is searched * @return a const reference to the array */ inline Array & operator()(const SupportType & type, GhostType ghost_type = _not_ghost); /*! insert data of a new type (not yet present) into the map. * @param type type of data (if this type is already present in the map, * an exception is thrown). * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @param vect the vector to include into the map * @return stored data corresponding to type. */ inline void setArray(const SupportType & type, GhostType ghost_type, const Array & vect); /*! frees all memory related to the data*/ inline void free(); inline void clear(); inline bool empty() const __attribute__((warn_unused_result)); /*! set all values in the ElementTypeMap to zero*/ inline void zero() { this->set(T()); } /*! set all values in the ElementTypeMap to value */ template inline void set(const ST & value); /*! deletes and reorders entries in the stored arrays * @param new_numbering a ElementTypeMapArray of new indices. UInt(-1) * indicates * deleted entries. */ inline void onElementsRemoved(const ElementTypeMapArray & new_numbering); /// text output helper void printself(std::ostream & stream, int indent = 0) const override; /*! set the id * @param id the new name */ inline void setID(const ID & id) { this->id = id; } /// return the id inline auto getID() const -> ID { return this->id; } ElementTypeMap getNbComponents(UInt dim = _all_dimensions, GhostType requested_ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const { ElementTypeMap nb_components; bool all_ghost_types = requested_ghost_type == _casper; for (auto ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) { continue; } for (auto & type : this->elementTypes(dim, ghost_type, kind)) { UInt nb_comp = (*this)(type, ghost_type).getNbComponent(); nb_components(type, ghost_type) = nb_comp; } } return nb_components; } /* ------------------------------------------------------------------------ */ /* more evolved allocators */ /* ------------------------------------------------------------------------ */ public: /// initialize the arrays in accordance to a functor template void initialize(const Func & f, const T & default_value, bool do_not_default); /// initialize with sizes and number of components in accordance of a mesh /// content template void initialize(const Mesh & mesh, pack &&... _pack); /// initialize with sizes and number of components in accordance of a fe /// engine content (aka integration points) template void initialize(const FEEngine & fe_engine, pack &&... _pack); /* ------------------------------------------------------------------------ */ /* Accesssors */ /* ------------------------------------------------------------------------ */ public: /// get the name of the internal field AKANTU_GET_MACRO(Name, name, ID); /** * get the size of the ElementTypeMapArray * @param[in] _pack * \parblock * optional arguments can be any of: * \li \c _spatial_dimension the dimension to consider (default: * _all_dimensions) * \li \c _ghost_type (default: _not_ghost) * \li \c _element_kind (default: _ek_not_defined) * \li \c _all_ghost_types (default: false) * \endparblock **/ template UInt size(pack &&... _pack) const; bool isNodal() const { return is_nodal; } void isNodal(bool is_nodal) { this->is_nodal = is_nodal; } private: UInt sizeImpl(UInt spatial_dimension, GhostType ghost_type, ElementKind kind) const; private: ID id; protected: /// name of the element type map: e.g. connectivity, grad_u ID name; /// Is the data stored by node of the element bool is_nodal{false}; }; /// to store data Array by element type using ElementTypeMapReal = ElementTypeMapArray; /// to store data Array by element type using ElementTypeMapInt = ElementTypeMapArray; /// to store data Array by element type using ElementTypeMapUInt = ElementTypeMapArray; } // namespace akantu #endif /* AKANTU_ELEMENT_TYPE_MAP_HH_ */ diff --git a/src/mesh_utils/cohesive_element_inserter.hh b/src/mesh_utils/cohesive_element_inserter.hh index 43bbb6624..3ca856ceb 100644 --- a/src/mesh_utils/cohesive_element_inserter.hh +++ b/src/mesh_utils/cohesive_element_inserter.hh @@ -1,175 +1,181 @@ /** * @file cohesive_element_inserter.hh * * @author Fabian Barras * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Dec 04 2013 * @date last modification: Tue Jul 21 2020 * * @brief Cohesive element inserter * * * @section LICENSE * * Copyright (©) 2014-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 "data_accessor.hh" #include "mesh_utils.hh" #include "parsable.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_COHESIVE_ELEMENT_INSERTER_HH_ #define AKANTU_COHESIVE_ELEMENT_INSERTER_HH_ namespace akantu { class GlobalIdsUpdater; class FacetSynchronizer; class SolidMechanicsModeslCohesivel; } // namespace akantu namespace akantu { class CohesiveElementInserter : public DataAccessor, public Parsable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: CohesiveElementInserter(Mesh & mesh, const ID & id = "cohesive_element_inserter"); ~CohesiveElementInserter() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// set range limitation for intrinsic cohesive element insertion void setLimit(SpatialDirection axis, Real first_limit, Real second_limit); /// insert intrinsic cohesive elements in a predefined range auto insertIntrinsicElements() -> UInt; /// insert extrinsic cohesive elements (returns the number of new /// cohesive elements) UInt insertElements(bool only_double_facets = false); /// limit check facets to match given insertion limits void limitCheckFacets(); + /// limit insertion to the surfaces + void addPhysicalSurface(const ID & surface_name); + + /// limit insertion to the volumes + void addPhysicalVolume(const ID & surface_name); + protected: void parseSection(const ParserSection & section) override; protected: /// internal version of limitCheckFacets void limitCheckFacets(ElementTypeMapArray & check_facets); /// update facet insertion arrays after facets doubling void updateInsertionFacets(); /// functions for parallel communications 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; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO_NOT_CONST(InsertionFacetsByElement, insertion_facets, ElementTypeMapArray &); AKANTU_GET_MACRO(InsertionFacetsByElement, insertion_facets, const ElementTypeMapArray &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(InsertionFacets, insertion_facets, bool); AKANTU_GET_MACRO(CheckFacets, check_facets, const ElementTypeMapArray &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(CheckFacets, check_facets, bool); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(CheckFacets, check_facets, bool); AKANTU_GET_MACRO(MeshFacets, mesh_facets, const Mesh &); AKANTU_GET_MACRO_NOT_CONST(MeshFacets, mesh_facets, Mesh &); AKANTU_SET_MACRO(IsExtrinsic, is_extrinsic, bool); public: friend class SolidMechanicsModelCohesive; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// object id ID id; /// main mesh where to insert cohesive elements Mesh & mesh; /// mesh containing facets Mesh & mesh_facets; /// list of facets where to insert elements ElementTypeMapArray insertion_facets; /// limits for element insertion Matrix insertion_limits; /// list of groups to consider for insertion, ignored if empty std::set physical_surfaces; /// list of groups in between which an inside which element are insterted std::set physical_zones; /// vector containing facets in which extrinsic cohesive elements can be /// inserted ElementTypeMapArray check_facets; /// global connectivity ids updater std::unique_ptr global_ids_updater; /// is this inserter used in extrinsic bool is_extrinsic{false}; }; class CohesiveNewNodesEvent : public NewNodesEvent { public: CohesiveNewNodesEvent(const std::string & origin) : NewNodesEvent(origin) {} ~CohesiveNewNodesEvent() override = default; AKANTU_GET_MACRO_NOT_CONST(OldNodesList, old_nodes, Array &); AKANTU_GET_MACRO(OldNodesList, old_nodes, const Array &); private: Array old_nodes; }; } // namespace akantu #include "cohesive_element_inserter_inline_impl.hh" #endif /* AKANTU_COHESIVE_ELEMENT_INSERTER_HH_ */ diff --git a/src/mesh_utils/cohesive_element_inserter_inline_impl.hh b/src/mesh_utils/cohesive_element_inserter_inline_impl.hh index a3afb8ae5..f671cc812 100644 --- a/src/mesh_utils/cohesive_element_inserter_inline_impl.hh +++ b/src/mesh_utils/cohesive_element_inserter_inline_impl.hh @@ -1,91 +1,103 @@ /** * @file cohesive_element_inserter_inline_impl.hh * * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Oct 13 2017 * @date last modification: Wed Nov 11 2020 * * @brief Cohesive element inserter inline functions * * * @section LICENSE * * Copyright (©) 2016-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 "cohesive_element_inserter.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_COHESIVE_ELEMENT_INSERTER_INLINE_IMPL_HH_ #define AKANTU_COHESIVE_ELEMENT_INSERTER_INLINE_IMPL_HH_ namespace akantu { +/* -------------------------------------------------------------------------- */ +inline void +CohesiveElementInserter::addPhysicalSurface(const ID & surface_name) { + physical_surfaces.insert(surface_name); +} + +/* -------------------------------------------------------------------------- */ +inline void +CohesiveElementInserter::addPhysicalVolume(const ID & surface_name) { + physical_zones.insert(surface_name); +} + /* -------------------------------------------------------------------------- */ inline UInt CohesiveElementInserter::getNbData(const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; if (tag == SynchronizationTag::_ce_groups) { size = elements.size() * sizeof(bool); } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void CohesiveElementInserter::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); if (tag == SynchronizationTag::_ce_groups) { for (const auto & el : elements) { const bool & data = insertion_facets(el); buffer << data; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void CohesiveElementInserter::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); if (tag == SynchronizationTag::_ce_groups) { for (const auto & el : elements) { bool & data = insertion_facets(el); buffer >> data; } } AKANTU_DEBUG_OUT(); } } // namespace akantu #endif /* AKANTU_COHESIVE_ELEMENT_INSERTER_INLINE_IMPL_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc index 816c1e240..e96efb926 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc @@ -1,569 +1,559 @@ /** * @file material_cohesive.cc * * @author Mauro Corrado * @author Nicolas Richart * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @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_cohesive.hh" #include "aka_random_generator.hh" #include "dof_synchronizer.hh" #include "fe_engine_template.hh" #include "integrator_gauss.hh" #include "shape_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ MaterialCohesive::MaterialCohesive(SolidMechanicsModel & model, const ID & id) : Material(model, id), facet_filter("facet_filter", id), fem_cohesive( model.getFEEngineClass("CohesiveFEEngine")), reversible_energy("reversible_energy", *this), total_energy("total_energy", *this), opening("opening", *this), tractions("tractions", *this), contact_tractions("contact_tractions", *this), contact_opening("contact_opening", *this), delta_max("delta max", *this), use_previous_delta_max(false), use_previous_opening(false), damage("damage", *this), sigma_c("sigma_c", *this), normal(0, spatial_dimension, "normal") { AKANTU_DEBUG_IN(); this->model = dynamic_cast(&model); this->registerParam("sigma_c", sigma_c, _pat_parsable | _pat_readable, "Critical stress"); this->registerParam("delta_c", delta_c, Real(0.), _pat_parsable | _pat_readable, "Critical displacement"); this->element_filter.initialize(this->model->getMesh(), _spatial_dimension = spatial_dimension, _element_kind = _ek_cohesive); // this->model->getMesh().initElementTypeMapArray( // this->element_filter, 1, spatial_dimension, false, _ek_cohesive); if (this->model->getIsExtrinsic()) { this->facet_filter.initialize(this->model->getMeshFacets(), _spatial_dimension = spatial_dimension - 1, _element_kind = _ek_regular); } // this->model->getMeshFacets().initElementTypeMapArray(facet_filter, 1, // spatial_dimension - // 1); this->reversible_energy.initialize(1); this->total_energy.initialize(1); this->tractions.initialize(spatial_dimension); this->tractions.initializeHistory(); this->contact_tractions.initialize(spatial_dimension); this->contact_opening.initialize(spatial_dimension); this->opening.initialize(spatial_dimension); this->opening.initializeHistory(); this->delta_max.initialize(1); this->damage.initialize(1); if (this->model->getIsExtrinsic()) { this->sigma_c.initialize(1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MaterialCohesive::~MaterialCohesive() = default; /* -------------------------------------------------------------------------- */ void MaterialCohesive::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); if (this->use_previous_delta_max) { this->delta_max.initializeHistory(); } if (this->use_previous_opening) { this->opening.initializeHistory(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); -#if defined(AKANTU_DEBUG_TOOLS) - debug::element_manager.printData(debug::_dm_material_cohesive, - "Cohesive Tractions", tractions); -#endif - auto & internal_force = const_cast &>(model->getInternalForce()); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { auto & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) { continue; } const auto & shapes = fem_cohesive.getShapes(type, ghost_type); auto & traction = tractions(type, ghost_type); auto & contact_traction = contact_tractions(type, ghost_type); UInt size_of_shapes = shapes.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem_cohesive.getNbIntegrationPoints(type, ghost_type); /// compute @f$t_i N_a@f$ auto * traction_cpy = new Array(nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes); auto traction_it = traction.begin(spatial_dimension, 1); auto contact_traction_it = contact_traction.begin(spatial_dimension, 1); auto shapes_filtered_begin = shapes.begin(1, size_of_shapes); auto traction_cpy_it = traction_cpy->begin(spatial_dimension, size_of_shapes); Matrix traction_tmp(spatial_dimension, 1); for (UInt el = 0; el < nb_element; ++el) { UInt current_quad = elem_filter(el) * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++traction_it, ++contact_traction_it, ++current_quad, ++traction_cpy_it) { const Matrix & shapes_filtered = shapes_filtered_begin[current_quad]; traction_tmp.copy(*traction_it); traction_tmp += *contact_traction_it; traction_cpy_it->mul(traction_tmp, shapes_filtered); } } /** * compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t * \mathbf{t}_q \overline w_q J_q@f$ */ auto * partial_int_t_N = new Array( nb_element, spatial_dimension * size_of_shapes, "int_t_N"); fem_cohesive.integrate(*traction_cpy, *partial_int_t_N, spatial_dimension * size_of_shapes, type, ghost_type, elem_filter); delete traction_cpy; auto * int_t_N = new Array( nb_element, 2 * spatial_dimension * size_of_shapes, "int_t_N"); Real * int_t_N_val = int_t_N->storage(); Real * partial_int_t_N_val = partial_int_t_N->storage(); for (UInt el = 0; el < nb_element; ++el) { std::copy_n(partial_int_t_N_val, size_of_shapes * spatial_dimension, int_t_N_val); std::copy_n(partial_int_t_N_val, size_of_shapes * spatial_dimension, int_t_N_val + size_of_shapes * spatial_dimension); for (UInt n = 0; n < size_of_shapes * spatial_dimension; ++n) { int_t_N_val[n] *= -1.; } int_t_N_val += nb_nodes_per_element * spatial_dimension; partial_int_t_N_val += size_of_shapes * spatial_dimension; } delete partial_int_t_N; /// assemble model->getDOFManager().assembleElementalArrayLocalArray( *int_t_N, internal_force, type, ghost_type, 1, elem_filter); delete int_t_N; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { UInt nb_quadrature_points = fem_cohesive.getNbIntegrationPoints(type, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Array & shapes = fem_cohesive.getShapes(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0U) { continue; } UInt size_of_shapes = shapes.getNbComponent(); auto * shapes_filtered = new Array(nb_element * nb_quadrature_points, size_of_shapes, "filtered shapes"); Real * shapes_filtered_val = shapes_filtered->storage(); UInt * elem_filter_val = elem_filter.storage(); for (UInt el = 0; el < nb_element; ++el) { auto * shapes_val = shapes.storage() + elem_filter_val[el] * size_of_shapes * nb_quadrature_points; memcpy(shapes_filtered_val, shapes_val, size_of_shapes * nb_quadrature_points * sizeof(Real)); shapes_filtered_val += size_of_shapes * nb_quadrature_points; } Matrix A(spatial_dimension * size_of_shapes, spatial_dimension * nb_nodes_per_element); for (UInt i = 0; i < spatial_dimension * size_of_shapes; ++i) { A(i, i) = 1; A(i, i + spatial_dimension * size_of_shapes) = -1; } /// get the tangent matrix @f$\frac{\partial{(t/\delta)}}{\partial{\delta}} /// @f$ auto * tangent_stiffness_matrix = new Array( nb_element * nb_quadrature_points, spatial_dimension * spatial_dimension, "tangent_stiffness_matrix"); // Array * normal = new Array(nb_element * // nb_quadrature_points, spatial_dimension, "normal"); normal.resize(nb_quadrature_points); computeNormal(model->getCurrentPosition(), normal, type, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ // computeOpening(model->getDisplacement(), opening(type, ghost_type), type, // ghost_type); tangent_stiffness_matrix->zero(); computeTangentTraction(type, *tangent_stiffness_matrix, normal, ghost_type); // delete normal; UInt size_at_nt_d_n_a = spatial_dimension * nb_nodes_per_element * spatial_dimension * nb_nodes_per_element; auto * at_nt_d_n_a = new Array(nb_element * nb_quadrature_points, size_at_nt_d_n_a, "A^t*N^t*D*N*A"); Array::iterator> shapes_filt_it = shapes_filtered->begin(size_of_shapes); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(spatial_dimension, spatial_dimension); Array::matrix_iterator At_Nt_D_N_A_it = at_nt_d_n_a->begin(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Array::matrix_iterator At_Nt_D_N_A_end = at_nt_d_n_a->end(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Matrix N(spatial_dimension, spatial_dimension * size_of_shapes); Matrix N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); Matrix D_N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); for (; At_Nt_D_N_A_it != At_Nt_D_N_A_end; ++At_Nt_D_N_A_it, ++D_it, ++shapes_filt_it) { N.zero(); /** * store the shapes in voigt notations matrix @f$\mathbf{N} = * \begin{array}{cccccc} N_0(\xi) & 0 & N_1(\xi) &0 & N_2(\xi) & 0 \\ * 0 & * N_0(\xi)& 0 &N_1(\xi)& 0 & N_2(\xi) \end{array} @f$ **/ for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt n = 0; n < size_of_shapes; ++n) { N(i, i + spatial_dimension * n) = (*shapes_filt_it)(n); } } /** * compute stiffness matrix @f$ \mathbf{K} = \delta \mathbf{U}^T * \int_{\Gamma_c} {\mathbf{P}^t \frac{\partial{\mathbf{t}}} *{\partial{\delta}} * \mathbf{P} d\Gamma \Delta \mathbf{U}} @f$ **/ N_A.mul(N, A); D_N_A.mul(*D_it, N_A); (*At_Nt_D_N_A_it).mul(D_N_A, N_A); } delete tangent_stiffness_matrix; delete shapes_filtered; auto * K_e = new Array(nb_element, size_at_nt_d_n_a, "K_e"); fem_cohesive.integrate(*at_nt_d_n_a, *K_e, size_at_nt_d_n_a, type, ghost_type, elem_filter); delete at_nt_d_n_a; model->getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _unsymmetric, elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- * * Compute traction from displacements * * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialCohesive::computeTraction(GhostType ghost_type) { AKANTU_DEBUG_IN(); -#if defined(AKANTU_DEBUG_TOOLS) - debug::element_manager.printData(debug::_dm_material_cohesive, - "Cohesive Openings", opening); -#endif - for (const auto & type : element_filter.elementTypes( spatial_dimension, ghost_type, _ek_cohesive)) { Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) { continue; } UInt nb_quadrature_points = nb_element * fem_cohesive.getNbIntegrationPoints(type, ghost_type); normal.resize(nb_quadrature_points); /// compute normals @f$\mathbf{n}@f$ computeNormal(model->getCurrentPosition(), normal, type, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ computeOpening(model->getDisplacement(), opening(type, ghost_type), type, ghost_type); /// compute traction @f$\mathbf{t}@f$ computeTraction(normal, type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeNormal(const Array & position, Array & normal, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem_cohesive = this->model->getFEEngineClass("CohesiveFEEngine"); normal.zero(); #define COMPUTE_NORMAL(type) \ fem_cohesive.getShapeFunctions() \ .computeNormalsOnIntegrationPoints( \ position, normal, ghost_type, element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMAL); #undef COMPUTE_NORMAL AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeOpening(const Array & displacement, Array & opening, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem_cohesive = this->model->getFEEngineClass("CohesiveFEEngine"); #define COMPUTE_OPENING(type) \ fem_cohesive.getShapeFunctions() \ .interpolateOnIntegrationPoints( \ displacement, opening, spatial_dimension, ghost_type, \ element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_OPENING); #undef COMPUTE_OPENING AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::updateEnergies(ElementType type) { AKANTU_DEBUG_IN(); if (Mesh::getKind(type) != _ek_cohesive) { return; } Vector b(spatial_dimension); Vector h(spatial_dimension); auto erev = reversible_energy(type).begin(); auto etot = total_energy(type).begin(); auto traction_it = tractions(type).begin(spatial_dimension); auto traction_old_it = tractions.previous(type).begin(spatial_dimension); auto opening_it = opening(type).begin(spatial_dimension); auto opening_old_it = opening.previous(type).begin(spatial_dimension); auto traction_end = tractions(type).end(spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++traction_old_it, ++opening_it, ++opening_old_it, ++erev, ++etot) { /// trapezoidal integration b = *opening_it; b -= *opening_old_it; h = *traction_old_it; h += *traction_it; *etot += .5 * b.dot(h); *erev = .5 * traction_it->dot(*opening_it); } /// update old values AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getReversibleEnergy() { AKANTU_DEBUG_IN(); Real erev = 0.; /// integrate reversible energy for each type of elements for (const auto & type : element_filter.elementTypes( spatial_dimension, _not_ghost, _ek_cohesive)) { erev += fem_cohesive.integrate(reversible_energy(type, _not_ghost), type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return erev; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getDissipatedEnergy() { AKANTU_DEBUG_IN(); Real edis = 0.; /// integrate dissipated energy for each type of elements for (const auto & type : element_filter.elementTypes( spatial_dimension, _not_ghost, _ek_cohesive)) { Array dissipated_energy(total_energy(type, _not_ghost)); dissipated_energy -= reversible_energy(type, _not_ghost); edis += fem_cohesive.integrate(dissipated_energy, type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return edis; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getContactEnergy() { AKANTU_DEBUG_IN(); Real econ = 0.; /// integrate contact energy for each type of elements for (const auto & type : element_filter.elementTypes( spatial_dimension, _not_ghost, _ek_cohesive)) { auto & el_filter = element_filter(type, _not_ghost); UInt nb_quad_per_el = fem_cohesive.getNbIntegrationPoints(type, _not_ghost); UInt nb_quad_points = el_filter.size() * nb_quad_per_el; Array contact_energy(nb_quad_points); auto contact_traction_it = contact_tractions(type, _not_ghost).begin(spatial_dimension); auto contact_opening_it = contact_opening(type, _not_ghost).begin(spatial_dimension); /// loop on each quadrature point for (UInt q = 0; q < nb_quad_points; ++contact_traction_it, ++contact_opening_it, ++q) { contact_energy(q) = .5 * contact_traction_it->dot(*contact_opening_it); } econ += fem_cohesive.integrate(contact_energy, type, _not_ghost, el_filter); } AKANTU_DEBUG_OUT(); return econ; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getEnergy(const std::string & type) { if (type == "reversible") { return getReversibleEnergy(); } if (type == "dissipated") { return getDissipatedEnergy(); } if (type == "cohesive contact") { return getContactEnergy(); } return 0.; } /* -------------------------------------------------------------------------- */ } // 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_ */