diff --git a/examples/c++/solid_mechanics_cohesive_model/my_first_test/material.dat b/examples/c++/solid_mechanics_cohesive_model/my_first_test/material.dat index 1f6fbb030..2f627470c 100644 --- a/examples/c++/solid_mechanics_cohesive_model/my_first_test/material.dat +++ b/examples/c++/solid_mechanics_cohesive_model/my_first_test/material.dat @@ -1,13 +1,13 @@ material elastic [ name = steel rho = 1 # density - E = 1 # young's modulus - nu = 0.1 # poisson's ratio + E = 3e10 # young's modulus + nu = 0. # poisson's ratio ] -material cohesive_damage [ +material cohesive_damage_intrinsic [ name = cohesive - sigma_c = 1.0 - k = 1.0 - G_c = 1.0 + sigma_c = 3e6 + k = 75e9 + G_c = 120 ] diff --git a/examples/c++/solid_mechanics_cohesive_model/my_first_test/my_first_test.py b/examples/c++/solid_mechanics_cohesive_model/my_first_test/my_first_test.py new file mode 100644 index 000000000..eda628221 --- /dev/null +++ b/examples/c++/solid_mechanics_cohesive_model/my_first_test/my_first_test.py @@ -0,0 +1,244 @@ +#!/usr/bin/env python3 +__copyright__ = ( + "Copyright (©) 2019-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)" + "Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)" +) +__license__ = "LGPLv3" + + +import sys +sys.path.insert(0,"/home/ble/buildAkantu/buildDebug/python") +print(sys.path) + +import akantu as aka +import numpy as np + +import scipy.optimize + +import matplotlib.pyplot as plt + +class degradation_function: + # g(d) + def val(self,d): + return + # dg_dd(d) + def deriv(self,d): + return + # 1+g(d) + def augmented_stiffness(self,d): + return + def augmented_compliance(self,d): + # 1/(1+g(d)) + return + def augmented_compliance_deriv(self,d): + # d()1/(1+g(d))_dd + return + +class degradation_function_linear_czm(degradation_function): + def val(self,d): + return 1./d-1. + def deriv(self,d): + return -1./d**2 + def augmented_stiffness(self,d): + return 1./d + def augmented_compliance(self,d): + return d + def augmented_compliance_deriv(self,d): + return 1. + +class softening_function: + def val(self,d): + return + def deriv(self,d): + return + +class softening_function_linear_czm: + def __init__(self,k,tc,Gc): + wc = 2.*Gc/tc + self.a = k*wc/tc + def val(self,d): + return d/(d+self.a*(1.-d)) + def deriv(self,d): + return self.a/(d+self.a*(1.-d))**2 + +class damage_updater: + def __init__(self,Gc, sigc, k, g, h): + self.Gc = Gc + self.sigc = sigc + self.k = k + self.g = g + self.h = h + def update_d(self,d_previous,lda): + dz = d_previous.flatten() + shape_d = d_previous.shape + lda_lda_on_k = (lda[:,0]**2 + lda[:,1]**2)/self.k + + def F(d): + return -0.5*np.dot(lda_lda_on_k,self.g.augmented_compliance(d))+self.Gc*np.sum(self.h.val(d)) + def dF(d): + return -0.5*self.g.augmented_compliance_deriv(d)*lda_lda_on_k+self.Gc*self.h.deriv(d) + + #dd = np.linspace(0.,1.,100) + #D = dz.copy() + #FF = [] + #for i in range(len(dd)): + # D[0] = dd[i] + # FF.append(F(D)) + #plt.plot(dd,FF) + #plt.show() + #print("dz = ",dz) + #print("lda_lda_on_k = ",lda_lda_on_k) + #print("F1 = ",-0.5*np.dot(lda_lda_on_k,self.g.augmented_compliance(dz))) + #print("F2 = ",self.Gc*np.sum(self.h.val(dz))) + #print("F(dz) = ",F(dz)) + #print("self.k = ",self.k) + + boundd = scipy.optimize.Bounds(dz, np.ones(len(dz))) + options = {'verbose':1} + resd = scipy.optimize.minimize(F, dz, bounds=boundd + , jac =dF + #, method = 'trust-constr' + #, options = options + #, tol = 1e-10 + ) + if not resd.success: + print("resd.message = ",resd.message) + #print("niterd = ",resd.nit) + return np.resize(resd.x,shape_d) + +def energy(model): + return model.getEnergy("potential") + model.getEnergy("dissipated") + +def set_dumpers(model): + model.setBaseName("cohesive") + 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") + + +def solve(material_file, mesh_file, increment): + aka.parseInput(material_file) + spatial_dimension = 2 + + # ------------------------------------------------------------------------- + # Initialization + # ------------------------------------------------------------------------- + mesh = aka.Mesh(spatial_dimension) + mesh.read(mesh_file) + + model = aka.SolidMechanicsModelCohesive(mesh) + model.getElementInserter().setLimit(aka._x, -0.1, 0.1); + + model.initFull(_analysis_method=aka._static, _is_extrinsic=False) + + set_dumpers(model) + + L = 0.4 + E = model.getMaterial("steel").getReal("E") + Gc = model.getMaterial("cohesive").getReal("G_c") + tc = model.getMaterial("cohesive").getReal("sigma_c") + wc = 2.*Gc/tc + k = model.getMaterial("cohesive").getReal("k") + g = degradation_function_linear_czm() + h = softening_function_linear_czm(k,tc,Gc) + updater = damage_updater(Gc,tc,k,g,h) + + imposed_disp = [0.] + reaction_force = [0.] + U = 0. + F = 0. + + # ------------------------------------------------------------------------- + # Boundary conditions + # ------------------------------------------------------------------------- + model.applyBC(aka.FixedValue(0.0, aka._x), "left") + model.applyBC(aka.FixedValue(0.0, aka._y), "point") + model.applyBC(aka.IncrementValue(increment, aka._x), "right") + U+=increment + imposed_disp.append(U) + + model.getExternalForce()[:] = 0 + + solver = model.getNonLinearSolver("static") + solver.set("max_iterations", 2) + solver.set("threshold", 1e-10) + solver.set("convergence_type", aka.SolveConvergenceCriteria.residual) + + model.getNewSolver("linear_static", aka.TimeStepSolverType.static, + aka.NonLinearSolverType.linear) + model.setIntegrationScheme("linear_static", "displacement", + aka.IntegrationSchemeType.pseudo_time) + model.setIntegrationScheme("linear_static", "lambda", + aka.IntegrationSchemeType.pseudo_time) + + d = model.getMaterial("cohesive").getInternalReal("czm_damage")(aka._cohesive_2d_4) + lda = model.getMaterial("cohesive").getInternalReal("lambda")(aka._cohesive_2d_4) + Fint = model.getInternalForce() + + nodes_right = mesh.getElementGroup("right").getNodeGroup().getNodes() + model.solveStep("linear_static") + F = -np.sum(Fint[nodes_right,0]) + reaction_force.append(F) + + model.dump() + #model.dump("cohesive elements") + + maxsteps = 8 + maxiter = 10 + tol = 1e-6 + + for i in range(0, maxsteps): + print("---- step {0}/{1}".format(i, maxsteps)) + model.applyBC(aka.IncrementValue(increment, aka._x), "right") + U+=increment + imposed_disp.append(U) + it = 0 + d_previous = model.getMaterial("cohesive").getInternalReal("czm_damage")(aka._cohesive_2d_4).copy() + W0 = energy(model) + while it < maxiter: + #aka.debug.setDebugLevel(aka.dblInfo); + model.solveStep("linear_static") + d[:] = updater.update_d(d_previous,lda) + W = energy(model) + err = abs(W-W0)/W + print("iter ",it,", err = ",err) + if err < tol: + break + else: + W0 = W + it+=1 + F = -np.sum(Fint[nodes_right,0]) + reaction_force.append(F) + + model.dump() + #model.dump("cohesive elements") + + plt.plot(imposed_disp,reaction_force,'-o') + plt.plot([0., (tc/E)*L, wc], [0., tc*L, 0.], 'k-o') + plt.show() + W = model.getEnergy("potential"); + D = model.getEnergy("dissipated"); + print("Elastic energy = ",W) + print("Dissipated energy = ",D) + +# ----------------------------------------------------------------------------- +# main +# ----------------------------------------------------------------------------- +def main(): + mesh_file = "triangle.msh" + material_file = "material.dat" + increment = 0.01e-3 + solve(material_file, mesh_file, increment) + + +# ----------------------------------------------------------------------------- +if __name__ == "__main__": + main() diff --git a/packages/cohesive_element.cmake b/packages/cohesive_element.cmake index 8914b039e..a735cb8ce 100644 --- a/packages/cohesive_element.cmake +++ b/packages/cohesive_element.cmake @@ -1,74 +1,79 @@ #=============================================================================== # Copyright (©) 2012-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # This file is part of Akantu # # 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 . # #=============================================================================== package_declare(cohesive_element DEFAULT ON DESCRIPTION "Use cohesive_element package of Akantu" DEPENDS solid_mechanics) package_declare_sources(cohesive_element fe_engine/cohesive_element.hh fe_engine/fe_engine_template_cohesive.cc fe_engine/shape_cohesive.hh fe_engine/shape_cohesive_inline_impl.hh mesh_utils/cohesive_element_inserter.cc mesh_utils/cohesive_element_inserter.hh mesh_utils/cohesive_element_inserter_inline_impl.hh mesh_utils/cohesive_element_inserter_parallel.cc mesh_utils/cohesive_element_inserter_helper.cc mesh_utils/cohesive_element_inserter_helper.hh model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.hh model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.cc model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/cohesive_internal_field.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/cohesive_internal_field_tmpl.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_fatigue.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_fatigue.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_uncoupled.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_uncoupled.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_inline_impl.hh model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh - model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_inline_impl.hh - + model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.cc + model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.hh + model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic_inline_impl.hh + model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.cc + model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.hh + model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic_inline_impl.hh + model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.hh model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_inline_impl.hh model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_parallel.cc ) diff --git a/python/py_material.cc b/python/py_material.cc index 2b56be2d4..49e4946dd 100644 --- a/python/py_material.cc +++ b/python/py_material.cc @@ -1,254 +1,260 @@ /** * Copyright (©) 2019-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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_constitutive_law.hh" /* -------------------------------------------------------------------------- */ #include #include #include #if defined(AKANTU_COHESIVE_ELEMENT) #include #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; void initMaterial() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, initMaterial, ); }; void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(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); } [[nodiscard]] Real getPushWaveSpeed(const Element & element) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(Real, _Material, getPushWaveSpeed, element); } [[nodiscard]] Real getShearWaveSpeed(const Element & element) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(Real, _Material, getShearWaveSpeed, element); } }; /* ------------------------------------------------------------------------ */ template void register_material_classes(py::module & mod, const std::string & name) { py::class_<_Material, Material, PyMaterial<_Material>>( mod, name.c_str(), py::multiple_inheritance()) .def(py::init()); } #if defined(AKANTU_COHESIVE_ELEMENT) // trampoline for the cohesive materials template > * = nullptr> class PyMaterialCohesive : public PyMaterial<_Material> { public: using Parent = PyMaterial<_Material>; /* Inherit the constructors */ using Parent::Parent; void checkInsertion(bool check_only) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, checkInsertion, check_only); } void computeTraction(ElementType el_type, GhostType ghost_type = _not_ghost) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, computeTraction, el_type, ghost_type); } void computeTangentTraction(ElementType el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, computeTangentTraction, el_type, tangent_matrix, ghost_type); } void computeStress(ElementType /*el_type*/, GhostType /*ghost_type*/ = _not_ghost) final {} void computeTangentModuli(ElementType /*el_type*/, Array & /*tangent_matrix*/, GhostType /*ghost_type*/ = _not_ghost) final {} }; template void register_material_cohesive_classes(py::module & mod, const std::string & name) { py::class_<_Material, MaterialCohesive, PyMaterialCohesive<_Material>>( mod, name.c_str(), py::multiple_inheritance()) .def(py::init()); } #endif } // namespace /* -------------------------------------------------------------------------- */ void register_material(py::module & mod) { register_constitutive_law(mod); py::class_, ConstitutiveLaw>(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, Element element) -> Real { return self.getPotentialEnergy(element); }, py::arg("element")) .def("getPotentialEnergy", [](Material & self) -> Real { return self.getPotentialEnergy(); }) .def("initMaterial", &Material::initMaterial) .def("getModel", [](Material & self) -> SolidMechanicsModel & { return self.getModel(); }) .def("getPushWaveSpeed", &Material::getPushWaveSpeed) .def("getShearWaveSpeed", &Material::getShearWaveSpeed); register_material_classes>(mod, "MaterialLinearElastic1D"); register_material_classes>(mod, "MaterialLinearElastic2D"); register_material_classes>(mod, "MaterialLinearElastic3D"); #if defined(AKANTU_COHESIVE_ELEMENT) /* ------------------------------------------------------------------------ */ py::class_>( mod, "MaterialCohesive", py::multiple_inheritance()) .def(py::init()) .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) + .def( + "getDamage", + [](MaterialCohesive & self, ElementType type, GhostType ghost_type) + -> decltype(auto) { return self.getDamage(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"); register_material_cohesive_classes>( mod, "MaterialCohesiveLinearFriction2D"); register_material_cohesive_classes>( mod, "MaterialCohesiveLinearFriction3D"); #endif 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](Int 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); } } // namespace akantu diff --git a/python/py_solid_mechanics_model_cohesive.cc b/python/py_solid_mechanics_model_cohesive.cc index d4d8c8db9..2ff63c793 100644 --- a/python/py_solid_mechanics_model_cohesive.cc +++ b/python/py_solid_mechanics_model_cohesive.cc @@ -1,119 +1,120 @@ /** * Copyright (©) 2020-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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 /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { #define def_function_nocopy(func_name) \ def( \ #func_name, \ - [](SolidMechanicsModel & self) -> decltype(auto) { \ + [](SolidMechanicsModelCohesive & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ - def(#func_name, [](SolidMechanicsModel & self) -> decltype(auto) { \ + def(#func_name, [](SolidMechanicsModelCohesive & self) -> decltype(auto) { \ return self.func_name(); \ }) void register_solid_mechanics_model_cohesive(py::module & mod) { py::class_(mod, "CohesiveElementInserter") .def("setLimit", &CohesiveElementInserter::setLimit) .def( "getCheckFacets", [](CohesiveElementInserter & self) -> decltype(auto) { return self.getCheckFacets(); }, py::return_value_policy::reference) .def( "getCheckFacets", [](CohesiveElementInserter & self, ElementType type, GhostType ghost_type) -> decltype(auto) { return self.getCheckFacets(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "getInsertionFacets", [](CohesiveElementInserter & self) -> decltype(auto) { return self.getInsertionFacetsByElement(); }, py::return_value_policy::reference) .def( "getInsertionFacets", [](CohesiveElementInserter & self, ElementType type, GhostType ghost_type) -> decltype(auto) { return self.getInsertionFacets(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def("addPhysicalSurface", &CohesiveElementInserter::addPhysicalSurface) .def("addPhysicalVolume", &CohesiveElementInserter::addPhysicalVolume); py::class_( mod, "SolidMechanicsModelCohesiveOptions") .def(py::init(), py::arg("analysis_method") = _explicit_lumped_mass, py::arg("is_extrinsic") = false); py::class_( mod, "SolidMechanicsModelCohesive") .def(py::init(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "solid_mechanics_model") .def( "initFull", [](SolidMechanicsModel & self, const AnalysisMethod & analysis_method, bool is_extrinsic) { self.initFull(_analysis_method = analysis_method, _is_extrinsic = is_extrinsic); }, py::arg("_analysis_method") = _explicit_lumped_mass, py::arg("_is_extrinsic") = false) .def("checkCohesiveStress", &SolidMechanicsModelCohesive::checkCohesiveStress) .def("getElementInserter", &SolidMechanicsModelCohesive::getElementInserter, py::return_value_policy::reference) .def("getStressOnFacets", &SolidMechanicsModelCohesive::getStressOnFacets, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def("getTangents", &SolidMechanicsModelCohesive::getTangents, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) + .def_function_nocopy(getLambda) .def("updateAutomaticInsertion", &SolidMechanicsModelCohesive::updateAutomaticInsertion); } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh index e93f23105..844b72b15 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,244 +1,244 @@ /** * Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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 "solid_mechanics_model_cohesive.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; + using MyFEEngineCohesiveType = + FEEngineTemplate; public: MaterialCohesive(SolidMechanicsModel & model, const ID & id = ""); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// 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) final{}; // add the facet to be handled by the material Idx addFacet(const Element & element); template inline decltype(auto) getArguments(ElementType element_type, GhostType ghost_type) { return zip( "normal"_n = make_view(this->normals(element_type, ghost_type)), "opening"_n = make_view(this->opening(element_type, ghost_type)), "traction"_n = make_view(this->tractions(element_type, ghost_type)), "contact_opening"_n = make_view(this->contact_opening(element_type, ghost_type)), "contact_traction"_n = make_view(this->contact_tractions(element_type, ghost_type)), "previous_opening"_n = make_view(this->opening.previous(element_type, ghost_type)), "previous_traction"_n = make_view(this->tractions.previous(element_type, ghost_type)), "delta_max"_n = this->delta_max(element_type, ghost_type), "damage"_n = this->damage(element_type, ghost_type)); } protected: virtual void computeTangentTraction(ElementType /*el_type*/, Array & /*tangent_matrix*/, 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(ElementType /*el_type*/, GhostType /*ghost_type*/ = _not_ghost) { AKANTU_TO_IMPLEMENT(); } /// parallelism functions [[nodiscard]] inline Int 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), Idx); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(FacetFilter, (*facet_filter), Idx); AKANTU_GET_MACRO_AUTO(FacetFilter, *facet_filter); AKANTU_GET_MACRO_AUTO_NOT_CONST(FacetFilter, *facet_filter); // AKANTU_GET_MACRO(ElementFilter, element_filter, const // ElementTypeMapArray &); /// compute reversible energy auto getReversibleEnergy() -> Real; /// compute dissipated energy auto getDissipatedEnergy() -> Real; /// compute contact energy auto getContactEnergy() -> Real; /// get energy auto getEnergy(const std::string & type) -> Real override; /// return the energy (identified by id) for the provided element // Real getEnergy(const std::string & energy_id, const Element & element // ) override { // return Material::getEnergy(energy_id, element); // } [[nodiscard]] virtual bool needLambda() const { return false; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// list of facets assigned to this material std::shared_ptr> 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; /// critical stress FacetRandomInternalField & sigma_c; /// pointer to the solid mechanics model for cohesive elements SolidMechanicsModelCohesive * model; /// critical displacement Real delta_c{0.}; /// array to temporarily store the normals CohesiveInternalField & normals; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "cohesive_internal_field_tmpl.hh" // NOLINT(unused-includes) #include "material_cohesive_inline_impl.hh" #endif /* AKANTU_MATERIAL_COHESIVE_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc index eb312510b..d5648e028 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc @@ -1,496 +1,70 @@ /** * Copyright (©) 2012-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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_damage.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" +//#include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ -template -MaterialCohesiveDamage::MaterialCohesiveDamage(SolidMechanicsModel & model, +MaterialCohesiveDamage::MaterialCohesiveDamage(SolidMechanicsModel & model, const ID & id) - : MaterialCohesive(model, id), - lambda(registerInternal("lambda", - spatial_dimension)), - err_openings(registerInternal( - "err_openings", spatial_dimension)), - czm_damage( - registerInternal("czm_damage", 1)) { + : MaterialCohesive(model, id) +{ AKANTU_DEBUG_IN(); this->registerParam("k", k, Real(0.), _pat_parsable | _pat_readable, - "Beta parameter"); + "Cohesive stiffness parameter"); this->registerParam("G_c", G_c, Real(0.), _pat_parsable | _pat_readable, "Mode I fracture energy"); + a = 2.*k*G_c/(sigma_c*sigma_c); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template void MaterialCohesiveDamage::initMaterial() { +void MaterialCohesiveDamage::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); AKANTU_DEBUG_OUT(); } - /* -------------------------------------------------------------------------- */ -template -void MaterialCohesiveDamage::assembleInternalForces(GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - auto & internal_force = const_cast &>(model->getInternalForce()); - - for (auto type : getElementFilter().elementTypes(spatial_dimension, - ghost_type, _ek_cohesive)) { - auto & elem_filter = getElementFilter(type, ghost_type); - auto 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 & err_opening = err_openings(type, ghost_type); - - auto size_of_shapes = shapes.getNbComponent(); - auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); - auto nb_quadrature_points = - fem_cohesive.getNbIntegrationPoints(type, ghost_type); - - /// compute @f$t_i N_a@f$ - - auto traction_cpy = std::make_shared>( - nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes); - - auto err_opening_cpy = std::make_shared>( - nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes); - - auto traction_it = traction.begin(spatial_dimension, 1); - auto err_opening_it = err_opening.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); - auto err_opening_cpy_it = - err_opening_cpy->begin(spatial_dimension, size_of_shapes); - - for (Int el = 0; el < nb_element; ++el) { - auto current_quad = elem_filter(el) * nb_quadrature_points; - - for (Int q = 0; q < nb_quadrature_points; ++q, ++traction_it, - ++err_opening_it, ++current_quad, ++traction_cpy_it, - ++err_opening_cpy_it) { - - auto && shapes_filtered = shapes_filtered_begin[current_quad]; - - *traction_cpy_it = (*traction_it) * shapes_filtered; - *err_opening_cpy_it = (*err_opening_it) * 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 = std::make_shared>( - 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); - - auto int_t_N = std::make_shared>( - nb_element, 2 * spatial_dimension * size_of_shapes, "int_t_N"); - - auto * int_t_N_val = int_t_N->data(); - auto * partial_int_t_N_val = partial_int_t_N->data(); - for (Int 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 (Int 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; - } - - /** - * compute @f$\int err \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t - * \mathbf{err}_q \overline w_q J_q@f$ - */ - auto int_err_N = std::make_shared>( - nb_element, spatial_dimension * size_of_shapes, "int_err_N"); - - fem_cohesive.integrate(*err_opening_cpy, *int_err_N, - spatial_dimension * size_of_shapes, type, ghost_type, - elem_filter); - - /// assemble - model->getDOFManager().assembleElementalArrayLocalArray( - "displacement", *int_t_N, internal_force, type, ghost_type, 1, - elem_filter); - - // auto lambda_connectivity = lambda_connectivities(type, ghost_type); - auto underlying_type = Mesh::getFacetType(type); - model->getDOFManager().assembleElementalArrayToResidual( - "lambda", *int_err_N, underlying_type, ghost_type, 1., elem_filter); - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void MaterialCohesiveDamage::assembleStiffnessMatrix( - GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - // auto & lambda_connectivities = - // model->getMesh().getElementalData("lambda_connectivities"); - - for (auto type : getElementFilter().elementTypes(spatial_dimension, - ghost_type, _ek_cohesive)) { - auto nb_quadrature_points = - fem_cohesive.getNbIntegrationPoints(type, ghost_type); - auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); - - const auto & elem_filter = getElementFilter(type, ghost_type); - auto nb_element = elem_filter.size(); - - if (nb_element == 0U) { - continue; - } - - const auto & shapes = fem_cohesive.getShapes(type, ghost_type); - auto size_of_shapes = shapes.getNbComponent(); - - auto shapes_filtered = std::make_shared>( - nb_element * nb_quadrature_points, size_of_shapes, "filtered shapes"); - - for (auto && data : - zip(filter(elem_filter, - make_view(shapes, size_of_shapes, nb_quadrature_points)), - make_view(*shapes_filtered, size_of_shapes, - nb_quadrature_points))) { - std::get<1>(data) = std::get<0>(data); - } - - Matrix A(spatial_dimension * size_of_shapes, - spatial_dimension * nb_nodes_per_element); - A.zero(); - for (Int 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$ - /// TODO : optimisation not to reassemble uu term, which does not change - /// during computation - auto tangent_stiffness_matrix_uu = std::make_unique>( - nb_element * nb_quadrature_points, - spatial_dimension * spatial_dimension, "tangent_stiffness_matrix_uu"); - - auto tangent_stiffness_matrix_ll = std::make_unique>( - nb_element * nb_quadrature_points, - spatial_dimension * spatial_dimension, "tangent_stiffness_matrix_ll"); - - computeNormal(model->getCurrentPosition(), normals(type, ghost_type), type, - ghost_type); - - /// compute openings @f$\mathbf{\delta}@f$ - // computeOpening(model->getDisplacement(), opening(type, ghost_type), type, - // ghost_type); - - tangent_stiffness_matrix_uu->zero(); - tangent_stiffness_matrix_ll->zero(); - - computeTangentTraction(type, *tangent_stiffness_matrix_uu, - *tangent_stiffness_matrix_ll, ghost_type); - - UInt size_at_nt_duu_n_a = spatial_dimension * nb_nodes_per_element * - spatial_dimension * nb_nodes_per_element; - auto at_nt_duu_n_a = - std::make_unique>(nb_element * nb_quadrature_points, - size_at_nt_duu_n_a, "A^t*N^t*Duu*N*A"); - - UInt size_nt_dll_n = - spatial_dimension * size_of_shapes * spatial_dimension * size_of_shapes; - auto nt_dll_n = std::make_unique>( - nb_element * nb_quadrature_points, size_nt_dll_n, "N^t*Dll*N"); - - UInt size_at_nt_dul_n = spatial_dimension * nb_nodes_per_element * - spatial_dimension * size_of_shapes; - auto at_nt_dul_n = std::make_unique>( - nb_element * nb_quadrature_points, size_at_nt_dul_n, "A^t*N^t*Dul*N"); - - Matrix N(spatial_dimension, spatial_dimension * size_of_shapes); - - for (auto && [At_Nt_Duu_N_A, Duu, Nt_Dll_N, Dll, At_Nt_Dul_N, shapes] : - zip(make_view(*at_nt_duu_n_a, spatial_dimension * nb_nodes_per_element, - spatial_dimension * nb_nodes_per_element), - make_view(*tangent_stiffness_matrix_uu, spatial_dimension, - spatial_dimension), - make_view(*nt_dll_n, spatial_dimension * size_of_shapes, - spatial_dimension * size_of_shapes), - make_view(*tangent_stiffness_matrix_ll, spatial_dimension, - spatial_dimension), - make_view(*at_nt_dul_n, spatial_dimension * nb_nodes_per_element, - spatial_dimension * size_of_shapes), - make_view(*shapes_filtered, size_of_shapes))) { - 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 (Int i = 0; i < spatial_dimension; ++i) { - for (Int n = 0; n < size_of_shapes; ++n) { - N(i, i + spatial_dimension * n) = shapes(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$ - **/ - auto && NA = N * A; - At_Nt_Duu_N_A = (Duu * NA).transpose() * NA; - Nt_Dll_N = (Dll * N).transpose() * N; - At_Nt_Dul_N = NA.transpose() * N; - } - - auto Kuu_e = - std::make_unique>(nb_element, size_at_nt_duu_n_a, "Kuu_e"); - - fem_cohesive.integrate(*at_nt_duu_n_a, *Kuu_e, size_at_nt_duu_n_a, type, - ghost_type, elem_filter); - - auto Kll_e = - std::make_unique>(nb_element, size_nt_dll_n, "Kll_e"); - - fem_cohesive.integrate(*nt_dll_n, *Kll_e, size_nt_dll_n, type, ghost_type, - elem_filter); - - auto Kul_e = - std::make_unique>(nb_element, size_at_nt_dul_n, "Kul_e"); - - fem_cohesive.integrate(*at_nt_dul_n, *Kul_e, size_at_nt_dul_n, type, - ghost_type, elem_filter); - - model->getDOFManager().assembleElementalMatricesToMatrix( - "K", "displacement", *Kuu_e, type, ghost_type, _symmetric, elem_filter); - - auto underlying_type = Mesh::getFacetType(type); - model->getDOFManager().assembleElementalMatricesToMatrix( - "K", "lambda", *Kll_e, underlying_type, ghost_type, _symmetric, - elem_filter); - - auto connectivity = model->getMesh().getConnectivity(type, ghost_type); - auto conn = make_view(connectivity, connectivity.getNbComponent()).begin(); - - auto lambda_connectivity = - model->getLambdaMesh().getConnectivity(underlying_type, ghost_type); - auto lambda_conn = - make_view(lambda_connectivity, lambda_connectivity.getNbComponent()) - .begin(); - - /// Assemble Kll_e - // TermsToAssemble term_ll("lambda", "lambda"); - // auto el_mat_it_ll = Kll_e->begin(spatial_dimension * size_of_shapes, - // spatial_dimension * size_of_shapes); - - // auto compute_ll = [&](const auto & el) { - // auto kll_e = *el_mat_it_ll; - // auto && lda_conn_el = lambda_conn[el]; - // auto N = lda_conn_el.rows(); - // for (Int m = 0; m < N; ++m) { - // auto ldai = lda_conn_el(m); - // for (Int n = m; n < N; ++n) { - // auto ldaj = lda_conn_el(n); - // for (Int k = 0; k < spatial_dimension; ++k) { - // for (Int l = 0; l < spatial_dimension; ++l) { - // auto && term_ll_ij = term_ll(ldai * spatial_dimension + k, - // ldaj * spatial_dimension + l); - // term_ll_ij = - // kll_e(m * spatial_dimension + k, n * spatial_dimension + - // l); - // } - // } - // } - // } - // ++el_mat_it_ll; - // }; - // for_each_element(nb_element, elem_filter, compute_ll); - - // model->getDOFManager().assemblePreassembledMatrix("K", term_ll); - // model->getDOFManager().getMatrix("K").saveMatrix("Kuu_terms.mtx"); - - /// Assemble Klu_e - TermsToAssemble term_ul("displacement", "lambda"); - auto el_mat_it_ul = Kul_e->begin(spatial_dimension * nb_nodes_per_element, - spatial_dimension * size_of_shapes); - - auto compute_ul = [&](const auto & el) { - auto kul_e = *el_mat_it_ul; - auto && u_conn_el = conn[el]; - auto && lda_conn_el = lambda_conn[el]; - auto M = u_conn_el.rows(); - auto N = lda_conn_el.rows(); - for (Int m = 0; m < M; ++m) { - for (Int n = 0; n < N; ++n) { - auto u = u_conn_el(m); - auto lda = lda_conn_el(n); - for (Int k = 0; k < spatial_dimension; ++k) { - for (Int l = 0; l < spatial_dimension; ++l) { - auto && term_ul_ij = term_ul(u * spatial_dimension + k, - lda * spatial_dimension + l); - term_ul_ij = - kul_e(m * spatial_dimension + k, n * spatial_dimension + l); - } - } - } - } - ++el_mat_it_ul; - }; - for_each_element(nb_element, elem_filter, compute_ul); - - model->getDOFManager().assemblePreassembledMatrix("K", term_ul); - } - - AKANTU_DEBUG_OUT(); +Real MaterialCohesiveDamage::stiffness(Real d) +{ + return k*(1./d-1.); } - -/* -------------------------------------------------------------------------- */ -template -void MaterialCohesiveDamage::computeLambdaOnQuad(ElementType type, - GhostType ghost_type) { - auto & fem_lambda = this->model->getFEEngine("LambdaFEEngine"); - const auto & lambda = this->model->getLambda(); - auto & lambda_on_quad = this->lambda(type, ghost_type); - -// std::cout << "lambda in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl; -// lambda.printself(std::cout << std::endl); -// ArrayPrintHelper::print_content(lambda,std::cout,0); - - - auto underlying_type = Mesh::getFacetType(type); - fem_lambda.interpolateOnIntegrationPoints( - lambda, lambda_on_quad, dim, underlying_type, ghost_type, - this->getElementFilter(type, ghost_type)); - -// std::cout << "lambda_on_quad in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl; -// lambda_on_quad.printself(std::cout << std::endl); -// ArrayPrintHelper::print_content(lambda_on_quad,std::cout,0); - -// std::cout << "shapes_lambda in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl; -// auto shapes_lambda = fem_lambda.getShapes(underlying_type,ghost_type,0); -// ArrayPrintHelper::print_content(shapes_lambda,std::cout,0); - -// std::cout << "shapes_cohesive in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl; -// auto & fem_cohesive = -// this->model->getFEEngineClass("CohesiveFEEngine"); -// auto shapes_cohesive = fem_cohesive.getShapes(type,ghost_type,0); -// ArrayPrintHelper::print_content(shapes_cohesive,std::cout,0); +Real MaterialCohesiveDamage::augmented_stiffness(Real d) +{ + return k*1./d; } - -/* -------------------------------------------------------------------------- */ -template -void MaterialCohesiveDamage::computeTraction(ElementType el_type, - GhostType ghost_type) { - - for (const auto & type : getElementFilter().elementTypes( - spatial_dimension, ghost_type, _ek_cohesive)) { - auto & elem_filter = getElementFilter(type, ghost_type); - auto nb_element = elem_filter.size(); - if (nb_element == 0) { - continue; - } - - /// compute normals @f$\mathbf{n}@f$ - computeNormal(model->getCurrentPosition(), normals(type, ghost_type), type, - ghost_type); - - /// compute openings @f$\mathbf{\delta}@f$ - computeOpening(model->getDisplacement(), opening(type, ghost_type), type, - ghost_type); - - computeLambdaOnQuad(el_type, ghost_type); - -// auto & traction = tractions(type, ghost_type); -// std::cout << "traction in : MaterialCohesiveDamage::computeTraction " << std::endl; -// ArrayPrintHelper::print_content(traction,std::cout,0); - -// auto & lambda_ = lambda(type, ghost_type); -// std::cout << "lambda in : MaterialCohesiveDamage::computeTraction " << std::endl; -// ArrayPrintHelper::print_content(lambda_,std::cout,0); - - for (auto && args : getArguments(el_type, ghost_type)) { - this->computeTractionOnQuad(args); - } - } +Real MaterialCohesiveDamage::augmented_compliance(Real d) +{ + return d/k; } -/* -------------------------------------------------------------------------- */ -template -void MaterialCohesiveDamage::computeTangentTraction( - ElementType el_type, Array & tangent_matrix_uu, - Array & tangent_matrix_ll, GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - for (auto && [args, tangent_uu, tangent_ll] : - zip(getArguments(el_type, ghost_type), - make_view(tangent_matrix_uu), - make_view(tangent_matrix_ll))) { - computeTangentTractionOnQuad(tangent_uu, tangent_ll, args); - } - - AKANTU_DEBUG_OUT(); } - -/* -------------------------------------------------------------------------- */ -template class MaterialCohesiveDamage<1>; -template class MaterialCohesiveDamage<2>; -template class MaterialCohesiveDamage<3>; -const bool material_is_alocated_cohesive_damage [[maybe_unused]] = - instantiateMaterial("cohesive_damage"); - -} // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh index 965cd4e85..538ad747d 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh @@ -1,126 +1,77 @@ /** * Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_ #define AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_ namespace akantu { /** - * Cohesive material linear damage for extrinsic case + * Cohesive material linear damage * * parameters in the material files : * - sigma_c : critical stress sigma_c (default: 0) - * - beta : weighting parameter for sliding and normal opening (default: - * 0) - * - G_cI : fracture energy for mode I (default: 0) - * - G_cII : fracture energy for mode II (default: 0) - * - penalty : stiffness in compression to prevent penetration + * - G_c : fracture energy for mode I (default: 0) + * - k : cohesive stiffness */ -template class MaterialCohesiveDamage : public MaterialCohesive { +class MaterialCohesiveDamage : public MaterialCohesive { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialCohesiveDamage(SolidMechanicsModel & model, const ID & id = ""); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material parameters void initMaterial() override; - /// assemble stiffness - void assembleStiffnessMatrix(GhostType ghost_type) override; + /// these functions should be put in a separate class + Real stiffness(Real d); + Real augmented_stiffness(Real d); + Real augmented_compliance(Real d); - /// assemble residual - void assembleInternalForces(GhostType ghost_type = _not_ghost) override; - -protected: - void computeLambdaOnQuad(ElementType type, GhostType ghost_type); - - inline decltype(auto) getArguments(ElementType element_type, - GhostType ghost_type) { - using namespace tuple; - return zip_append( - MaterialCohesive::getArguments(element_type, ghost_type), - "lambda"_n = make_view(this->lambda(element_type, ghost_type)), - "err_opening"_n = - make_view(this->err_openings(element_type, ghost_type))); - } - - /// constitutive law - void computeTraction(ElementType el_type, - GhostType ghost_type = _not_ghost) override; - - /// compute tangent stiffness matrix - /// WARNING : override was removed, not sure it should have - void computeTangentTraction(ElementType el_type, - Array & tangent_matrix_uu, - Array & tangent_matrix_ll, - GhostType ghost_type); - - /// compute the traction for a given quadrature point - template inline void computeTractionOnQuad(Args && args); - - template - inline void - computeTangentTractionOnQuad(Eigen::MatrixBase & tangent_uu, - Eigen::MatrixBase & tangent_ll, - Args && args); - - [[nodiscard]] bool needLambda() const override { return true; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// cohesive stiffness Real k; /// mode I fracture energy Real G_c; - /// augmented lagrange multiplier - CohesiveInternalField & lambda; - - /// target opening - CohesiveInternalField & err_openings; - - /// cohesive damage - CohesiveInternalField & czm_damage; - - Vector normal_opening; - Real normal_opening_norm{0.}; + /// parameter for softening function + Real a; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ } // namespace akantu -#include "material_cohesive_damage_inline_impl.hh" - #endif /* AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.cc new file mode 100644 index 000000000..abdf02aba --- /dev/null +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.cc @@ -0,0 +1,169 @@ +/** + * Copyright (©) 2012-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * This file is part of Akantu + * + * 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_damage_extrinsic.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" +/* -------------------------------------------------------------------------- */ +#include +#include +/* -------------------------------------------------------------------------- */ + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +template +MaterialCohesiveDamageExtrinsic::MaterialCohesiveDamageExtrinsic(SolidMechanicsModel & model, + const ID & id) + : MaterialCohesiveDamage(model, id), + lambda(registerInternal("lambda", + spatial_dimension)), + err_openings(registerInternal( + "err_openings", spatial_dimension)), + czm_damage( + registerInternal("czm_damage", 1)) +{ + if(not this->model->getIsExtrinsic()) + { + AKANTU_EXCEPTION( + "MaterialCohesiveDamageExtrinsic can be used only with extrinsic cohesive elements"); + } +} + +/* -------------------------------------------------------------------------- */ +template // NOLINTNEXTLINE(readability-function-cognitive-complexity) +void MaterialCohesiveDamageExtrinsic::checkInsertion(bool check_only) { + AKANTU_DEBUG_IN(); + +// const auto & mesh_facets = model->getMeshFacets(); +// auto & inserter = model->getElementInserter(); + +// for (const auto & type_facet : mesh_facets.elementTypes(dim - 1)) { +// auto type_cohesive = FEEngine::getCohesiveElementType(type_facet); +// auto & f_insertion = inserter.getInsertionFacets(type_facet); +// auto & damage = czm_damage(type_cohesive); + +// const auto & facets_check = inserter.getCheckFacets(type_facet); +// const auto & facet_filter_array = getFacetFilter(type_facet); + +// auto nb_quad_facet = +// model->getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet); + +//#ifndef AKANTU_NDEBUG +// auto nb_quad_cohesive = model->getFEEngine("CohesiveFEEngine") +// .getNbIntegrationPoints(type_cohesive); + +// AKANTU_DEBUG_ASSERT(nb_quad_cohesive == nb_quad_facet, +// "The cohesive element and the corresponding facet do " +// "not have the same numbers of integration points"); +//#endif + +// // loop over each facet belonging to this material +// for (auto && [facet] : facet_filter_array) { +// // skip facets where check shouldn't be realized +// if (!facets_check(facet)) { +// continue; +// } + +// Vector d_check(nb_quad_facet); +// // compute the effective norm on each quadrature point of the facet +// for (Int q = 0; q < nb_quad_facet; ++q) { +// auto current_quad = facet * nb_quad_facet + q; +// auto && normal = normal_begin[current_quad]; +// auto && tangent = tangent_begin[current_quad]; +// auto && facet_stress = facet_stress_begin[current_quad]; + +// // compute average stress on the current quadrature point +// auto && d_1 = facet_stress(0); +// auto && d_2 = facet_stress(1); + +// auto && stress_avg = (stress_1 + stress_2) / 2.; + +// // compute normal and effective stress +// stress_check(q) = computeEffectiveNorm(stress_avg, normal, tangent, +// normal_traction(q)); +// } + +// // verify if the effective stress overcomes the threshold +// auto final_stress = stress_check.mean(); +// if (max_quad_stress_insertion) { +// final_stress = +// *std::max_element(stress_check.begin(), stress_check.end()); +// } + +// if (final_stress > sigma_limit) { +// f_insertion(facet) = true; + +// if (check_only) { +// continue; +// } +// } +// } +// } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialCohesiveDamageExtrinsic::computeLambdaOnQuad(ElementType type, + GhostType ghost_type) { + auto & fem_lambda = model->getFEEngine("LambdaFEEngine"); + const auto & lambda = model->getLambda(); + auto & lambda_on_quad = this->lambda(type, ghost_type); + + auto underlying_type = Mesh::getFacetType(type); + fem_lambda.interpolateOnIntegrationPoints( + lambda, lambda_on_quad, dim, underlying_type, ghost_type, + this->getElementFilter(type, ghost_type)); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialCohesiveDamageExtrinsic::computeTraction(ElementType el_type, + GhostType ghost_type) { + + for (const auto & type : getElementFilter().elementTypes( + spatial_dimension, ghost_type, _ek_cohesive)) { + auto & elem_filter = getElementFilter(type, ghost_type); + auto nb_element = elem_filter.size(); + if (nb_element == 0) { + continue; + } + + for (auto && args : getArguments(el_type, ghost_type)) { + computeTractionOnQuad(args); + } + } +} + +/* -------------------------------------------------------------------------- */ +template class MaterialCohesiveDamageExtrinsic<1>; +template class MaterialCohesiveDamageExtrinsic<2>; +template class MaterialCohesiveDamageExtrinsic<3>; +const bool material_is_alocated_cohesive_damage [[maybe_unused]] = + instantiateMaterial("cohesive_damage_extrinsic"); + +} // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.hh similarity index 62% copy from src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh copy to src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.hh index 965cd4e85..17d69e9da 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.hh @@ -1,126 +1,98 @@ /** * Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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 "material_cohesive_damage.hh" /* -------------------------------------------------------------------------- */ -#ifndef AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_ -#define AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_ +#ifndef AKANTU_MATERIAL_COHESIVE_DAMAGE_EXTRINSIC_HH_ +#define AKANTU_MATERIAL_COHESIVE_DAMAGE_EXTRINSIC_HH_ namespace akantu { /** - * Cohesive material linear damage for extrinsic case + * Cohesive material linear damage for extrinsinc case * * parameters in the material files : * - sigma_c : critical stress sigma_c (default: 0) - * - beta : weighting parameter for sliding and normal opening (default: - * 0) - * - G_cI : fracture energy for mode I (default: 0) - * - G_cII : fracture energy for mode II (default: 0) - * - penalty : stiffness in compression to prevent penetration + * - G_c : fracture energy for mode I (default: 0) + * - k : cohesive stiffness */ -template class MaterialCohesiveDamage : public MaterialCohesive { +template class MaterialCohesiveDamageExtrinsic : public MaterialCohesiveDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - MaterialCohesiveDamage(SolidMechanicsModel & model, const ID & id = ""); + MaterialCohesiveDamageExtrinsic(SolidMechanicsModel & model, const ID & id = ""); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ -public: - /// initialize the material parameters - void initMaterial() override; - - /// assemble stiffness - void assembleStiffnessMatrix(GhostType ghost_type) override; - /// assemble residual - void assembleInternalForces(GhostType ghost_type = _not_ghost) override; + // /// check stress for cohesive elements' insertion + void checkInsertion(bool check_only = false) override; protected: void computeLambdaOnQuad(ElementType type, GhostType ghost_type); inline decltype(auto) getArguments(ElementType element_type, GhostType ghost_type) { using namespace tuple; return zip_append( MaterialCohesive::getArguments(element_type, ghost_type), - "lambda"_n = make_view(this->lambda(element_type, ghost_type)), - "err_opening"_n = - make_view(this->err_openings(element_type, ghost_type))); + "czm_damage"_n = this->czm_damage(element_type, ghost_type)); } /// constitutive law void computeTraction(ElementType el_type, GhostType ghost_type = _not_ghost) override; - /// compute tangent stiffness matrix - /// WARNING : override was removed, not sure it should have - void computeTangentTraction(ElementType el_type, - Array & tangent_matrix_uu, - Array & tangent_matrix_ll, - GhostType ghost_type); - /// compute the traction for a given quadrature point template inline void computeTractionOnQuad(Args && args); - template - inline void - computeTangentTractionOnQuad(Eigen::MatrixBase & tangent_uu, - Eigen::MatrixBase & tangent_ll, - Args && args); + [[nodiscard]] bool needLambda() const override { return false; } - [[nodiscard]] bool needLambda() const override { return true; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: - /// cohesive stiffness - Real k; - - /// mode I fracture energy - Real G_c; /// augmented lagrange multiplier CohesiveInternalField & lambda; /// target opening CohesiveInternalField & err_openings; /// cohesive damage CohesiveInternalField & czm_damage; Vector normal_opening; Real normal_opening_norm{0.}; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ } // namespace akantu -#include "material_cohesive_damage_inline_impl.hh" +#include "material_cohesive_damage_extrinsic_inline_impl.hh" -#endif /* AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_ */ +#endif /* AKANTU_MATERIAL_COHESIVE_DAMAGE_EXTRINSIC_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic_inline_impl.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic_inline_impl.hh new file mode 100644 index 000000000..1914d8552 --- /dev/null +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic_inline_impl.hh @@ -0,0 +1,47 @@ +/** + * Copyright (©) 2015-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * This file is part of Akantu + * + * 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_damage_extrinsic.hh" + +/* -------------------------------------------------------------------------- */ +#ifndef AKANTU_MATERIAL_COHESIVE_DAMAGE_EXTRINSIC_INLINE_IMPL_HH_ +#define AKANTU_MATERIAL_COHESIVE_DAMAGE_EXTRINSIC_INLINE_IMPL_HH_ + +/* -------------------------------------------------------------------------- */ + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +template +template +inline void MaterialCohesiveDamageExtrinsic::computeTractionOnQuad(Args && args) { + auto && d = args["czm_damage"_n]; + auto && opening = args["opening"_n]; + auto && traction = args["traction"_n]; + + Real A = this->augmented_stiffness(d); + traction = k*A*opening; +} + +} // namespace akantu + +/* -------------------------------------------------------------------------- */ +#endif // AKANTU_MATERIAL_COHESIVE_DAMAGE_EXTRINSIC_INLINE_IMPL_HH_ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.cc similarity index 89% copy from src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc copy to src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.cc index eb312510b..8f823be64 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.cc @@ -1,496 +1,484 @@ /** * Copyright (©) 2012-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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_damage.hh" +#include "material_cohesive_damage_intrinsic.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" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template -MaterialCohesiveDamage::MaterialCohesiveDamage(SolidMechanicsModel & model, +MaterialCohesiveDamageIntrinsic::MaterialCohesiveDamageIntrinsic(SolidMechanicsModel & model, const ID & id) - : MaterialCohesive(model, id), + : MaterialCohesiveDamage(model, id), lambda(registerInternal("lambda", spatial_dimension)), err_openings(registerInternal( "err_openings", spatial_dimension)), czm_damage( registerInternal("czm_damage", 1)) { - AKANTU_DEBUG_IN(); - - this->registerParam("k", k, Real(0.), _pat_parsable | _pat_readable, - "Beta parameter"); - - this->registerParam("G_c", G_c, Real(0.), _pat_parsable | _pat_readable, - "Mode I fracture energy"); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template void MaterialCohesiveDamage::initMaterial() { - AKANTU_DEBUG_IN(); - - MaterialCohesive::initMaterial(); - - AKANTU_DEBUG_OUT(); + if(this->model->getIsExtrinsic()) + { + AKANTU_EXCEPTION( + "MaterialCohesiveDamageIntrinsic can be used only with intrinsic cohesive elements"); + } } /* -------------------------------------------------------------------------- */ template -void MaterialCohesiveDamage::assembleInternalForces(GhostType ghost_type) { +void MaterialCohesiveDamageIntrinsic::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & internal_force = const_cast &>(model->getInternalForce()); for (auto type : getElementFilter().elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { auto & elem_filter = getElementFilter(type, ghost_type); auto 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 & err_opening = err_openings(type, ghost_type); auto size_of_shapes = shapes.getNbComponent(); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem_cohesive.getNbIntegrationPoints(type, ghost_type); /// compute @f$t_i N_a@f$ auto traction_cpy = std::make_shared>( nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes); auto err_opening_cpy = std::make_shared>( nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes); auto traction_it = traction.begin(spatial_dimension, 1); auto err_opening_it = err_opening.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); auto err_opening_cpy_it = err_opening_cpy->begin(spatial_dimension, size_of_shapes); for (Int el = 0; el < nb_element; ++el) { auto current_quad = elem_filter(el) * nb_quadrature_points; for (Int q = 0; q < nb_quadrature_points; ++q, ++traction_it, ++err_opening_it, ++current_quad, ++traction_cpy_it, ++err_opening_cpy_it) { auto && shapes_filtered = shapes_filtered_begin[current_quad]; *traction_cpy_it = (*traction_it) * shapes_filtered; *err_opening_cpy_it = (*err_opening_it) * 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 = std::make_shared>( 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); auto int_t_N = std::make_shared>( nb_element, 2 * spatial_dimension * size_of_shapes, "int_t_N"); auto * int_t_N_val = int_t_N->data(); auto * partial_int_t_N_val = partial_int_t_N->data(); for (Int 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 (Int 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; } /** * compute @f$\int err \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t * \mathbf{err}_q \overline w_q J_q@f$ */ auto int_err_N = std::make_shared>( nb_element, spatial_dimension * size_of_shapes, "int_err_N"); fem_cohesive.integrate(*err_opening_cpy, *int_err_N, spatial_dimension * size_of_shapes, type, ghost_type, elem_filter); /// assemble model->getDOFManager().assembleElementalArrayLocalArray( "displacement", *int_t_N, internal_force, type, ghost_type, 1, elem_filter); // auto lambda_connectivity = lambda_connectivities(type, ghost_type); auto underlying_type = Mesh::getFacetType(type); model->getDOFManager().assembleElementalArrayToResidual( "lambda", *int_err_N, underlying_type, ghost_type, 1., elem_filter); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template -void MaterialCohesiveDamage::assembleStiffnessMatrix( +void MaterialCohesiveDamageIntrinsic::assembleStiffnessMatrix( GhostType ghost_type) { AKANTU_DEBUG_IN(); // auto & lambda_connectivities = // model->getMesh().getElementalData("lambda_connectivities"); for (auto type : getElementFilter().elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { auto nb_quadrature_points = fem_cohesive.getNbIntegrationPoints(type, ghost_type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const auto & elem_filter = getElementFilter(type, ghost_type); auto nb_element = elem_filter.size(); if (nb_element == 0U) { continue; } const auto & shapes = fem_cohesive.getShapes(type, ghost_type); auto size_of_shapes = shapes.getNbComponent(); auto shapes_filtered = std::make_shared>( nb_element * nb_quadrature_points, size_of_shapes, "filtered shapes"); for (auto && data : zip(filter(elem_filter, make_view(shapes, size_of_shapes, nb_quadrature_points)), make_view(*shapes_filtered, size_of_shapes, nb_quadrature_points))) { std::get<1>(data) = std::get<0>(data); } Matrix A(spatial_dimension * size_of_shapes, spatial_dimension * nb_nodes_per_element); A.zero(); for (Int 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$ /// TODO : optimisation not to reassemble uu term, which does not change /// during computation auto tangent_stiffness_matrix_uu = std::make_unique>( nb_element * nb_quadrature_points, spatial_dimension * spatial_dimension, "tangent_stiffness_matrix_uu"); auto tangent_stiffness_matrix_ll = std::make_unique>( nb_element * nb_quadrature_points, spatial_dimension * spatial_dimension, "tangent_stiffness_matrix_ll"); computeNormal(model->getCurrentPosition(), normals(type, ghost_type), type, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ // computeOpening(model->getDisplacement(), opening(type, ghost_type), type, // ghost_type); tangent_stiffness_matrix_uu->zero(); tangent_stiffness_matrix_ll->zero(); computeTangentTraction(type, *tangent_stiffness_matrix_uu, *tangent_stiffness_matrix_ll, ghost_type); UInt size_at_nt_duu_n_a = spatial_dimension * nb_nodes_per_element * spatial_dimension * nb_nodes_per_element; auto at_nt_duu_n_a = std::make_unique>(nb_element * nb_quadrature_points, size_at_nt_duu_n_a, "A^t*N^t*Duu*N*A"); UInt size_nt_dll_n = spatial_dimension * size_of_shapes * spatial_dimension * size_of_shapes; auto nt_dll_n = std::make_unique>( nb_element * nb_quadrature_points, size_nt_dll_n, "N^t*Dll*N"); UInt size_at_nt_dul_n = spatial_dimension * nb_nodes_per_element * spatial_dimension * size_of_shapes; auto at_nt_dul_n = std::make_unique>( nb_element * nb_quadrature_points, size_at_nt_dul_n, "A^t*N^t*Dul*N"); Matrix N(spatial_dimension, spatial_dimension * size_of_shapes); for (auto && [At_Nt_Duu_N_A, Duu, Nt_Dll_N, Dll, At_Nt_Dul_N, shapes] : zip(make_view(*at_nt_duu_n_a, spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element), make_view(*tangent_stiffness_matrix_uu, spatial_dimension, spatial_dimension), make_view(*nt_dll_n, spatial_dimension * size_of_shapes, spatial_dimension * size_of_shapes), make_view(*tangent_stiffness_matrix_ll, spatial_dimension, spatial_dimension), make_view(*at_nt_dul_n, spatial_dimension * nb_nodes_per_element, spatial_dimension * size_of_shapes), make_view(*shapes_filtered, size_of_shapes))) { 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 (Int i = 0; i < spatial_dimension; ++i) { for (Int n = 0; n < size_of_shapes; ++n) { N(i, i + spatial_dimension * n) = shapes(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$ **/ auto && NA = N * A; At_Nt_Duu_N_A = (Duu * NA).transpose() * NA; Nt_Dll_N = (Dll * N).transpose() * N; At_Nt_Dul_N = NA.transpose() * N; } auto Kuu_e = std::make_unique>(nb_element, size_at_nt_duu_n_a, "Kuu_e"); fem_cohesive.integrate(*at_nt_duu_n_a, *Kuu_e, size_at_nt_duu_n_a, type, ghost_type, elem_filter); auto Kll_e = std::make_unique>(nb_element, size_nt_dll_n, "Kll_e"); fem_cohesive.integrate(*nt_dll_n, *Kll_e, size_nt_dll_n, type, ghost_type, elem_filter); auto Kul_e = std::make_unique>(nb_element, size_at_nt_dul_n, "Kul_e"); fem_cohesive.integrate(*at_nt_dul_n, *Kul_e, size_at_nt_dul_n, type, ghost_type, elem_filter); model->getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *Kuu_e, type, ghost_type, _symmetric, elem_filter); auto underlying_type = Mesh::getFacetType(type); model->getDOFManager().assembleElementalMatricesToMatrix( "K", "lambda", *Kll_e, underlying_type, ghost_type, _symmetric, elem_filter); auto connectivity = model->getMesh().getConnectivity(type, ghost_type); auto conn = make_view(connectivity, connectivity.getNbComponent()).begin(); auto lambda_connectivity = model->getLambdaMesh().getConnectivity(underlying_type, ghost_type); auto lambda_conn = make_view(lambda_connectivity, lambda_connectivity.getNbComponent()) .begin(); /// Assemble Kll_e // TermsToAssemble term_ll("lambda", "lambda"); // auto el_mat_it_ll = Kll_e->begin(spatial_dimension * size_of_shapes, // spatial_dimension * size_of_shapes); // auto compute_ll = [&](const auto & el) { // auto kll_e = *el_mat_it_ll; // auto && lda_conn_el = lambda_conn[el]; // auto N = lda_conn_el.rows(); // for (Int m = 0; m < N; ++m) { // auto ldai = lda_conn_el(m); // for (Int n = m; n < N; ++n) { // auto ldaj = lda_conn_el(n); // for (Int k = 0; k < spatial_dimension; ++k) { // for (Int l = 0; l < spatial_dimension; ++l) { // auto && term_ll_ij = term_ll(ldai * spatial_dimension + k, // ldaj * spatial_dimension + l); // term_ll_ij = // kll_e(m * spatial_dimension + k, n * spatial_dimension + // l); // } // } // } // } // ++el_mat_it_ll; // }; // for_each_element(nb_element, elem_filter, compute_ll); // model->getDOFManager().assemblePreassembledMatrix("K", term_ll); // model->getDOFManager().getMatrix("K").saveMatrix("Kuu_terms.mtx"); /// Assemble Klu_e TermsToAssemble term_ul("displacement", "lambda"); auto el_mat_it_ul = Kul_e->begin(spatial_dimension * nb_nodes_per_element, spatial_dimension * size_of_shapes); auto compute_ul = [&](const auto & el) { auto kul_e = *el_mat_it_ul; auto && u_conn_el = conn[el]; auto && lda_conn_el = lambda_conn[el]; auto M = u_conn_el.rows(); auto N = lda_conn_el.rows(); for (Int m = 0; m < M; ++m) { for (Int n = 0; n < N; ++n) { auto u = u_conn_el(m); auto lda = lda_conn_el(n); for (Int k = 0; k < spatial_dimension; ++k) { for (Int l = 0; l < spatial_dimension; ++l) { auto && term_ul_ij = term_ul(u * spatial_dimension + k, lda * spatial_dimension + l); term_ul_ij = kul_e(m * spatial_dimension + k, n * spatial_dimension + l); } } } } ++el_mat_it_ul; }; for_each_element(nb_element, elem_filter, compute_ul); model->getDOFManager().assemblePreassembledMatrix("K", term_ul); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template -void MaterialCohesiveDamage::computeLambdaOnQuad(ElementType type, +void MaterialCohesiveDamageIntrinsic::computeLambdaOnQuad(ElementType type, GhostType ghost_type) { - auto & fem_lambda = this->model->getFEEngine("LambdaFEEngine"); - const auto & lambda = this->model->getLambda(); + auto & fem_lambda = model->getFEEngine("LambdaFEEngine"); + const auto & lambda = model->getLambda(); auto & lambda_on_quad = this->lambda(type, ghost_type); -// std::cout << "lambda in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl; +// std::cout << "lambda in : MaterialCohesiveDamageIntrinsic::computeLambdaOnQuad " << std::endl; // lambda.printself(std::cout << std::endl); // ArrayPrintHelper::print_content(lambda,std::cout,0); auto underlying_type = Mesh::getFacetType(type); fem_lambda.interpolateOnIntegrationPoints( lambda, lambda_on_quad, dim, underlying_type, ghost_type, this->getElementFilter(type, ghost_type)); -// std::cout << "lambda_on_quad in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl; +// std::cout << "lambda_on_quad in : MaterialCohesiveDamageIntrinsic::computeLambdaOnQuad " << std::endl; // lambda_on_quad.printself(std::cout << std::endl); // ArrayPrintHelper::print_content(lambda_on_quad,std::cout,0); -// std::cout << "shapes_lambda in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl; +// std::cout << "shapes_lambda in : MaterialCohesiveDamageIntrinsic::computeLambdaOnQuad " << std::endl; // auto shapes_lambda = fem_lambda.getShapes(underlying_type,ghost_type,0); // ArrayPrintHelper::print_content(shapes_lambda,std::cout,0); -// std::cout << "shapes_cohesive in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl; +// std::cout << "shapes_cohesive in : MaterialCohesiveDamageIntrinsic::computeLambdaOnQuad " << std::endl; // auto & fem_cohesive = -// this->model->getFEEngineClass("CohesiveFEEngine"); +// model->getFEEngineClass("CohesiveFEEngine"); // auto shapes_cohesive = fem_cohesive.getShapes(type,ghost_type,0); // ArrayPrintHelper::print_content(shapes_cohesive,std::cout,0); } /* -------------------------------------------------------------------------- */ template -void MaterialCohesiveDamage::computeTraction(ElementType el_type, +void MaterialCohesiveDamageIntrinsic::computeTraction(ElementType el_type, GhostType ghost_type) { for (const auto & type : getElementFilter().elementTypes( spatial_dimension, ghost_type, _ek_cohesive)) { auto & elem_filter = getElementFilter(type, ghost_type); auto nb_element = elem_filter.size(); if (nb_element == 0) { continue; } /// compute normals @f$\mathbf{n}@f$ computeNormal(model->getCurrentPosition(), normals(type, ghost_type), type, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ computeOpening(model->getDisplacement(), opening(type, ghost_type), type, ghost_type); computeLambdaOnQuad(el_type, ghost_type); // auto & traction = tractions(type, ghost_type); -// std::cout << "traction in : MaterialCohesiveDamage::computeTraction " << std::endl; +// std::cout << "traction in : MaterialCohesiveDamageIntrinsic::computeTraction " << std::endl; // ArrayPrintHelper::print_content(traction,std::cout,0); // auto & lambda_ = lambda(type, ghost_type); -// std::cout << "lambda in : MaterialCohesiveDamage::computeTraction " << std::endl; +// std::cout << "lambda in : MaterialCohesiveDamageIntrinsic::computeTraction " << std::endl; // ArrayPrintHelper::print_content(lambda_,std::cout,0); +// lambda.printself(std::cout); for (auto && args : getArguments(el_type, ghost_type)) { - this->computeTractionOnQuad(args); + computeTractionOnQuad(args); } } } /* -------------------------------------------------------------------------- */ template -void MaterialCohesiveDamage::computeTangentTraction( +void MaterialCohesiveDamageIntrinsic::computeTangentTraction( ElementType el_type, Array & tangent_matrix_uu, Array & tangent_matrix_ll, GhostType ghost_type) { AKANTU_DEBUG_IN(); for (auto && [args, tangent_uu, tangent_ll] : zip(getArguments(el_type, ghost_type), make_view(tangent_matrix_uu), make_view(tangent_matrix_ll))) { computeTangentTractionOnQuad(tangent_uu, tangent_ll, args); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template class MaterialCohesiveDamage<1>; -template class MaterialCohesiveDamage<2>; -template class MaterialCohesiveDamage<3>; +template class MaterialCohesiveDamageIntrinsic<1>; +template class MaterialCohesiveDamageIntrinsic<2>; +template class MaterialCohesiveDamageIntrinsic<3>; const bool material_is_alocated_cohesive_damage [[maybe_unused]] = - instantiateMaterial("cohesive_damage"); + instantiateMaterial("cohesive_damage_intrinsic"); } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.hh similarity index 82% copy from src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh copy to src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.hh index 965cd4e85..ce3112be6 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.hh @@ -1,126 +1,117 @@ /** * Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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 "material_cohesive_damage.hh" /* -------------------------------------------------------------------------- */ -#ifndef AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_ -#define AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_ +#ifndef AKANTU_MATERIAL_COHESIVE_DAMAGE_INTRINSIC_HH_ +#define AKANTU_MATERIAL_COHESIVE_DAMAGE_INTRINSIC_HH_ namespace akantu { /** - * Cohesive material linear damage for extrinsic case + * Cohesive material linear damage for intrinsinc case * * parameters in the material files : * - sigma_c : critical stress sigma_c (default: 0) - * - beta : weighting parameter for sliding and normal opening (default: - * 0) - * - G_cI : fracture energy for mode I (default: 0) - * - G_cII : fracture energy for mode II (default: 0) - * - penalty : stiffness in compression to prevent penetration + * - G_c : fracture energy for mode I (default: 0) + * - k : cohesive stiffness */ -template class MaterialCohesiveDamage : public MaterialCohesive { +template class MaterialCohesiveDamageIntrinsic : public MaterialCohesiveDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - MaterialCohesiveDamage(SolidMechanicsModel & model, const ID & id = ""); + MaterialCohesiveDamageIntrinsic(SolidMechanicsModel & model, const ID & id = ""); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - /// initialize the material parameters - void initMaterial() override; - /// assemble stiffness void assembleStiffnessMatrix(GhostType ghost_type) override; /// assemble residual void assembleInternalForces(GhostType ghost_type = _not_ghost) override; protected: void computeLambdaOnQuad(ElementType type, GhostType ghost_type); inline decltype(auto) getArguments(ElementType element_type, GhostType ghost_type) { using namespace tuple; return zip_append( MaterialCohesive::getArguments(element_type, ghost_type), "lambda"_n = make_view(this->lambda(element_type, ghost_type)), + "czm_damage"_n = this->czm_damage(element_type, ghost_type), "err_opening"_n = make_view(this->err_openings(element_type, ghost_type))); } /// constitutive law void computeTraction(ElementType el_type, GhostType ghost_type = _not_ghost) override; /// compute tangent stiffness matrix /// WARNING : override was removed, not sure it should have void computeTangentTraction(ElementType el_type, Array & tangent_matrix_uu, Array & tangent_matrix_ll, GhostType ghost_type); /// compute the traction for a given quadrature point template inline void computeTractionOnQuad(Args && args); template inline void computeTangentTractionOnQuad(Eigen::MatrixBase & tangent_uu, Eigen::MatrixBase & tangent_ll, Args && args); [[nodiscard]] bool needLambda() const override { return true; } + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: - /// cohesive stiffness - Real k; - - /// mode I fracture energy - Real G_c; /// augmented lagrange multiplier CohesiveInternalField & lambda; /// target opening CohesiveInternalField & err_openings; /// cohesive damage CohesiveInternalField & czm_damage; Vector normal_opening; Real normal_opening_norm{0.}; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ } // namespace akantu -#include "material_cohesive_damage_inline_impl.hh" +#include "material_cohesive_damage_intrinsic_inline_impl.hh" -#endif /* AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_ */ +#endif /* AKANTU_MATERIAL_COHESIVE_DAMAGE_INTRINSIC_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_inline_impl.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic_inline_impl.hh similarity index 57% rename from src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_inline_impl.hh rename to src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic_inline_impl.hh index b2d9bfb24..d62cd4753 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_inline_impl.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic_inline_impl.hh @@ -1,78 +1,63 @@ /** * Copyright (©) 2015-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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_damage.hh" -#include "aka_static_if.hh" -#include "solid_mechanics_model_cohesive.hh" -/* -------------------------------------------------------------------------- */ +#include "material_cohesive_damage_intrinsic.hh" /* -------------------------------------------------------------------------- */ -#ifndef AKANTU_MATERIAL_COHESIVE_DAMAGE_INLINE_IMPL_HH_ -#define AKANTU_MATERIAL_COHESIVE_DAMAGE_INLINE_IMPL_HH_ +#ifndef AKANTU_MATERIAL_COHESIVE_DAMAGE_INTRINSIC_INLINE_IMPL_HH_ +#define AKANTU_MATERIAL_COHESIVE_DAMAGE_INTRINSIC_INLINE_IMPL_HH_ /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template template -inline void MaterialCohesiveDamage::computeTractionOnQuad(Args && args) { +inline void MaterialCohesiveDamageIntrinsic::computeTractionOnQuad(Args && args) { auto && lambda = args["lambda"_n]; + auto && d = args["czm_damage"_n]; auto && err_opening = args["err_opening"_n]; auto && opening = args["opening"_n]; auto && traction = args["traction"_n]; -// traction = lambda - (opening * k); - traction = lambda; - /// TODO : COMPUTE augmented_compliance -// Real d(0.0); -// Real augmented_compliance = d / k; -// err_opening = opening - lambda * augmented_compliance; - err_opening = opening; -// std::cout << "--- " << std::endl; -// std::cout << "k = " << k << std::endl; -// std::cout << "opening = " << opening << std::endl; -// std::cout << "lambda = " << lambda << std::endl; -// std::cout << "traction = " << traction << std::endl; -// std::cout << "err_opening = " << err_opening << std::endl; + traction = lambda - (opening * this->k); + Real A = this->augmented_compliance(d); + err_opening = opening - lambda * A; } /* -------------------------------------------------------------------------- */ template template -inline void MaterialCohesiveDamage::computeTangentTractionOnQuad( +inline void MaterialCohesiveDamageIntrinsic::computeTangentTractionOnQuad( Eigen::MatrixBase & tangent_uu, - Eigen::MatrixBase & tangent_ll, Args && /*args*/) { - /// TODO : COMPUTE augmented_compliance + Eigen::MatrixBase & tangent_ll, Args && args) { /// TODO : check basis - Real d(0.0); -// Real augmented_compliance = d / k; - for (Int i = 0; i < dim; ++i) { -// tangent_uu(i, i) = -k; -// tangent_ll(i, i) = -augmented_compliance; - } + auto && d = args["czm_damage"_n]; + Real A = this->augmented_compliance(d); + tangent_uu = Eigen::Matrix::Identity()*(-this->k); + tangent_ll = Eigen::Matrix::Identity()*(-A); } /* -------------------------------------------------------------------------- */ } // namespace akantu /* -------------------------------------------------------------------------- */ -#endif // AKANTU_MATERIAL_COHESIVE_DAMAGE_INLINE_IMPL_HH_ +#endif // AKANTU_MATERIAL_COHESIVE_DAMAGE_INTRINSIC_INLINE_IMPL_HH_ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc index 023f2f468..aa5a9889d 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc @@ -1,808 +1,873 @@ /** * Copyright (©) 2012-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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 "solid_mechanics_model_cohesive.hh" #include "aka_iterators.hh" #include "cohesive_element_inserter.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "fe_engine_template.hh" #include "global_ids_updater.hh" #include "integrator_gauss.hh" #include "material_cohesive.hh" #include "mesh_accessor.hh" #include "mesh_global_data_updater.hh" #include "parser.hh" #include "shape_cohesive.hh" /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class CohesiveMeshGlobalDataUpdater : public MeshGlobalDataUpdater { public: CohesiveMeshGlobalDataUpdater(SolidMechanicsModelCohesive & model) : model(model), mesh(model.getMesh()), global_ids_updater(model.getMesh(), model.cohesive_synchronizer.get()) { } /* ------------------------------------------------------------------------ */ std::tuple updateData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) override { auto * cohesive_nodes_event = dynamic_cast(&nodes_event); if (cohesive_nodes_event == nullptr) { return std::make_tuple(nodes_event.getList().size(), elements_event.getList().size()); } /// update nodes type auto & new_nodes = cohesive_nodes_event->getList(); auto & old_nodes = cohesive_nodes_event->getOldNodesList(); auto local_nb_new_nodes = new_nodes.size(); auto nb_new_nodes = local_nb_new_nodes; if (mesh.isDistributed()) { MeshAccessor mesh_accessor(mesh); auto & nodes_flags = mesh_accessor.getNodesFlags(); auto nb_old_nodes = nodes_flags.size(); nodes_flags.resize(nb_old_nodes + local_nb_new_nodes); for (auto && [old_node, new_node] : zip(old_nodes, new_nodes)) { nodes_flags(new_node) = nodes_flags(old_node); } model.updateCohesiveSynchronizers(elements_event); nb_new_nodes = global_ids_updater.updateGlobalIDs(new_nodes.size()); } auto nb_new_elements = elements_event.getList().size(); const auto & comm = mesh.getCommunicator(); comm.allReduce(nb_new_elements, SynchronizerOperation::_sum); if (nb_new_elements > 0) { mesh.sendEvent(elements_event); } if (nb_new_nodes > 0) { mesh.sendEvent(nodes_event); } return std::make_tuple(nb_new_nodes, nb_new_elements); } private: SolidMechanicsModelCohesive & model; Mesh & mesh; GlobalIdsUpdater global_ids_updater; }; /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive( Mesh & mesh, Int dim, const ID & id, const std::shared_ptr & dof_manager) : SolidMechanicsModel(mesh, dim, id, dof_manager, ModelType::_solid_mechanics_model_cohesive), tangents("tangents", id), facet_stress("facet_stress", id), facet_material("facet_material", id) { AKANTU_DEBUG_IN(); registerFEEngineObject("CohesiveFEEngine", mesh, Model::spatial_dimension); auto && tmp_material_selector = std::make_shared(*this); tmp_material_selector->setFallback(this->getConstitutiveLawSelector()); this->setConstitutiveLawSelector(tmp_material_selector); this->mesh.registerDumper("cohesive elements", id); this->mesh.addDumpMeshToDumper("cohesive elements", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); if (this->mesh.isDistributed()) { /// create the distributed synchronizer for cohesive elements this->cohesive_synchronizer = std::make_unique( mesh, "cohesive_distributed_synchronizer"); auto & synchronizer = mesh.getElementSynchronizer(); this->cohesive_synchronizer->split(synchronizer, [](auto && el) { return Mesh::getKind(el.type) == _ek_cohesive; }); this->registerSynchronizer(*cohesive_synchronizer, SynchronizationTag::_constitutive_law_id); this->registerSynchronizer(*cohesive_synchronizer, SynchronizationTag::_smm_stress); this->registerSynchronizer(*cohesive_synchronizer, SynchronizationTag::_smm_boundary); } this->inserter = std::make_unique( this->mesh, id + ":cohesive_element_inserter"); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::setTimeStep(Real time_step, const ID & solver_id) { SolidMechanicsModel::setTimeStep(time_step, solver_id); this->mesh.getDumper("cohesive elements").setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); const auto & smmc_options = aka::as_type(options); this->is_extrinsic = smmc_options.is_extrinsic; inserter->setIsExtrinsic(is_extrinsic); if (mesh.isDistributed()) { auto & mesh_facets = inserter->getMeshFacets(); auto & synchronizer = aka::as_type(mesh_facets.getElementSynchronizer()); // synchronizeGhostFacetsConnectivity(); /// create the facet synchronizer for extrinsic simulations if (is_extrinsic) { facet_stress_synchronizer = std::make_unique( synchronizer, id + ":facet_stress_synchronizer"); facet_stress_synchronizer->swapSendRecv(); this->registerSynchronizer(*facet_stress_synchronizer, SynchronizationTag::_smmc_facets_stress); } } MeshAccessor mesh_accessor(mesh); mesh_accessor.registerGlobalDataUpdater( std::make_unique(*this)); auto && [section, is_empty] = this->getParserSection(); if (not is_empty) { auto inserter_section = section.getSubSections(ParserType::_cohesive_inserter); if (inserter_section.begin() != inserter_section.end()) { inserter->parseSection(*inserter_section.begin()); } } SolidMechanicsModel::initFullImpl(options); AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initConstitutiveLaws() { AKANTU_DEBUG_IN(); // make sure the material are instantiated instantiateMaterials(); auto & material_selector = this->getConstitutiveLawSelector(); // set the facet information in the material in case of dynamic insertion // to know what material to call for stress checks const Mesh & mesh_facets = inserter->getMeshFacets(); facet_material.initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _default_value = DefaultMaterialCohesiveSelector::getDefaultCohesiveMaterial(*this)); for_each_element( mesh_facets, [&](auto && element) { auto mat_index = material_selector(element); if (not mat_index) { return; } auto & mat = aka::as_type(this->getConstitutiveLaw(mat_index)); facet_material(element) = mat_index; if (is_extrinsic) { mat.addFacet(element); } }, _spatial_dimension = spatial_dimension - 1, _ghost_type = _not_ghost); SolidMechanicsModel::initConstitutiveLaws(); auto & initial_nodes = mesh.getNodalData("initial_nodes_match"); initial_nodes.resize(mesh.getNbNodes()); for (auto && [node, init_node] : enumerate(initial_nodes)) { init_node = node; } auto need_lambda = false; this->for_each_constitutive_law([&need_lambda](auto && material) { if (aka::is_of_type(material)) { need_lambda |= aka::as_type(material).needLambda(); } }); if (need_lambda) { lambda = std::make_unique>(0, spatial_dimension, "cohesive lambda"); previous_lambda = std::make_unique>(0, spatial_dimension, "previous lambda"); lambda_increment = std::make_unique>(0, spatial_dimension, "lambda increment"); lambda_blocked_dofs = std::make_unique>(0, spatial_dimension, "lambda_blocked_dofs"); lambda_mesh = std::make_unique(spatial_dimension, mesh.getCommunicator(), "lambda_mesh"); registerFEEngineObject("LambdaFEEngine", *lambda_mesh, Model::spatial_dimension - 1); auto & nodes_to_lambda = mesh.getNodalData("nodes_to_lambda"); nodes_to_lambda.resize(mesh.getNbNodes(), -1); } if (is_extrinsic) { this->initAutomaticInsertion(); } else { this->insertIntrinsicElements(); } if (lambda) { auto & dof_manager = this->getDOFManager(); if (!dof_manager.hasDOFs("lambda")) { dof_manager.registerDOFs("lambda", *this->lambda, *lambda_mesh); dof_manager.registerDOFsIncrement("lambda", *this->lambda_increment); dof_manager.registerDOFsPrevious("lambda", *this->previous_lambda); dof_manager.registerBlockedDOFs("lambda", *this->lambda_blocked_dofs); } } AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModelCohesive::initModel() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initModel(); /// add cohesive type connectivity ElementType type = _not_defined; for (auto && type_ghost : ghost_types) { for (const auto & tmp_type : mesh.elementTypes(spatial_dimension, type_ghost)) { const auto & connectivity = mesh.getConnectivity(tmp_type, type_ghost); if (connectivity.empty()) { continue; } type = tmp_type; auto type_facet = Mesh::getFacetType(type); auto type_cohesive = FEEngine::getCohesiveElementType(type_facet); AKANTU_DEBUG_ASSERT(Mesh::getKind(type_cohesive) == _ek_cohesive, "The element type " << type_cohesive << " is not a cohesive type"); mesh.addConnectivityType(type_cohesive, type_ghost); } } AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertIntrinsicElements() { AKANTU_DEBUG_IN(); inserter->insertIntrinsicElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); this->resizeFacetStress(); /// compute normals on facets this->computeNormals(); this->initStressInterpolation(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateAutomaticInsertion() { AKANTU_DEBUG_IN(); this->inserter->limitCheckFacets(); this->updateFacetStressSynchronizer(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initStressInterpolation() { Mesh & mesh_facets = inserter->getMeshFacets(); /// compute quadrature points coordinates on facets Array & position = mesh.getNodes(); ElementTypeMapArray quad_facets("quad_facets", id); quad_facets.initialize(mesh_facets, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension - 1); getFEEngine("FacetsFEEngine") .interpolateOnIntegrationPoints(position, quad_facets); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure ElementTypeMapArray elements_quad_facets("elements_quad_facets", id); elements_quad_facets.initialize( mesh, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension); for (auto elem_gt : ghost_types) { for (const auto & type : mesh.elementTypes(Model::spatial_dimension, elem_gt)) { auto nb_element = mesh.getNbElement(type, elem_gt); if (nb_element == 0) { continue; } /// compute elements' quadrature points and list of facet /// quadrature points positions by element const auto & facet_to_element = mesh_facets.getSubelementToElement(type, elem_gt); auto & el_q_facet = elements_quad_facets(type, elem_gt); auto facet_type = Mesh::getFacetType(type); auto nb_quad_per_facet = getFEEngine("FacetsFEEngine").getNbIntegrationPoints(facet_type); auto nb_facet_per_elem = facet_to_element.getNbComponent(); // small hack in the loop to skip boundary elements, they are silently // initialized to NaN to see if this causes problems el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet, std::numeric_limits::quiet_NaN()); for (auto && data : zip(make_view(facet_to_element), make_view(el_q_facet, spatial_dimension, nb_quad_per_facet))) { const auto & global_facet = std::get<0>(data); auto & el_q = std::get<1>(data); if (global_facet == ElementNull) { continue; } auto && quad_f = make_view(quad_facets(global_facet.type, global_facet.ghost_type), spatial_dimension, nb_quad_per_facet) .begin()[global_facet.element]; el_q = quad_f; } } } /// loop over non cohesive materials this->for_each_constitutive_law([&](auto && material) { if (aka::is_of_type(material)) { return; } /// initialize the interpolation function material.initElementalFieldInterpolation(elements_quad_facets); }); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::assembleInternalForces() { AKANTU_DEBUG_IN(); // f_int += f_int_cohe this->for_each_constitutive_law([&](auto && material) { try { auto & mat = aka::as_type(material); mat.computeTraction(_not_ghost); } catch (std::bad_cast & bce) { } }); // ArrayPrintHelper::ArrayPrintHelper::print_content(*(this->blocked_dofs),std::cout,0); SolidMechanicsModel::assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeNormals() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = this->inserter->getMeshFacets(); this->getFEEngine("FacetsFEEngine") .computeNormalsOnIntegrationPoints(_not_ghost); /** * @todo store tangents while computing normals instead of * recomputing them as follows: */ /* ------------------------------------------------------------------------ */ UInt tangent_components = Model::spatial_dimension * (Model::spatial_dimension - 1); tangents.initialize(mesh_facets, _nb_component = tangent_components, _spatial_dimension = Model::spatial_dimension - 1); for (auto facet_type : mesh_facets.elementTypes(Model::spatial_dimension - 1)) { const Array & normals = this->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(facet_type); Array & tangents = this->tangents(facet_type); Math::compute_tangents(normals, tangents); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::interpolateStress() { ElementTypeMapArray by_elem_result("temporary_stress_by_facets", id); this->for_each_constitutive_law([&](auto && material) { if (not aka::is_of_type(material)) { /// interpolate stress on facet quadrature points positions material.interpolateStressOnFacets(facet_stress, by_elem_result); } }); this->synchronize(SynchronizationTag::_smmc_facets_stress); } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::checkCohesiveStress() { AKANTU_DEBUG_IN(); if (not is_extrinsic) { AKANTU_EXCEPTION( "This function can only be used for extrinsic cohesive elements"); } interpolateStress(); this->for_each_constitutive_law([&](auto && material) { if (aka::is_of_type(material)) { /// check which not ghost cohesive elements are to be created auto & mat_cohesive = aka::as_type(material); mat_cohesive.checkInsertion(); } }); /// insert cohesive elements auto nb_new_elements = inserter->insertElements(); AKANTU_DEBUG_OUT(); return nb_new_elements; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateLambdaMesh() { +// std::cout << " SolidMechanicsModelCohesive::updateLambdaMesh " << std::endl; + const auto & connectivities = mesh.getConnectivities(); const auto & initial_nodes = mesh.getNodalData("initial_nodes_match"); auto & nodes_to_lambda = mesh.getNodalData("nodes_to_lambda"); nodes_to_lambda.resize(mesh.getNbNodes(), -1); auto & lambda_connectivities = MeshAccessor(*lambda_mesh).getConnectivities(); auto lambda_id = lambda->size(); std::vector new_nodes; NewNodesEvent node_event(*lambda_mesh, AKANTU_CURRENT_FUNCTION); auto & node_list = node_event.getList(); +// std::cout << "node_list before = " << std::endl ; +// ArrayPrintHelper::print_content(node_list,std::cout,0); NewElementsEvent element_event(*lambda_mesh, AKANTU_CURRENT_FUNCTION); auto & element_list = element_event.getList(); +// std::cout << "element_list before = " << std::endl ; +// ArrayPrintHelper::print_content(element_list,std::cout,0); for (auto ghost_type : ghost_types) { for (auto type : connectivities.elementTypes(_element_kind = _ek_cohesive, _ghost_type = ghost_type)) { auto underlying_type = Mesh::getFacetType(type); auto & connectivity = connectivities(type, ghost_type); auto & lambda_connectivity = lambda_connectivities(underlying_type, ghost_type); +// std::cout << "connectivity = " << std::endl ; +// ArrayPrintHelper::print_content(connectivity,std::cout,0); + +// std::cout << "lambda_connectivity = " << std::endl ; +// ArrayPrintHelper::print_content(lambda_connectivity,std::cout,0); + + auto nb_new_elements = connectivity.size() - lambda_connectivity.size(); auto nb_old_elements = lambda_connectivity.size(); lambda_connectivity.resize(connectivity.size()); auto && view_iconn = make_view(connectivity, connectivity.getNbComponent()); auto && view_conn = make_view(lambda_connectivity, lambda_connectivity.getNbComponent()); +// std::cout << " connectivity.getNbComponent() " << connectivity.getNbComponent() << std::endl ; +// std::cout << " lambda_connectivity.getNbComponent() " << lambda_connectivity.getNbComponent() << std::endl ; + +// for (auto && [conn, lambda_conn] : +// zip(range(view_iconn.begin() + nb_old_elements, view_iconn.end()), +// range(view_conn.begin() + nb_old_elements, view_conn.end()))) { +// std::cout << "conn = " << conn << std::endl ; +// std::cout << "lambda_conn = " << lambda_conn << std::endl ; + +// for (auto && [n, lambda_node] : enumerate(lambda_conn)) { +// std::cout << "n = " << n << std::endl ; +// auto node = conn[n]; +// std::cout << "node conn= " << node << std::endl ; + +// node = initial_nodes(node); +// std::cout << "node initial_nodes = " << node << std::endl ; + +// auto & ntl = nodes_to_lambda(node); +// std::cout << "ntl = " << ntl << std::endl ; + +// if (ntl == -1) { +// ntl = lambda_id; +// new_nodes.push_back(node); +// node_list.push_back(lambda_id); +// ++lambda_id; +// } + +// lambda_node = ntl; +// } +// } + for (auto && [conn, lambda_conn] : zip(range(view_iconn.begin() + nb_old_elements, view_iconn.end()), range(view_conn.begin() + nb_old_elements, view_conn.end()))) { +// std::cout << "conn = " << conn << std::endl ; +// std::cout << "lambda_conn = " << lambda_conn << std::endl ; + for (auto && [n, lambda_node] : enumerate(lambda_conn)) { +// std::cout << "n = " << n << std::endl ; auto node = conn[n]; +// std::cout << "node conn= " << node << std::endl ; + node = initial_nodes(node); +// std::cout << "node initial_nodes = " << node << std::endl ; + auto & ntl = nodes_to_lambda(node); - if (ntl == -1) { +// std::cout << "ntl = " << ntl << std::endl ; + +// if (ntl == -1) { ntl = lambda_id; new_nodes.push_back(node); node_list.push_back(lambda_id); ++lambda_id; - } +// } lambda_node = ntl; } } for (auto && el : arange(nb_old_elements, nb_new_elements)) { +// std::cout << "el = " << el << std::endl ; element_list.push_back({underlying_type, el, ghost_type}); } } } + +// std::cout << "node_list after = " << std::endl ; +// ArrayPrintHelper::print_content(node_list,std::cout,0); + +// std::cout << "element_list after = " << std::endl ; +// ArrayPrintHelper::print_content(element_list,std::cout,0); + lambda_mesh->copyNodes(mesh, new_nodes); lambda->resize(lambda_id, 0.); lambda_increment->resize(lambda_id, 0.); previous_lambda->resize(lambda_id, 0.); lambda_mesh->sendEvent(element_event); lambda_mesh->sendEvent(node_event); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); +// std::cout << " SolidMechanicsModelCohesive::onElementsAdded " << std::endl; + SolidMechanicsModel::onElementsAdded(element_list, event); if (lambda) { auto & lambda_connectivities = MeshAccessor(*lambda_mesh).getConnectivities(); for (auto ghost_type : ghost_types) { for (auto type : mesh.elementTypes(_element_kind = _ek_cohesive, _ghost_type = ghost_type)) { auto underlying_type = Mesh::getFacetType(type); if (not lambda_connectivities.exists(underlying_type, ghost_type)) { auto nb_nodes_per_element = Mesh::getNbNodesPerElement(underlying_type); lambda_connectivities.alloc(0, nb_nodes_per_element, underlying_type, ghost_type, -1); } } } } if (is_extrinsic) { resizeFacetStress(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array & new_nodes, const NewNodesEvent & event) { SolidMechanicsModel::onNodesAdded(new_nodes, event); auto & initial_nodes = mesh.getNodalData("initial_nodes_match"); + auto old_max_nodes = initial_nodes.size(); initial_nodes.resize(mesh.getNbNodes(), -1); const auto * cohesive_event = dynamic_cast(&event); if (cohesive_event == nullptr) { for (auto && node : new_nodes) { initial_nodes(node) = node; } return; } auto old_nodes = cohesive_event->getOldNodesList(); // getting a corrected old_nodes for multiple doubling for (auto & onode : old_nodes) { while (onode >= old_max_nodes) { auto nnode = new_nodes.find(onode); AKANTU_DEBUG_ASSERT(nnode != -1, "If the old node is bigger than old_max_nodes it " "should also be a new node"); onode = nnode; } } auto copy = [&new_nodes, &old_nodes](auto & arr) { auto it = make_view(arr, arr.getNbComponent()).begin(); for (auto [new_node, old_node] : zip(new_nodes, old_nodes)) { it[new_node] = it[old_node]; } }; for (auto && ptr : {displacement.get(), velocity.get(), acceleration.get(), current_position.get(), previous_displacement.get(), displacement_increment.get()}) { if (ptr != nullptr) { copy(*ptr); } } copy(*blocked_dofs); copy(getDOFManager().getSolution("displacement")); copy(initial_nodes); // correct connectivities if (lambda) { lambda_blocked_dofs->resize(mesh.getNbNodes(), false); copy(*lambda_blocked_dofs); updateLambdaMesh(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::afterSolveStep(bool converged) { /* * This is required because the Cauchy stress is the stress measure that * is used to check the insertion of cohesive elements */ if (converged) { this->for_each_constitutive_law([](auto && mat) { if (mat.isFiniteDeformation()) { mat.computeAllCauchyStresses(_not_ghost); } }); } SolidMechanicsModel::afterSolveStep(converged); } /* -------------------------------------------------------------------------- */ ModelSolverOptions SolidMechanicsModelCohesive::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options = SolidMechanicsModel::getDefaultSolverOptions(type); - // if (lambda) { - if (true) { + if (lambda) { +// if (true) { switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["lambda"] = IntegrationSchemeType::_central_difference; options.solution_type["lambda"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["lambda"] = IntegrationSchemeType::_pseudo_time; options.solution_type["lambda"] = IntegrationScheme::_not_defined; break; } case TimeStepSolverType::_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["lambda"] = IntegrationSchemeType::_central_difference; options.solution_type["lambda"] = IntegrationScheme::_acceleration; } else { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["lambda"] = IntegrationSchemeType::_trapezoidal_rule_2; options.solution_type["lambda"] = IntegrationScheme::_displacement; } break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } } return options; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "SolidMechanicsModelCohesive [" << "\n"; SolidMechanicsModel::printself(stream, indent + 2); stream << space << "]\n"; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::resizeFacetStress() { AKANTU_DEBUG_IN(); this->facet_stress.initialize(getFEEngine("FacetsFEEngine"), _nb_component = 2 * spatial_dimension * spatial_dimension, _spatial_dimension = spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::addDumpGroupFieldToDumper( const std::string & dumper_name, const std::string & field_id, const std::string & group_name, ElementKind element_kind, bool padding_flag) { AKANTU_DEBUG_IN(); Int spatial_dimension = Model::spatial_dimension; ElementKind _element_kind = element_kind; if (dumper_name == "cohesive elements") { _element_kind = _ek_cohesive; } else if (dumper_name == "facets") { spatial_dimension = Model::spatial_dimension - 1; } SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, spatial_dimension, _element_kind, padding_flag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onDump() { this->flattenAllRegisteredInternals(_ek_cohesive); SolidMechanicsModel::onDump(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.hh index aadbdd329..b6fe43f03 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.hh @@ -1,324 +1,346 @@ /** * Copyright (©) 2012-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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" #include "material_selector_cohesive.hh" #include "random_internal_field.hh" // included to have the specialization of // ParameterTyped::operator Real() #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH_ #define AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH_ /* -------------------------------------------------------------------------- */ namespace akantu { class FacetSynchronizer; class FacetStressSynchronizer; class ElementSynchronizer; } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ struct FacetsCohesiveIntegrationOrderFunctor { template ::cohesive_type> struct _helper { static constexpr int get() { - return ElementClassProperty::polynomial_degree; +// return 1; + return ElementClassProperty::polynomial_degree; } }; template struct _helper { static constexpr int get() { - return ElementClassProperty::polynomial_degree; - } + // return 1; + return ElementClassProperty::polynomial_degree; } + }; + + template static inline constexpr int getOrder() { + return _helper::get(); + } +}; + +/* -------------------------------------------------------------------------- */ +struct CohesiveIntegrationOrderFunctor { + template ::cohesive_type> + struct _helper { + static constexpr int get() { + // return 1; + return ElementClassProperty::polynomial_degree; } + }; + + template struct _helper { + static constexpr int get() { + // return 1; + return ElementClassProperty::polynomial_degree; } }; template static inline constexpr int getOrder() { return _helper::get(); } }; /* -------------------------------------------------------------------------- */ /* Solid Mechanics Model for Cohesive elements */ /* -------------------------------------------------------------------------- */ class SolidMechanicsModelCohesive : public SolidMechanicsModel, public SolidMechanicsModelEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewCohesiveNodesEvent : public NewNodesEvent { public: AKANTU_GET_MACRO_NOT_CONST(OldNodesList, old_nodes, Array &); AKANTU_GET_MACRO(OldNodesList, old_nodes, const Array &); protected: Array old_nodes; }; using MyFEEngineCohesiveType = FEEngineTemplate; using MyFEEngineFacetType = FEEngineTemplate; using MyFEEngineLambdaType = FEEngineTemplate; SolidMechanicsModelCohesive( Mesh & mesh, Int dim = _all_dimensions, const ID & id = "solid_mechanics_model_cohesive", const std::shared_ptr & dof_manager = nullptr); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize the cohesive model void initFullImpl(const ModelOptions & options) override; public: /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// assemble the residual for the explicit scheme void assembleInternalForces() override; /// function to perform a stress check on each facet and insert /// cohesive elements if needed (returns the number of new cohesive /// elements) UInt checkCohesiveStress(); /// interpolate stress on facets void interpolateStress(); /// update automatic insertion after a change in the element inserter void updateAutomaticInsertion(); /// insert intrinsic cohesive elements void insertIntrinsicElements(); // template bool solveStepCohesive(Real tolerance, Real & error, UInt // max_iteration = 100, // bool load_reduction = false, // Real tol_increase_factor = 1.0, // bool do_not_factorize = false); protected: /// initialize stress interpolation void initStressInterpolation(); /// initialize the model void initModel() override; /// initialize cohesive material void initConstitutiveLaws() override; /// init facet filters for cohesive materials void initFacetFilter(); /// function to print the contain of the class void printself(std::ostream & stream, int indent = 0) const override; private: /// insert cohesive elements along a given physical surface of the mesh void insertElementsFromMeshData(const std::string & physical_name); /// initialize completely the model for extrinsic elements void initAutomaticInsertion(); /// compute facets' normals void computeNormals(); /// resize facet stress void resizeFacetStress(); /// init facets_check array void initFacetsCheck(); /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & new_nodes, const NewNodesEvent & event) override; void onElementsAdded(const Array & element_list, const NewElementsEvent & event) override; /* ------------------------------------------------------------------------ */ /* SolidMechanicsModelEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: void afterSolveStep(bool converged = true) override; [[nodiscard]] ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: void onDump() override; void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, ElementKind element_kind, bool padding_flag) override; protected: void updateCohesiveSynchronizers(NewElementsEvent & elements_event); void updateFacetStressSynchronizer(); friend class CohesiveElementInserter; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: [[nodiscard]] Int getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; protected: [[nodiscard]] Int getNbQuadsForFacetCheck(const Array & elements) const; template void packFacetStressDataHelper(const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements) const; template void unpackFacetStressDataHelper(ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements) const; template void packUnpackFacetStressDataHelper( std::conditional_t, ElementTypeMapArray> & data_to_pack, CommunicationBuffer & buffer, const Array & element) const; void updateLambdaMesh(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get facet mesh AKANTU_GET_MACRO_AUTO(MeshFacets, mesh.getMeshFacets()); /// get stress on facets vector AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(StressOnFacets, facet_stress, Real); /// get facet material AKANTU_GET_MACRO_BY_ELEMENT_TYPE(FacetMaterial, facet_material, Idx); /// get facet material AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FacetMaterial, facet_material, Idx); /// get facet material AKANTU_GET_MACRO_AUTO(FacetMaterial, facet_material); /// @todo THIS HAS TO BE CHANGED AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Tangents, tangents, Real); /// get element inserter AKANTU_GET_MACRO_NOT_CONST(ElementInserter, *inserter, CohesiveElementInserter &); /// get is_extrinsic boolean AKANTU_GET_MACRO(IsExtrinsic, is_extrinsic, bool); /// get cohesive elements synchronizer AKANTU_GET_MACRO_NOT_CONST(CohesiveSynchronizer, *cohesive_synchronizer, ElementSynchronizer &); - /// get the SolidMechanicsModel::displacement array + /// get the SolidMechanicsModel::lambda array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Lambda, lambda); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR(Lambda, lambda); /// get lambda mesh const Mesh & getLambdaMesh() { if (not lambda_mesh) { AKANTU_EXCEPTION("This model does not have a lambda mesh"); } return *lambda_mesh; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: friend class CohesiveMeshGlobalDataUpdater; /// lambda array std::unique_ptr> lambda; /// lambda array at the previous time step std::unique_ptr> previous_lambda; /// increment of lambda std::unique_ptr> lambda_increment; /// array specifing if a lambda degree of freedom is blocked or not std::unique_ptr> lambda_blocked_dofs; /// @todo store tangents when normals are computed: ElementTypeMapArray tangents; /// stress on facets on the two sides by quadrature point ElementTypeMapArray facet_stress; /// material to use if a cohesive element is created on a facet ElementTypeMapArray facet_material; bool is_extrinsic{false}; /// cohesive element inserter std::unique_ptr inserter; /// facet stress synchronizer std::unique_ptr facet_stress_synchronizer; /// cohesive elements synchronizer std::unique_ptr cohesive_synchronizer; std::unique_ptr lambda_mesh; }; } // namespace akantu #include "solid_mechanics_model_cohesive_inline_impl.hh" #endif /* AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH_ */