diff --git a/examples/c++/solid_mechanics_cohesive_model/cohesive_damage_extrinsic/cohesive_damage_extrinsic.cc b/examples/c++/solid_mechanics_cohesive_model/cohesive_damage_extrinsic/cohesive_damage_extrinsic.cc
index c21fe8761..4dfaa7ae2 100644
--- a/examples/c++/solid_mechanics_cohesive_model/cohesive_damage_extrinsic/cohesive_damage_extrinsic.cc
+++ b/examples/c++/solid_mechanics_cohesive_model/cohesive_damage_extrinsic/cohesive_damage_extrinsic.cc
@@ -1,112 +1,113 @@
/**
* 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 <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#include "element_group.hh"
#include "mesh_iterators.hh"
#include "non_linear_solver.hh"
-#include "solid_mechanics_model_cohesive.hh"
+#include "solid_mechanics_model_cohesive_damage.hh"
/* -------------------------------------------------------------------------- */
#include <iostream>
/* -------------------------------------------------------------------------- */
using namespace akantu;
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
const Int spatial_dimension = 2;
const Int max_steps = 1;
Real L = 0.4;
Real eps0 = 1e-1;
Mesh mesh(spatial_dimension);
mesh.read("triangle.msh");
- SolidMechanicsModelCohesive model(mesh);
+ SolidMechanicsModelCohesiveDamage model(mesh);
/// model initialization
model.initFull(_analysis_method = _explicit_lumped_mass
,_is_extrinsic = true);
Real time_step = model.getStableTimeStep() * 0.05;
model.setTimeStep(time_step);
std::cout << "Time step: " << time_step << std::endl;
model.getElementInserter().setLimit(_x, 0.5*L-0.1, 0.5*L+0.1);
model.updateAutomaticInsertion();
model.applyBC(BC::Dirichlet::FixedValue(0., _x), "left");
model.applyBC(BC::Dirichlet::FixedValue(0., _y), "point");
Int nb_nodes = mesh.getNbNodes();
Array<Real> & position = mesh.getNodes();
Array<Real> & velocity = model.getVelocity();
for (Int n = 0; n < nb_nodes; ++n) {
velocity(n,0) = eps0*position(n,0);
}
model.setBaseName("cohesive");
model.addDumpFieldVector("displacement");
// model.addDumpField("velocity");
// model.addDumpField("acceleration");
model.addDumpField("stress");
model.addDumpField("grad_u");
model.addDumpField("external_force");
model.addDumpField("internal_force");
model.dump();
/// Main loop
for (Int s = 1; s <= max_steps; ++s) {
model.applyBC(BC::Dirichlet::IncrementValue(time_step*eps0*L, _x), "right");
// debug::setDebugLevel(dblDump);
debug::setDebugLevel(dblInfo);
model.checkCohesiveInsertion(); /// TODO : replace by a function which does not interpolate stress
model.solveStep();
model.computeLagrangeMultiplier();
+ model.computeDamage();
model.checkCohesiveStress();
debug::setDebugLevel(dblWarning);
if (s % 1 == 0) {
model.dump();
std::cout << "passing step " << s << "/" << max_steps << std::endl;
}
}
Real E = model.getEnergy("potential");
std::cout << "Elastic energy : " << E << std::endl;
// if (Ed < Edt * 0.999 || Ed > Edt * 1.001 || std::isnan(Ed)) {
// std::cout << "The dissipated energy is incorrect" << std::endl;
// return EXIT_FAILURE;
// }
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/c++/solid_mechanics_cohesive_model/cohesive_damage_extrinsic/cohesive_damage_extrinsic.py b/examples/c++/solid_mechanics_cohesive_model/cohesive_damage_extrinsic/cohesive_damage_extrinsic.py
index f23324eb0..4ef34df3d 100644
--- a/examples/c++/solid_mechanics_cohesive_model/cohesive_damage_extrinsic/cohesive_damage_extrinsic.py
+++ b/examples/c++/solid_mechanics_cohesive_model/cohesive_damage_extrinsic/cohesive_damage_extrinsic.py
@@ -1,203 +1,204 @@
#!/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/developpement/akantu/build/buildCZMDamageDebug/python")
print(sys.path)
import akantu as aka
import numpy as np
import matplotlib.pyplot as plt
from czm_damage import *
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")
class FixedDisplacement (aka.DirichletFunctor):
'''
Fix the displacement at its current value
'''
def __init__(self, axis, vel):
super().__init__(axis)
self.axis = axis
self.time = 0
self.vel = vel
def set_time(self, t):
self.time = t
def get_imposed_disp(self):
return self.vel*self.time
def __call__(self, node, flags, disp, coord):
# sets the blocked dofs vector to true in the desired axis
flags[int(self.axis)] = True
disp[int(self.axis)] = self.get_imposed_disp()
def solve(material_file, mesh_file):
aka.parseInput(material_file)
spatial_dimension = 2
L = 0.4
# -------------------------------------------------------------------------
# Initialization
# -------------------------------------------------------------------------
mesh = aka.Mesh(spatial_dimension)
mesh.read(mesh_file)
- model = aka.SolidMechanicsModelCohesive(mesh)
+ model = aka.SolidMechanicsModelCohesiveDamage(mesh)
model.getElementInserter().setLimit(aka._x, 0.5*L-0.1, 0.5*L+0.1);
model.initFull(_analysis_method=aka._static, _is_extrinsic=True)
# Initilize a new solver (explicit Newmark with lumped mass)
model.initNewSolver(aka._explicit_lumped_mass)
# Dynamic insertion of cohesive elements
model.updateAutomaticInsertion()
set_dumpers(model)
E = model.getMaterial("steel").getReal("E")
Gc = model.getMaterial("cohesive").getReal("G_c")
sigc = model.getMaterial("cohesive").getReal("sigma_c")
wc = 2.*Gc/sigc
w0 = (sigc/E)*L
k = model.getMaterial("cohesive").getReal("k")
g = degradation_function_linear_czm()
h = softening_function_linear_czm(k,sigc,Gc)
updater = damage_updater(Gc,sigc,k,g,h)
imposed_disp = [0.]
reaction_force = [0.]
E_pot = []
E_kin = []
E_dis = []
E_rev = []
E_con = []
U = 0.
F = 0.
# -------------------------------------------------------------------------
# Initial and boundary conditions
# -------------------------------------------------------------------------
eps0dot = 5e-1
tc = w0/(eps0dot*L)
functor_r = FixedDisplacement(aka._x, L*eps0dot)
functor_r.set_time(tc)
model.applyBC(functor_r, 'right')
model.applyBC(aka.FixedValue(0.0, aka._x), "left")
model.applyBC(aka.FixedValue(0.0, aka._y), "point")
nodes = model.getMesh().getNodes()
disp_field = np.zeros(nodes.shape)
disp_field[:, 0] = nodes[:, 0]*(w0/L)
model.getDisplacement()[:] = disp_field
vel_field = np.zeros(nodes.shape)
vel_field[:, 0] = nodes[:, 0]*eps0dot
model.getVelocity()[:] = vel_field
model.getExternalForce()[:] = 0
model.assembleInternalForces()
d = model.getMaterial("cohesive").getInternalReal("czm_damage")(aka._segment_2)
lda = model.getMaterial("cohesive").getInternalReal("lambda")(aka._segment_2)
Fint = model.getInternalForce()
nodes_right = mesh.getElementGroup("right").getNodeGroup().getNodes()
U = functor_r.get_imposed_disp()
imposed_disp.append(U)
F = -np.sum(Fint[nodes_right,0])
reaction_force.append(F)
Ep = model.getEnergy("potential")
Ek = model.getEnergy("kinetic")
Ed = model.getEnergy("dissipated")
Er = model.getEnergy("reversible")
E_pot.append(Ep)
E_kin.append(Ek)
E_dis.append(Ed)
E_rev.append(Er)
model.dump()
#model.dump("cohesive elements")
- maxsteps = 2
+ maxsteps = 1
dt = model.getStableTimeStep()*0.1
# choose the timestep
model.setTimeStep(dt)
time = tc
for i in range(0, maxsteps):
print("---- step {0}/{1}".format(i, maxsteps))
time+=dt
functor_r.set_time(time)
# fix displacements of the right boundary
model.applyBC(functor_r, 'right')
d_previous = model.getMaterial("cohesive").getInternalReal("czm_damage")(aka._segment_2).copy()
#print("lda = ",lda)
model.solveStep('explicit_lumped')
model.computeLagrangeMultiplier()
- d[:] = updater.update_d(d_previous,lda)
- #print("d = ",d)
+ #d[:] = updater.update_d(d_previous,lda)
+ model.computeDamage()
+ # print("d = ",d)
model.checkCohesiveInsertion()
U = functor_r.get_imposed_disp()
imposed_disp.append(U)
F = -np.sum(Fint[nodes_right,0])
reaction_force.append(F)
Ep = model.getEnergy("potential")
Ek = model.getEnergy("kinetic")
Ed = model.getEnergy("dissipated")
Er = model.getEnergy("reversible")
E_pot.append(Ep)
E_kin.append(Ek)
E_dis.append(Ed)
E_rev.append(Er)
model.dump()
#model.dump("cohesive elements")
plt.plot(imposed_disp,reaction_force,'-o')
plt.plot([0., w0, wc], [0., sigc*L, 0.], 'k-o')
plt.show()
# -----------------------------------------------------------------------------
# main
# -----------------------------------------------------------------------------
def main():
mesh_file = "triangle.msh"
material_file = "material.dat"
solve(material_file, mesh_file)
# -----------------------------------------------------------------------------
if __name__ == "__main__":
main()
diff --git a/packages/cohesive_element.cmake b/packages/cohesive_element.cmake
index 2f21f8004..f498a5601 100644
--- a/packages/cohesive_element.cmake
+++ b/packages/cohesive_element.cmake
@@ -1,86 +1,90 @@
-#===============================================================================
+ #===============================================================================
# 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 <http://www.gnu.org/licenses/>.
#
#===============================================================================
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/materials/material_bulk_damage.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_bulk_damage.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_bulk_damage_inline_impl.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_clip.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_clip.hh
- model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_clip_inline_impl.hh
+ #model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_clip_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
+
+ model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_damage.cc
+ model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_damage.hh
+
)
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index 603658d7a..e32b27a2b 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -1,139 +1,140 @@
#===============================================================================
# Copyright (©) 2014-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 <http://www.gnu.org/licenses/>.
#
#===============================================================================
if(NOT SKBUILD)
package_get_all_include_directories(
AKANTU_LIBRARY_INCLUDE_DIRS
)
package_get_all_external_informations(
PRIVATE_INCLUDE AKANTU_PRIVATE_EXTERNAL_INCLUDE_DIR
INTERFACE_INCLUDE AKANTU_INTERFACE_EXTERNAL_INCLUDE_DIR
LIBRARIES AKANTU_EXTERNAL_LIBRARIES
)
endif()
set(PYAKANTU_SRCS
py_aka_common.cc
py_aka_error.cc
py_akantu.cc
py_boundary_conditions.cc
py_constitutive_law.cc
py_constitutive_law_selector.cc
py_dof_manager.cc
py_dumpable.cc
py_fe_engine.cc
py_group_manager.cc
py_integration_scheme.cc
py_mesh.cc
py_model.cc
py_parser.cc
py_solver.cc
)
package_is_activated(solid_mechanics _is_activated)
if (_is_activated)
list(APPEND PYAKANTU_SRCS
py_solid_mechanics_model.cc
py_material.cc
)
endif()
package_is_activated(cohesive_element _is_activated)
if (_is_activated)
list(APPEND PYAKANTU_SRCS
py_solid_mechanics_model_cohesive.cc
+ py_solid_mechanics_model_cohesive_damage.cc
py_fragment_manager.cc
)
endif()
package_is_activated(diffusion _is_activated)
if (_is_activated)
list(APPEND PYAKANTU_SRCS
py_heat_transfer_model.cc
)
endif()
package_is_activated(contact_mechanics _is_activated)
if(_is_activated)
list(APPEND PYAKANTU_SRCS
py_contact_mechanics_model.cc
py_model_couplers.cc
)
endif()
package_is_activated(phase_field _is_activated)
if (_is_activated)
list(APPEND PYAKANTU_SRCS
py_phase_field_model.cc
)
endif()
package_is_activated(structural_mechanics _is_activated)
if (_is_activated)
list(APPEND PYAKANTU_SRCS
py_structural_mechanics_model.cc
)
endif()
pybind11_add_module(py11_akantu ${PYAKANTU_SRCS})
# to avoid compilation warnings from pybind11
target_include_directories(py11_akantu
SYSTEM BEFORE
PRIVATE ${PYBIND11_INCLUDE_DIR}
PRIVATE ${pybind11_INCLUDE_DIR}
PRIVATE ${Python_INCLUDE_DIRS})
target_link_libraries(py11_akantu PUBLIC akantu)
set_target_properties(py11_akantu PROPERTIES
DEBUG_POSTFIX ""
LIBRARY_OUTPUT_DIRECTORY akantu)
if (CMAKE_CXX_COMPILER_ID STREQUAL "Clang")
target_compile_options(py11_akantu PUBLIC -fsized-deallocation)
endif()
file(COPY akantu DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
if(NOT Python_MAJOR)
set(Python_VERSION_MAJOR ${PYTHON_VERSION_MAJOR})
set(Python_VERSION_MINOR ${PYTHON_VERSION_MINOR})
endif()
if(NOT SKBUILD)
set(_python_install_dir
${CMAKE_INSTALL_LIBDIR}/python${Python_VERSION_MAJOR}.${Python_VERSION_MINOR}/site-packages/akantu)
else()
set(_python_install_dir python/akantu)
endif()
install(TARGETS py11_akantu
LIBRARY DESTINATION ${_python_install_dir})
if(NOT SKBUILD)
install(DIRECTORY akantu
DESTINATION ${_python_install_dir}
FILES_MATCHING PATTERN "*.py")
endif()
diff --git a/python/py_akantu.cc b/python/py_akantu.cc
index 245d02bd6..0d11240ac 100644
--- a/python/py_akantu.cc
+++ b/python/py_akantu.cc
@@ -1,169 +1,171 @@
/**
* Copyright (©) 2018-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 <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#include "aka_config.hh"
// for NLSNotConvergedException
#include "non_linear_solver.hh"
/* -------------------------------------------------------------------------- */
#include "py_aka_common.hh"
#include "py_aka_error.hh"
#include "py_boundary_conditions.hh"
#include "py_constitutive_law.hh"
#include "py_constitutive_law_selector.hh"
#include "py_dof_manager.hh"
#include "py_dumpable.hh"
#include "py_fe_engine.hh"
#include "py_group_manager.hh"
#include "py_integration_scheme.hh"
#include "py_mesh.hh"
#include "py_model.hh"
#include "py_parser.hh"
#include "py_solver.hh"
#if defined(AKANTU_SOLID_MECHANICS)
#include "py_material.hh"
#include "py_solid_mechanics_model.hh"
#endif
#if defined(AKANTU_DIFFUSION)
#include "py_heat_transfer_model.hh"
#endif
#if defined(AKANTU_COHESIVE_ELEMENT)
#include "py_fragment_manager.hh"
#include "py_solid_mechanics_model_cohesive.hh"
+#include "py_solid_mechanics_model_cohesive_damage.hh"
#endif
#if defined(AKANTU_CONTACT_MECHANICS)
#include "py_contact_mechanics_model.hh"
#include "py_model_couplers.hh"
#endif
#if defined(AKANTU_PHASE_FIELD)
#include "py_phase_field_model.hh"
#endif
#if defined(AKANTU_STRUCTURAL_MECHANICS)
#include "py_structural_mechanics_model.hh"
#endif
/* -------------------------------------------------------------------------- */
#include <aka_error.hh>
/* -------------------------------------------------------------------------- */
#include <pybind11/pybind11.h>
/* -------------------------------------------------------------------------- */
#include <iostream>
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
namespace akantu {
void register_all(pybind11::module & mod) {
register_initialize(mod);
register_enums(mod);
register_error(mod);
register_functions(mod);
register_parser(mod);
register_solvers(mod);
register_dumpable(mod);
register_group_manager(mod);
register_mesh(mod);
register_fe_engine(mod);
register_integration_schemes(mod);
register_dof_manager(mod);
register_boundary_conditions(mod);
register_model(mod);
register_constitutive_law_selector(mod);
register_constitutive_law_internal_handler(mod);
#if defined(AKANTU_DIFFUSION)
register_heat_transfer_model(mod);
#endif
#if defined(AKANTU_SOLID_MECHANICS)
register_solid_mechanics_model(mod);
register_material(mod);
#endif
#if defined(AKANTU_COHESIVE_ELEMENT)
register_solid_mechanics_model_cohesive(mod);
+ register_solid_mechanics_model_cohesive_damage(mod);
register_fragment_manager(mod);
#endif
#if defined(AKANTU_STRUCTURAL_MECHANICS)
register_structural_mechanics_model(mod);
#endif
#if defined(AKANTU_CONTACT_MECHANICS)
register_contact_mechanics_model(mod);
register_model_couplers(mod);
#endif
#if defined(AKANTU_PHASE_FIELD)
register_phase_field_model(mod);
register_phase_field_coupler(mod);
#endif
}
} // namespace akantu
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
PYBIND11_MODULE(py11_akantu, mod) {
mod.doc() = "Akantu python interface";
static py::exception<akantu::debug::Exception> akantu_exception(mod,
"Exception");
static py::exception<akantu::debug::NLSNotConvergedException>
akantu_exception_nls_not_converged(mod, "NLSNotConvergedException");
py::register_exception_translator([](std::exception_ptr ptr) {
try {
if (ptr) {
std::rethrow_exception(ptr);
}
} catch (akantu::debug::NLSNotConvergedException & e) {
akantu_exception_nls_not_converged(e.info().c_str());
} catch (akantu::debug::Exception & e) {
if (akantu::debug::debugger.printBacktrace()) {
akantu::debug::printBacktrace();
}
akantu_exception(e.info().c_str());
}
});
akantu::register_all(mod);
mod.def("has_mpi",
[]() {
#if defined(AKANTU_USE_MPI)
return true;
#else
return false;
#endif
})
.def("getVersion", &akantu::getVersion);
} // Module akantu
diff --git a/python/py_solid_mechanics_model_cohesive.cc b/python/py_solid_mechanics_model_cohesive.cc
index 17ddf99a6..68293f042 100644
--- a/python/py_solid_mechanics_model_cohesive.cc
+++ b/python/py_solid_mechanics_model_cohesive.cc
@@ -1,126 +1,119 @@
/**
* 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 <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#include "py_aka_array.hh"
/* -------------------------------------------------------------------------- */
#include <element_synchronizer.hh>
#include <non_linear_solver.hh>
#include <solid_mechanics_model_cohesive.hh>
/* -------------------------------------------------------------------------- */
#include <pybind11/pybind11.h>
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
/* -------------------------------------------------------------------------- */
namespace akantu {
#define def_function_nocopy(func_name) \
def( \
#func_name, \
[](SolidMechanicsModelCohesive & self) -> decltype(auto) { \
return self.func_name(); \
}, \
py::return_value_policy::reference)
#define def_function(func_name) \
def(#func_name, [](SolidMechanicsModelCohesive & self) -> decltype(auto) { \
return self.func_name(); \
})
void register_solid_mechanics_model_cohesive(py::module & mod) {
py::class_<CohesiveElementInserter>(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_<SolidMechanicsModelCohesiveOptions, SolidMechanicsModelOptions>(
mod, "SolidMechanicsModelCohesiveOptions")
.def(py::init<AnalysisMethod, bool>(),
py::arg("analysis_method") = _explicit_lumped_mass,
py::arg("is_extrinsic") = false);
py::class_<SolidMechanicsModelCohesive, SolidMechanicsModel>(
mod, "SolidMechanicsModelCohesive")
.def(py::init<Mesh &, Int, const ID &>(), 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("checkCohesiveInsertion",
- &SolidMechanicsModelCohesive::checkCohesiveInsertion)
- .def("computeLagrangeMultiplier",
- &SolidMechanicsModelCohesive::computeLagrangeMultiplier)
- .def("computeDamage",
- &SolidMechanicsModelCohesive::computeDamage)
.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/python/py_solid_mechanics_model_cohesive.cc b/python/py_solid_mechanics_model_cohesive_damage.cc
similarity index 54%
copy from python/py_solid_mechanics_model_cohesive.cc
copy to python/py_solid_mechanics_model_cohesive_damage.cc
index 17ddf99a6..8fab3c850 100644
--- a/python/py_solid_mechanics_model_cohesive.cc
+++ b/python/py_solid_mechanics_model_cohesive_damage.cc
@@ -1,126 +1,87 @@
/**
* 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 <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#include "py_aka_array.hh"
/* -------------------------------------------------------------------------- */
#include <element_synchronizer.hh>
#include <non_linear_solver.hh>
-#include <solid_mechanics_model_cohesive.hh>
+#include <solid_mechanics_model_cohesive_damage.hh>
/* -------------------------------------------------------------------------- */
#include <pybind11/pybind11.h>
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
/* -------------------------------------------------------------------------- */
namespace akantu {
#define def_function_nocopy(func_name) \
def( \
#func_name, \
- [](SolidMechanicsModelCohesive & self) -> decltype(auto) { \
+ [](SolidMechanicsModelCohesiveDamage & self) -> decltype(auto) { \
return self.func_name(); \
}, \
py::return_value_policy::reference)
#define def_function(func_name) \
- def(#func_name, [](SolidMechanicsModelCohesive & self) -> decltype(auto) { \
+ def(#func_name, [](SolidMechanicsModelCohesiveDamage & self) -> decltype(auto) { \
return self.func_name(); \
})
-void register_solid_mechanics_model_cohesive(py::module & mod) {
- py::class_<CohesiveElementInserter>(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_<SolidMechanicsModelCohesiveOptions, SolidMechanicsModelOptions>(
- mod, "SolidMechanicsModelCohesiveOptions")
- .def(py::init<AnalysisMethod, bool>(),
- py::arg("analysis_method") = _explicit_lumped_mass,
- py::arg("is_extrinsic") = false);
-
- py::class_<SolidMechanicsModelCohesive, SolidMechanicsModel>(
- mod, "SolidMechanicsModelCohesive")
+void register_solid_mechanics_model_cohesive_damage(py::module & mod) {
+
+ py::class_<SolidMechanicsModelCohesiveDamage, SolidMechanicsModel>(
+ mod, "SolidMechanicsModelCohesiveDamage")
.def(py::init<Mesh &, Int, const ID &>(), 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)
+ &SolidMechanicsModelCohesiveDamage::checkCohesiveStress)
.def("checkCohesiveInsertion",
- &SolidMechanicsModelCohesive::checkCohesiveInsertion)
+ &SolidMechanicsModelCohesiveDamage::checkCohesiveInsertion)
.def("computeLagrangeMultiplier",
- &SolidMechanicsModelCohesive::computeLagrangeMultiplier)
+ &SolidMechanicsModelCohesiveDamage::computeLagrangeMultiplier)
.def("computeDamage",
- &SolidMechanicsModelCohesive::computeDamage)
+ &SolidMechanicsModelCohesiveDamage::computeDamage)
.def("getElementInserter",
- &SolidMechanicsModelCohesive::getElementInserter,
+ &SolidMechanicsModelCohesiveDamage::getElementInserter,
py::return_value_policy::reference)
- .def("getStressOnFacets", &SolidMechanicsModelCohesive::getStressOnFacets,
+ .def("getStressOnFacets", &SolidMechanicsModelCohesiveDamage::getStressOnFacets,
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
- .def("getTangents", &SolidMechanicsModelCohesive::getTangents,
+ .def("getTangents", &SolidMechanicsModelCohesiveDamage::getTangents,
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
.def_function_nocopy(getLambda)
.def("updateAutomaticInsertion",
- &SolidMechanicsModelCohesive::updateAutomaticInsertion);
+ &SolidMechanicsModelCohesiveDamage::updateAutomaticInsertion);
}
} // namespace akantu
diff --git a/python/py_solid_mechanics_model_cohesive_damage.hh b/python/py_solid_mechanics_model_cohesive_damage.hh
new file mode 100644
index 000000000..83436997f
--- /dev/null
+++ b/python/py_solid_mechanics_model_cohesive_damage.hh
@@ -0,0 +1,33 @@
+/**
+ * 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 <http://www.gnu.org/licenses/>.
+ */
+
+/* -------------------------------------------------------------------------- */
+#include <pybind11/pybind11.h>
+
+#ifndef AKANTU_PY_SOLID_MECHANICS_MODEL_COHESIVE_DAMAGE_HH_
+#define AKANTU_PY_SOLID_MECHANICS_MODEL_COHESIVE_DAMAGE_HH_
+
+namespace akantu {
+
+void register_solid_mechanics_model_cohesive_damage(pybind11::module & mod);
+
+} // namespace akantu
+
+#endif // AKANTU_PY_SOLID_MECHANICS_MODEL_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
index 6b18b7e9d..6097e6bbf 100644
--- 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
@@ -1,332 +1,333 @@
/**
* 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 <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#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 "solid_mechanics_model_cohesive_damage.hh"
#include "sparse_matrix.hh"
/* -------------------------------------------------------------------------- */
#include <algorithm>
#include <numeric>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
template <Int dim>
MaterialCohesiveDamageExtrinsic<dim>::MaterialCohesiveDamageExtrinsic(SolidMechanicsModel & model,
const ID & id)
: MaterialCohesiveDamage(model, id),
lambda(registerInternal<Real, FacetInternalField>("lambda",
spatial_dimension)),
czm_damage(
registerInternal<Real, FacetInternalField>("czm_damage", 1)),
insertion_stress(registerInternal<Real, CohesiveInternalField>(
"insertion_stress", dim)),
is_new_crack(registerInternal<Real, CohesiveInternalField>(
"is_new_crack", 1))
{
if(not this->model->getIsExtrinsic())
{
AKANTU_EXCEPTION(
"MaterialCohesiveDamageExtrinsic can be used only with extrinsic cohesive elements");
}
// czm_damage.setDefaultValue(0.1);
is_new_crack.setDefaultValue(0.);
}
/* -------------------------------------------------------------------------- */
template <Int dim> // NOLINTNEXTLINE(readability-function-cognitive-complexity)
void MaterialCohesiveDamageExtrinsic<dim>::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 & lambda = this->lambda(type_facet);
auto & damage = this->czm_damage(type_facet);
auto & ins_stress = insertion_stress(type_cohesive);
auto & is_new_crack = this->is_new_crack(type_cohesive);
Int old_nb_quad_points = is_new_crack.size();
Int new_nb_quad_points = 0;
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
auto damage_begin = damage.begin();
auto lambda_begin = make_view<dim>(lambda).begin();
std::vector<Vector<Real>> new_normal_traction;
// 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;
}
// compute the effective norm on each quadrature point of the facet
Real d(0.);
for (Int q = 0; q < nb_quad_facet; ++q) {
auto current_quad = facet * nb_quad_facet + q;
d+=damage_begin[current_quad];
}
d/=nb_quad_facet;
// verify if the effective stress overcomes the threshold
if (not Math::are_float_equal(d, 0.)) {
f_insertion(facet) = true;
if (check_only) {
continue;
}
for (Int q = 0; q < nb_quad_facet; ++q) {
auto current_quad = facet * nb_quad_facet + q;
new_normal_traction.emplace_back(lambda_begin[current_quad]);
new_nb_quad_points++;
}
}
}
// update material data for the new elements
ins_stress.resize(old_nb_quad_points + new_nb_quad_points);
is_new_crack.resize(old_nb_quad_points + new_nb_quad_points);
for (Int q = 0; q < new_nb_quad_points; ++q) {
is_new_crack(old_nb_quad_points + q) = 1;
for (Int d = 0; d < dim; ++d) {
ins_stress(old_nb_quad_points + q, d) = new_normal_traction[q](d);
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialCohesiveDamageExtrinsic<dim>::computeLambdaOnQuad(ElementType type,
GhostType ghost_type) {
+ SolidMechanicsModelCohesiveDamage * model_lambda = &aka::as_type<SolidMechanicsModelCohesiveDamage>(*model);
auto & fem_lambda = model->getFEEngine("LambdaFEEngine");
- const auto & lambda = model->getLambda();
+ const auto & lambda = model_lambda->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 <Int dim>
void MaterialCohesiveDamageExtrinsic<dim>::computeTraction(ElementType el_type,
GhostType ghost_type) {
std::cout << " In MaterialCohesiveDamageExtrinsic<dim>::computeTraction " << std::endl ;
const auto & mesh_facets = model->getMeshFacets();
for (const auto & type : getElementFilter().elementTypes(
spatial_dimension, ghost_type, _ek_cohesive)) {
auto type_facet = Mesh::getFacetType(type);
auto nb_quad_facet =
model->getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet);
#ifndef AKANTU_NDEBUG
auto nb_quad_cohesive = model->getFEEngine("CohesiveFEEngine")
.getNbIntegrationPoints(type);
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
auto & elem_filter = this->getElementFilter(type, ghost_type);
auto nb_element = elem_filter.size();
if (nb_element == 0) {
continue;
}
const auto & element_to_facet =
mesh_facets.getSubelementToElement(el_type,ghost_type);
auto insertion_stress_begin = make_view<dim>(insertion_stress(el_type,ghost_type)).begin();
auto is_new_crack_begin = is_new_crack(el_type,ghost_type).begin();
auto damage_begin = czm_damage(type_facet,ghost_type).begin();
auto opening_begin = make_view<dim>(opening(el_type,ghost_type)).begin();
auto traction_begin = make_view<dim>(tractions(el_type,ghost_type)).begin();
for (auto && elem : elem_filter) {
auto && facet = -1;
auto && facet_1 = element_to_facet(elem,0);
auto && facet_2 = element_to_facet(elem,1);
facet = facet_1<facet_2?facet_1.element:facet_2.element;
for (Int q = 0; q < nb_quad_facet; ++q) {
auto current_quad = elem * nb_quad_facet + q;
auto && is_new = is_new_crack_begin[current_quad];
auto && t = traction_begin[current_quad];
if(is_new > 0)
{
t = insertion_stress_begin[current_quad];
is_new = 0;
}
else
{
auto current_facet_quad = facet * nb_quad_facet + q;
auto && w = opening_begin[current_quad];
auto && d = damage_begin[current_facet_quad];
std::cout << "facet = " << facet << std::endl ;
std::cout << " w = " << w << std::endl ;
std::cout << " d = " << d << std::endl;
t = k*stiffness(d)*w;
std::cout << " t = " << t << std::endl;
}
}
}
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialCohesiveDamageExtrinsic<dim>::computeLambda(GhostType ghost_type)
{
std::cout << " In MaterialCohesiveDamageExtrinsic<dim>::computeLambda " << std::endl ;
const auto & mesh_facets = model->getMeshFacets();
for (const auto & type_facet : mesh_facets.elementTypes(dim - 1)) {
auto type_cohesive = FEEngine::getCohesiveElementType(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
const auto & element_to_subelement =
mesh_facets.getElementToSubelement(type_facet, ghost_type);
const auto & normals = model->getFEEngine("FacetsFEEngine")
.getNormalsOnIntegrationPoints(type_facet);
const auto & f_stress = model->getStressOnFacets(type_facet);
auto facet_stress_begin = make_view(f_stress, dim, dim, 2).begin();
auto normal_begin = make_view<dim>(normals).begin();
auto damage_begin = czm_damage(type_facet,ghost_type).begin();
auto opening_begin = make_view<dim>(opening(type_cohesive,ghost_type)).begin();
for (auto data :
enumerate(make_view(lambda(type_facet), dim,nb_quad_facet))) {
auto && f = std::get<0>(data);
auto && lda = std::get<1>(data);
auto && elem = element_to_subelement(f);
auto && e = -1;
for(auto ec : elem)
{
if(ec!=ElementNull)
{
if(ec.kind() == _ek_cohesive)
{
e = ec.element;
continue;
}
}
}
if(e < 0)
{
for (Int q = 0; q < nb_quad_facet; ++q) {
auto current_quad = f * nb_quad_facet + q;
auto && normal = normal_begin[current_quad];
auto && facet_stress = facet_stress_begin[current_quad];
auto && stress_1 = facet_stress(0);
auto && stress_2 = facet_stress(1);
auto && stress_avg = (stress_1 + stress_2) / 2.;
lda(q) = stress_avg*normal;
}
}
else
{
for (Int q = 0; q < nb_quad_facet; ++q) {
auto current_quad = f * nb_quad_facet + q;
auto current_elem_quad = e * nb_quad_facet + q;
auto && w = opening_begin[current_elem_quad];
auto && d = damage_begin[current_quad];
lda(q) = k*augmented_stiffness(d)*w;
std::cout << "f = " << f << std::endl ;
std::cout << " w = " << w << std::endl ;
std::cout << " d = " << d << std::endl ;
std::cout << " lda = " << lda(q) << std::endl ;
}
}
}
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialCohesiveDamageExtrinsic<dim>::computeDamage(GhostType ghost_type)
{
const auto & mesh_facets = model->getMeshFacets();
for (const auto & type : mesh_facets.elementTypes(dim - 1)) {
auto & elem_filter = getFacetFilter(type);
auto nb_element = elem_filter.size();
if (nb_element == 0) {
continue;
}
for (auto && args : getArguments(type, ghost_type)) {
computeDamageOnQuad(args);
}
}
}
/* -------------------------------------------------------------------------- */
template class MaterialCohesiveDamageExtrinsic<1>;
template class MaterialCohesiveDamageExtrinsic<2>;
template class MaterialCohesiveDamageExtrinsic<3>;
const bool material_is_alocated_cohesive_damage [[maybe_unused]] =
instantiateMaterial<MaterialCohesiveDamageExtrinsic>("cohesive_damage_extrinsic");
} // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.cc
index a1a62f154..7fe1fc28e 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.cc
@@ -1,488 +1,490 @@
/**
* 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 <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#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 "solid_mechanics_model_cohesive_damage.hh"
#include "sparse_matrix.hh"
/* -------------------------------------------------------------------------- */
#include <algorithm>
#include <numeric>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
template <Int dim>
MaterialCohesiveDamageIntrinsic<dim>::MaterialCohesiveDamageIntrinsic(SolidMechanicsModel & model,
const ID & id)
: MaterialCohesiveDamage(model, id),
lambda(registerInternal<Real, CohesiveInternalField>("lambda",
spatial_dimension)),
err_openings(registerInternal<Real, CohesiveInternalField>(
"err_openings", spatial_dimension)),
czm_damage(
registerInternal<Real, CohesiveInternalField>("czm_damage", 1)) {
if(this->model->getIsExtrinsic())
{
AKANTU_EXCEPTION(
"MaterialCohesiveDamageIntrinsic can be used only with intrinsic cohesive elements");
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialCohesiveDamageIntrinsic<dim>::assembleInternalForces(GhostType ghost_type) {
AKANTU_DEBUG_IN();
auto & internal_force = const_cast<Array<Real> &>(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<Array<Real>>(
nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes);
auto err_opening_cpy = std::make_shared<Array<Real>>(
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<Array<Real>>(
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<Array<Real>>(
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<Array<Real>>(
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 <Int dim>
void MaterialCohesiveDamageIntrinsic<dim>::assembleStiffnessMatrix(
GhostType ghost_type) {
AKANTU_DEBUG_IN();
// auto & lambda_connectivities =
// model->getMesh().getElementalData<Idx>("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<Array<Real>>(
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<Real> 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<Array<Real>>(
nb_element * nb_quadrature_points,
spatial_dimension * spatial_dimension, "tangent_stiffness_matrix_uu");
auto tangent_stiffness_matrix_ll = std::make_unique<Array<Real>>(
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<Array<Real>>(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<Array<Real>>(
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<Array<Real>>(
nb_element * nb_quadrature_points, size_at_nt_dul_n, "A^t*N^t*Dul*N");
Matrix<Real> 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<Array<Real>>(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<Array<Real>>(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<Array<Real>>(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();
+ SolidMechanicsModelCohesiveDamage * model_lambda = &aka::as_type<SolidMechanicsModelCohesiveDamage>(*model);
auto lambda_connectivity =
- model->getLambdaMesh().getConnectivity(underlying_type, ghost_type);
+ model_lambda->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 <Int dim>
void MaterialCohesiveDamageIntrinsic<dim>::computeDamage(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(type, ghost_type)) {
computeDamageOnQuad(args);
}
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialCohesiveDamageIntrinsic<dim>::computeLambdaOnQuad(ElementType type,
GhostType ghost_type) {
+ SolidMechanicsModelCohesiveDamage * model_lambda = &aka::as_type<SolidMechanicsModelCohesiveDamage>(*model);
auto & fem_lambda = model->getFEEngine("LambdaFEEngine");
- const auto & lambda = model->getLambda();
+ const auto & lambda = model_lambda->getLambda();
auto & lambda_on_quad = this->lambda(type, ghost_type);
// std::cout << "lambda in : MaterialCohesiveDamageIntrinsic<dim>::computeLambdaOnQuad " << std::endl;
// lambda.printself(std::cout << std::endl);
// ArrayPrintHelper<true>::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 : MaterialCohesiveDamageIntrinsic<dim>::computeLambdaOnQuad " << std::endl;
// lambda_on_quad.printself(std::cout << std::endl);
// ArrayPrintHelper<true>::print_content(lambda_on_quad,std::cout,0);
// std::cout << "shapes_lambda in : MaterialCohesiveDamageIntrinsic<dim>::computeLambdaOnQuad " << std::endl;
// auto shapes_lambda = fem_lambda.getShapes(underlying_type,ghost_type,0);
// ArrayPrintHelper<true>::print_content(shapes_lambda,std::cout,0);
// std::cout << "shapes_cohesive in : MaterialCohesiveDamageIntrinsic<dim>::computeLambdaOnQuad " << std::endl;
// auto & fem_cohesive =
// model->getFEEngineClass<MyFEEngineCohesiveType>("CohesiveFEEngine");
// auto shapes_cohesive = fem_cohesive.getShapes(type,ghost_type,0);
// ArrayPrintHelper<true>::print_content(shapes_cohesive,std::cout,0);
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialCohesiveDamageIntrinsic<dim>::computeTraction(ElementType el_type,
GhostType ghost_type) {
auto & elem_filter = getElementFilter(el_type, ghost_type);
auto nb_element = elem_filter.size();
if (nb_element > 0) {
computeLambdaOnQuad(el_type, ghost_type);
// auto & traction = tractions(type, ghost_type);
// std::cout << "traction in : MaterialCohesiveDamageIntrinsic<dim>::computeTraction " << std::endl;
// ArrayPrintHelper<true>::print_content(traction,std::cout,0);
// auto & lambda_ = lambda(type, ghost_type);
// std::cout << "lambda in : MaterialCohesiveDamageIntrinsic<dim>::computeTraction " << std::endl;
// ArrayPrintHelper<true>::print_content(lambda_,std::cout,0);
// lambda.printself(std::cout);
for (auto && args : getArguments(el_type, ghost_type)) {
computeTractionOnQuad(args);
}
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialCohesiveDamageIntrinsic<dim>::computeTangentTraction(
ElementType el_type, Array<Real> & tangent_matrix_uu,
Array<Real> & tangent_matrix_ll, GhostType ghost_type) {
AKANTU_DEBUG_IN();
for (auto && [args, tangent_uu, tangent_ll] :
zip(getArguments(el_type, ghost_type),
make_view<dim, dim>(tangent_matrix_uu),
make_view<dim, dim>(tangent_matrix_ll))) {
computeTangentTractionOnQuad(tangent_uu, tangent_ll, args);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template class MaterialCohesiveDamageIntrinsic<1>;
template class MaterialCohesiveDamageIntrinsic<2>;
template class MaterialCohesiveDamageIntrinsic<3>;
const bool material_is_alocated_cohesive_damage [[maybe_unused]] =
instantiateMaterial<MaterialCohesiveDamageIntrinsic>("cohesive_damage_intrinsic");
} // namespace akantu
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 b45a736ba..a2f304e7d 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,935 +1,623 @@
/**
* 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 <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#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_damage.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 <algorithm>
/* -------------------------------------------------------------------------- */
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<UInt, UInt>
updateData(NewNodesEvent & nodes_event,
NewElementsEvent & elements_event) override {
auto * cohesive_nodes_event =
dynamic_cast<CohesiveNewNodesEvent *>(&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<DOFManager> & 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<MyFEEngineCohesiveType>("CohesiveFEEngine", mesh,
Model::spatial_dimension);
auto && tmp_material_selector =
std::make_shared<DefaultMaterialCohesiveSelector>(*this);
tmp_material_selector->setFallback(this->getConstitutiveLawSelector());
this->setConstitutiveLawSelector(tmp_material_selector);
this->mesh.registerDumper<DumperParaview>("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<ElementSynchronizer>(
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<CohesiveElementInserter>(
this->mesh, id + ":cohesive_element_inserter");
registerFEEngineObject<MyFEEngineFacetType>(
"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<SolidMechanicsModelCohesiveOptions>(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<FacetSynchronizer>(mesh_facets.getElementSynchronizer());
// synchronizeGhostFacetsConnectivity();
/// create the facet synchronizer for extrinsic simulations
if (is_extrinsic) {
facet_stress_synchronizer = std::make_unique<ElementSynchronizer>(
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<CohesiveMeshGlobalDataUpdater>(*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<MaterialCohesive>(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<Idx>("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<MaterialCohesive>(material)) {
- need_lambda |= aka::as_type<MaterialCohesive>(material).needLambda();
- }
- });
-
- if (need_lambda) {
- lambda =
- std::make_unique<Array<Real>>(0, spatial_dimension, "cohesive lambda");
- previous_lambda =
- std::make_unique<Array<Real>>(0, spatial_dimension, "previous lambda");
- lambda_increment =
- std::make_unique<Array<Real>>(0, spatial_dimension, "lambda increment");
- lambda_blocked_dofs = std::make_unique<Array<bool>>(0, spatial_dimension,
- "lambda_blocked_dofs");
-
- lambda_mesh = std::make_unique<Mesh>(spatial_dimension,
- mesh.getCommunicator(), "lambda_mesh");
- registerFEEngineObject<MyFEEngineLambdaType>("LambdaFEEngine", *lambda_mesh,
- Model::spatial_dimension - 1);
-
- auto & nodes_to_lambda = mesh.getNodalData<Idx>("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<Real> & position = mesh.getNodes();
ElementTypeMapArray<Real> 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<Real> 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<Real>::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<MaterialCohesive>(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<MaterialCohesive>(material);
mat.computeTraction(_not_ghost);
} catch (std::bad_cast & bce) {
}
});
- // ArrayPrintHelper<true>::ArrayPrintHelper::print_content<bool>(*(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<Real> & normals =
this->getFEEngine("FacetsFEEngine")
.getNormalsOnIntegrationPoints(facet_type);
Array<Real> & tangents = this->tangents(facet_type);
Math::compute_tangents(normals, tangents);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::interpolateStress() {
ElementTypeMapArray<Real> by_elem_result("temporary_stress_by_facets", id);
this->for_each_constitutive_law([&](auto && material) {
if (not aka::is_of_type<MaterialCohesive>(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<MaterialCohesive>(material)) {
/// check which not ghost cohesive elements are to be created
auto & mat_cohesive = aka::as_type<MaterialCohesive>(material);
mat_cohesive.checkInsertion();
}
});
/// insert cohesive elements
auto nb_new_elements = inserter->insertElements();
AKANTU_DEBUG_OUT();
return nb_new_elements;
}
-/* -------------------------------------------------------------------------- */
-UInt SolidMechanicsModelCohesive::checkCohesiveInsertion() {
- AKANTU_DEBUG_IN();
-
- if (not is_extrinsic) {
- AKANTU_EXCEPTION(
- "This function can only be used for extrinsic cohesive elements");
- }
-
- this->for_each_constitutive_law([&](auto && material) {
- if (aka::is_of_type<MaterialCohesive>(material)) {
- /// check which not ghost cohesive elements are to be created
- auto & mat_cohesive = aka::as_type<MaterialCohesive>(material);
- mat_cohesive.checkInsertion();
- }
- });
-
- /// insert cohesive elements
- auto nb_new_elements = inserter->insertElements();
-
- AKANTU_DEBUG_OUT();
-
- return nb_new_elements;
-}
-
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModelCohesive::computeLagrangeMultiplier() {
- 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<MaterialCohesive>(material)) {
- /// check which not ghost cohesive elements are to be created
- auto & mat_cohesive = aka::as_type<MaterialCohesive>(material);
- mat_cohesive.computeLambda();
- }
- });
-
- AKANTU_DEBUG_OUT();
-}
-
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModelCohesive::computeDamage() {
- AKANTU_DEBUG_IN();
-
- this->for_each_constitutive_law([&](auto && material) {
- if (aka::is_of_type<MaterialCohesiveDamage>(material)) {
- /// check which not ghost cohesive elements are to be created
- auto & mat_cohesive = aka::as_type<MaterialCohesiveDamage>(material);
- mat_cohesive.computeDamage();
- }
- });
-
- AKANTU_DEBUG_OUT();
-}
-
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModelCohesive::updateLambdaMesh() {
-// std::cout << " SolidMechanicsModelCohesive::updateLambdaMesh " << std::endl;
-
- const auto & connectivities = mesh.getConnectivities();
- const auto & initial_nodes = mesh.getNodalData<Idx>("initial_nodes_match");
- auto & nodes_to_lambda = mesh.getNodalData<Idx>("nodes_to_lambda");
- nodes_to_lambda.resize(mesh.getNbNodes(), -1);
-
- auto & lambda_connectivities = MeshAccessor(*lambda_mesh).getConnectivities();
-
- auto lambda_id = lambda->size();
-
- std::vector<Idx> new_nodes;
-
- NewNodesEvent node_event(*lambda_mesh, AKANTU_CURRENT_FUNCTION);
- auto & node_list = node_event.getList();
-// std::cout << "node_list before = " << std::endl ;
-// ArrayPrintHelper<true>::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<true>::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<true>::print_content(connectivity,std::cout,0);
-
-// std::cout << "lambda_connectivity = " << std::endl ;
-// ArrayPrintHelper<true>::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);
-// 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<true>::print_content(node_list,std::cout,0);
-
-// std::cout << "element_list after = " << std::endl ;
-// ArrayPrintHelper<true>::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> & 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<Idx> & new_nodes,
const NewNodesEvent & event) {
- SolidMechanicsModel::onNodesAdded(new_nodes, event);
-
- auto & initial_nodes = mesh.getNodalData<Idx>("initial_nodes_match");
+ AKANTU_DEBUG_IN();
- auto old_max_nodes = initial_nodes.size();
- initial_nodes.resize(mesh.getNbNodes(), -1);
+ SolidMechanicsModel::onNodesAdded(new_nodes, event);
const auto * cohesive_event =
dynamic_cast<const CohesiveNewNodesEvent *>(&event);
if (cohesive_event == nullptr) {
- for (auto && node : new_nodes) {
- initial_nodes(node) = node;
- }
return;
}
- auto old_nodes = cohesive_event->getOldNodesList();
+ const 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();
+ auto copy = [this, &new_nodes, &old_nodes](auto & arr) {
+ auto it = make_view(arr, spatial_dimension).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(*displacement);
+ copy(*blocked_dofs);
+
+ if (velocity) {
+ copy(*velocity);
}
- copy(*blocked_dofs);
- copy(getDOFManager().getSolution("displacement"));
- copy(initial_nodes);
+ if (acceleration) {
+ copy(*acceleration);
+ }
+
+ if (current_position) {
+ copy(*current_position);
+ }
+
+ if (previous_displacement) {
+ copy(*previous_displacement);
+ }
- // correct connectivities
- if (lambda) {
- lambda_blocked_dofs->resize(mesh.getNbNodes(), false);
- copy(*lambda_blocked_dofs);
- updateLambdaMesh();
+ if (displacement_increment) {
+ copy(*displacement_increment);
}
+
+ copy(getDOFManager().getSolution("displacement"));
+
+ AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
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) {
- 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 fadac8a49..59bb582da 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,357 +1,292 @@
/**
* 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 <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#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 <ElementType type, ElementType cohesive_type =
CohesiveFacetProperty<type>::cohesive_type>
struct _helper {
static constexpr int get() {
// return 1;
return ElementClassProperty<cohesive_type>::polynomial_degree;
}
};
template <ElementType type> struct _helper<type, _not_defined> {
static constexpr int get() {
// return 1;
return ElementClassProperty<type>::polynomial_degree; }
};
template <ElementType type> static inline constexpr int getOrder() {
return _helper<type>::get();
}
};
-/* -------------------------------------------------------------------------- */
-struct CohesiveIntegrationOrderFunctor {
- template <ElementType type, ElementType cohesive_type =
- CohesiveFacetProperty<type>::cohesive_type>
- struct _helper {
- static constexpr int get() {
- // return 1;
- return ElementClassProperty<cohesive_type>::polynomial_degree; }
- };
-
- template <ElementType type> struct _helper<type, _not_defined> {
- static constexpr int get() {
- // return 1;
- return ElementClassProperty<type>::polynomial_degree; }
- };
-
- template <ElementType type> static inline constexpr int getOrder() {
- return _helper<type>::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<UInt> &);
AKANTU_GET_MACRO(OldNodesList, old_nodes, const Array<UInt> &);
protected:
Array<UInt> old_nodes;
};
using MyFEEngineCohesiveType =
FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_cohesive>;
using MyFEEngineFacetType =
FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_regular,
FacetsCohesiveIntegrationOrderFunctor>;
- using MyFEEngineLambdaType =
- FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_regular,
- FacetsCohesiveIntegrationOrderFunctor>;
SolidMechanicsModelCohesive(
Mesh & mesh, Int dim = _all_dimensions,
const ID & id = "solid_mechanics_model_cohesive",
const std::shared_ptr<DOFManager> & 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();
- /// function to perform an insertion check on each facet and insert
- /// cohesive elements if needed (returns the number of new cohesive
- /// elements)
- UInt checkCohesiveInsertion();
-
/// interpolate stress on facets
void interpolateStress();
/// update automatic insertion after a change in the element inserter
void updateAutomaticInsertion();
/// insert intrinsic cohesive elements
void insertIntrinsicElements();
- /// compute Lagrange multiplier
- void computeLagrangeMultiplier();
-
- /// compute damage variable
- void computeDamage();
-
// template <SolveConvergenceMethod cmethod, SolveConvergenceCriteria
// criteria> 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:
+protected:
/// 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<Idx> & new_nodes,
const NewNodesEvent & event) override;
void onElementsAdded(const Array<Element> & 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 synchronizeGhostFacetsConnectivity();
+
void updateCohesiveSynchronizers(NewElementsEvent & elements_event);
void updateFacetStressSynchronizer();
friend class CohesiveElementInserter;
/* ------------------------------------------------------------------------ */
/* Data Accessor inherited members */
/* ------------------------------------------------------------------------ */
public:
[[nodiscard]] Int getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const override;
void packData(CommunicationBuffer & buffer, const Array<Element> & elements,
const SynchronizationTag & tag) const override;
void unpackData(CommunicationBuffer & buffer, const Array<Element> & elements,
const SynchronizationTag & tag) override;
protected:
[[nodiscard]] Int
getNbQuadsForFacetCheck(const Array<Element> & elements) const;
template <typename T>
void packFacetStressDataHelper(const ElementTypeMapArray<T> & data_to_pack,
CommunicationBuffer & buffer,
const Array<Element> & elements) const;
template <typename T>
void unpackFacetStressDataHelper(ElementTypeMapArray<T> & data_to_unpack,
CommunicationBuffer & buffer,
const Array<Element> & elements) const;
template <typename T, bool pack_helper>
void packUnpackFacetStressDataHelper(
std::conditional_t<pack_helper, const ElementTypeMapArray<T>,
ElementTypeMapArray<T>> & data_to_pack,
CommunicationBuffer & buffer, const Array<Element> & 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::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:
+protected:
friend class CohesiveMeshGlobalDataUpdater;
- /// lambda array
- std::unique_ptr<Array<Real>> lambda;
-
- /// lambda array at the previous time step
- std::unique_ptr<Array<Real>> previous_lambda;
-
- /// increment of lambda
- std::unique_ptr<Array<Real>> lambda_increment;
-
- /// array specifing if a lambda degree of freedom is blocked or not
- std::unique_ptr<Array<bool>> lambda_blocked_dofs;
-
/// @todo store tangents when normals are computed:
ElementTypeMapArray<Real> tangents;
/// stress on facets on the two sides by quadrature point
ElementTypeMapArray<Real> facet_stress;
/// material to use if a cohesive element is created on a facet
ElementTypeMapArray<Idx> facet_material;
bool is_extrinsic{false};
/// cohesive element inserter
std::unique_ptr<CohesiveElementInserter> inserter;
/// facet stress synchronizer
std::unique_ptr<ElementSynchronizer> facet_stress_synchronizer;
/// cohesive elements synchronizer
std::unique_ptr<ElementSynchronizer> cohesive_synchronizer;
-
- std::unique_ptr<Mesh> lambda_mesh;
};
} // namespace akantu
#include "solid_mechanics_model_cohesive_inline_impl.hh"
#endif /* AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH_ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_damage.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_damage.cc
new file mode 100644
index 000000000..dc465e0aa
--- /dev/null
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_damage.cc
@@ -0,0 +1,505 @@
+/**
+ * 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 <http://www.gnu.org/licenses/>.
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "solid_mechanics_model_cohesive_damage.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_damage.hh"
+#include "mesh_accessor.hh"
+#include "mesh_global_data_updater.hh"
+#include "parser.hh"
+#include "shape_cohesive.hh"
+/* -------------------------------------------------------------------------- */
+#include "dumper_iohelper_paraview.hh"
+/* -------------------------------------------------------------------------- */
+#include <algorithm>
+/* -------------------------------------------------------------------------- */
+
+namespace akantu {
+
+/* -------------------------------------------------------------------------- */
+SolidMechanicsModelCohesiveDamage::SolidMechanicsModelCohesiveDamage(
+ Mesh & mesh, Int dim, const ID & id,
+ const std::shared_ptr<DOFManager> & dof_manager)
+ : SolidMechanicsModelCohesive(mesh, dim, id, dof_manager)
+{}
+
+/* -------------------------------------------------------------------------- */
+void SolidMechanicsModelCohesiveDamage::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<MaterialCohesive>(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<Idx>("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<MaterialCohesive>(material)) {
+ need_lambda |= aka::as_type<MaterialCohesive>(material).needLambda();
+ }
+ });
+
+ if (need_lambda) {
+ lambda =
+ std::make_unique<Array<Real>>(0, spatial_dimension, "cohesive lambda");
+ previous_lambda =
+ std::make_unique<Array<Real>>(0, spatial_dimension, "previous lambda");
+ lambda_increment =
+ std::make_unique<Array<Real>>(0, spatial_dimension, "lambda increment");
+ lambda_blocked_dofs = std::make_unique<Array<bool>>(0, spatial_dimension,
+ "lambda_blocked_dofs");
+
+ lambda_mesh = std::make_unique<Mesh>(spatial_dimension,
+ mesh.getCommunicator(), "lambda_mesh");
+ registerFEEngineObject<MyFEEngineLambdaType>("LambdaFEEngine", *lambda_mesh,
+ Model::spatial_dimension - 1);
+
+ auto & nodes_to_lambda = mesh.getNodalData<Idx>("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
+
+/* -------------------------------------------------------------------------- */
+UInt SolidMechanicsModelCohesiveDamage::checkCohesiveInsertion() {
+ AKANTU_DEBUG_IN();
+
+ if (not is_extrinsic) {
+ AKANTU_EXCEPTION(
+ "This function can only be used for extrinsic cohesive elements");
+ }
+
+ this->for_each_constitutive_law([&](auto && material) {
+ if (aka::is_of_type<MaterialCohesive>(material)) {
+ /// check which not ghost cohesive elements are to be created
+ auto & mat_cohesive = aka::as_type<MaterialCohesive>(material);
+ mat_cohesive.checkInsertion();
+ }
+ });
+
+ /// insert cohesive elements
+ auto nb_new_elements = inserter->insertElements();
+
+ AKANTU_DEBUG_OUT();
+
+ return nb_new_elements;
+}
+
+/* -------------------------------------------------------------------------- */
+void SolidMechanicsModelCohesiveDamage::computeLagrangeMultiplier() {
+ 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<MaterialCohesive>(material)) {
+ /// check which not ghost cohesive elements are to be created
+ auto & mat_cohesive = aka::as_type<MaterialCohesive>(material);
+ mat_cohesive.computeLambda();
+ }
+ });
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+void SolidMechanicsModelCohesiveDamage::computeDamage() {
+ AKANTU_DEBUG_IN();
+
+ this->for_each_constitutive_law([&](auto && material) {
+ if (aka::is_of_type<MaterialCohesiveDamage>(material)) {
+ /// check which not ghost cohesive elements are to be created
+ auto & mat_cohesive = aka::as_type<MaterialCohesiveDamage>(material);
+ mat_cohesive.computeDamage();
+ }
+ });
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+void SolidMechanicsModelCohesiveDamage::updateLambdaMesh() {
+// std::cout << " SolidMechanicsModelCohesiveDamage::updateLambdaMesh " << std::endl;
+
+ const auto & connectivities = mesh.getConnectivities();
+ const auto & initial_nodes = mesh.getNodalData<Idx>("initial_nodes_match");
+ auto & nodes_to_lambda = mesh.getNodalData<Idx>("nodes_to_lambda");
+ nodes_to_lambda.resize(mesh.getNbNodes(), -1);
+
+ auto & lambda_connectivities = MeshAccessor(*lambda_mesh).getConnectivities();
+
+ auto lambda_id = lambda->size();
+
+ std::vector<Idx> new_nodes;
+
+ NewNodesEvent node_event(*lambda_mesh, AKANTU_CURRENT_FUNCTION);
+ auto & node_list = node_event.getList();
+// std::cout << "node_list before = " << std::endl ;
+// ArrayPrintHelper<true>::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<true>::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<true>::print_content(connectivity,std::cout,0);
+
+// std::cout << "lambda_connectivity = " << std::endl ;
+// ArrayPrintHelper<true>::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);
+// 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<true>::print_content(node_list,std::cout,0);
+
+// std::cout << "element_list after = " << std::endl ;
+// ArrayPrintHelper<true>::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 SolidMechanicsModelCohesiveDamage::onElementsAdded(
+ const Array<Element> & element_list, const NewElementsEvent & event) {
+ AKANTU_DEBUG_IN();
+
+// std::cout << " SolidMechanicsModelCohesiveDamage::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 SolidMechanicsModelCohesiveDamage::onNodesAdded(const Array<Idx> & new_nodes,
+ const NewNodesEvent & event) {
+ SolidMechanicsModel::onNodesAdded(new_nodes, event);
+
+ auto & initial_nodes = mesh.getNodalData<Idx>("initial_nodes_match");
+
+ auto old_max_nodes = initial_nodes.size();
+ initial_nodes.resize(mesh.getNbNodes(), -1);
+
+ const auto * cohesive_event =
+ dynamic_cast<const CohesiveNewNodesEvent *>(&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(*displacement);
+ copy(*blocked_dofs);
+
+ if (velocity) {
+ copy(*velocity);
+ }
+
+ if (acceleration) {
+ copy(*acceleration);
+ }
+
+ if (current_position) {
+ copy(*current_position);
+ }
+
+ if (previous_displacement) {
+ copy(*previous_displacement);
+ }
+
+ if (displacement_increment) {
+ copy(*displacement_increment);
+ }
+
+ copy(getDOFManager().getSolution("displacement"));
+
+ // correct connectivities
+ if (lambda) {
+ lambda_blocked_dofs->resize(mesh.getNbNodes(), false);
+ copy(*lambda_blocked_dofs);
+ updateLambdaMesh();
+ }
+}
+
+/* -------------------------------------------------------------------------- */
+ModelSolverOptions SolidMechanicsModelCohesiveDamage::getDefaultSolverOptions(
+ const TimeStepSolverType & type) const {
+ ModelSolverOptions options =
+ SolidMechanicsModel::getDefaultSolverOptions(type);
+
+ 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 SolidMechanicsModelCohesiveDamage::printself(std::ostream & stream,
+ int indent) const {
+ std::string space(indent, AKANTU_INDENT);
+
+ stream << space << "SolidMechanicsModelCohesiveDamage ["
+ << "\n";
+ SolidMechanicsModelCohesive::printself(stream, indent + 2);
+ stream << space << "]\n";
+}
+
+/* -------------------------------------------------------------------------- */
+
+} // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_damage.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_damage.hh
new file mode 100644
index 000000000..0737a6000
--- /dev/null
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_damage.hh
@@ -0,0 +1,157 @@
+/**
+ * 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 <http://www.gnu.org/licenses/>.
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "solid_mechanics_model_cohesive.hh"
+/* -------------------------------------------------------------------------- */
+
+#ifndef AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_DAMAGE_HH_
+#define AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_DAMAGE_HH_
+
+/* -------------------------------------------------------------------------- */
+namespace akantu {
+class FacetSynchronizer;
+class FacetStressSynchronizer;
+class ElementSynchronizer;
+} // namespace akantu
+
+namespace akantu {
+
+/* -------------------------------------------------------------------------- */
+struct CohesiveIntegrationOrderFunctor {
+ template <ElementType type, ElementType cohesive_type =
+ CohesiveFacetProperty<type>::cohesive_type>
+ struct _helper {
+ static constexpr int get() {
+ // return 1;
+ return ElementClassProperty<cohesive_type>::polynomial_degree; }
+ };
+
+ template <ElementType type> struct _helper<type, _not_defined> {
+ static constexpr int get() {
+ // return 1;
+ return ElementClassProperty<type>::polynomial_degree; }
+ };
+
+ template <ElementType type> static inline constexpr int getOrder() {
+ return _helper<type>::get();
+ }
+};
+
+/* -------------------------------------------------------------------------- */
+/* Solid Mechanics Model for Cohesive damage elements */
+/* -------------------------------------------------------------------------- */
+class SolidMechanicsModelCohesiveDamage : public SolidMechanicsModelCohesive {
+
+ using MyFEEngineLambdaType =
+ FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_regular,
+ FacetsCohesiveIntegrationOrderFunctor>;
+
+ /* ------------------------------------------------------------------------ */
+ /* Constructors/Destructors */
+ /* ------------------------------------------------------------------------ */
+public:
+ SolidMechanicsModelCohesiveDamage(
+ Mesh & mesh, Int dim = _all_dimensions,
+ const ID & id = "solid_mechanics_model_cohesive_damage",
+ const std::shared_ptr<DOFManager> & dof_manager = nullptr);
+
+ /* ------------------------------------------------------------------------ */
+ /* Methods */
+ /* ------------------------------------------------------------------------ */
+
+public:
+
+ /// function to perform an insertion check on each facet and insert
+ /// cohesive elements if needed (returns the number of new cohesive
+ /// elements)
+ UInt checkCohesiveInsertion();
+
+ /// compute Lagrange multiplier
+ void computeLagrangeMultiplier();
+
+ /// compute damage variable
+ void computeDamage();
+
+protected:
+
+ /// initialize cohesive material
+ void initConstitutiveLaws() override;
+
+ /// function to print the contain of the class
+ void printself(std::ostream & stream, int indent = 0) const override;
+
+ /* ------------------------------------------------------------------------ */
+ /* SolidMechanicsModelEventHandler inherited members */
+ /* ------------------------------------------------------------------------ */
+public:
+
+ [[nodiscard]] ModelSolverOptions
+ getDefaultSolverOptions(const TimeStepSolverType & type) const override;
+
+ void updateLambdaMesh();
+
+ /* ------------------------------------------------------------------------ */
+ /* Accessors */
+ /* ------------------------------------------------------------------------ */
+public:
+
+ /// 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;
+ }
+
+protected:
+ void onNodesAdded(const Array<Idx> & new_nodes,
+ const NewNodesEvent & event) override;
+ void onElementsAdded(const Array<Element> & element_list,
+ const NewElementsEvent & event) override;
+
+ /* ------------------------------------------------------------------------ */
+ /* Class Members */
+ /* ------------------------------------------------------------------------ */
+private:
+
+ /// lambda array
+ std::unique_ptr<Array<Real>> lambda;
+
+ /// lambda array at the previous time step
+ std::unique_ptr<Array<Real>> previous_lambda;
+
+ /// increment of lambda
+ std::unique_ptr<Array<Real>> lambda_increment;
+
+ /// array specifing if a lambda degree of freedom is blocked or not
+ std::unique_ptr<Array<bool>> lambda_blocked_dofs;
+
+ std::unique_ptr<Mesh> lambda_mesh;
+};
+
+} // namespace akantu
+
+#endif /* AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_DAMAGE_HH_ */