diff --git a/examples/c++/solid_mechanics_cohesive_model/my_first_test/material.dat b/examples/c++/solid_mechanics_cohesive_model/my_first_test/material.dat
index 1f6fbb030..2f627470c 100644
--- a/examples/c++/solid_mechanics_cohesive_model/my_first_test/material.dat
+++ b/examples/c++/solid_mechanics_cohesive_model/my_first_test/material.dat
@@ -1,13 +1,13 @@
material elastic [
name = steel
rho = 1 # density
- E = 1 # young's modulus
- nu = 0.1 # poisson's ratio
+ E = 3e10 # young's modulus
+ nu = 0. # poisson's ratio
]
-material cohesive_damage [
+material cohesive_damage_intrinsic [
name = cohesive
- sigma_c = 1.0
- k = 1.0
- G_c = 1.0
+ sigma_c = 3e6
+ k = 75e9
+ G_c = 120
]
diff --git a/examples/c++/solid_mechanics_cohesive_model/my_first_test/my_first_test.py b/examples/c++/solid_mechanics_cohesive_model/my_first_test/my_first_test.py
new file mode 100644
index 000000000..eda628221
--- /dev/null
+++ b/examples/c++/solid_mechanics_cohesive_model/my_first_test/my_first_test.py
@@ -0,0 +1,244 @@
+#!/usr/bin/env python3
+__copyright__ = (
+ "Copyright (©) 2019-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)"
+ "Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)"
+)
+__license__ = "LGPLv3"
+
+
+import sys
+sys.path.insert(0,"/home/ble/buildAkantu/buildDebug/python")
+print(sys.path)
+
+import akantu as aka
+import numpy as np
+
+import scipy.optimize
+
+import matplotlib.pyplot as plt
+
+class degradation_function:
+ # g(d)
+ def val(self,d):
+ return
+ # dg_dd(d)
+ def deriv(self,d):
+ return
+ # 1+g(d)
+ def augmented_stiffness(self,d):
+ return
+ def augmented_compliance(self,d):
+ # 1/(1+g(d))
+ return
+ def augmented_compliance_deriv(self,d):
+ # d()1/(1+g(d))_dd
+ return
+
+class degradation_function_linear_czm(degradation_function):
+ def val(self,d):
+ return 1./d-1.
+ def deriv(self,d):
+ return -1./d**2
+ def augmented_stiffness(self,d):
+ return 1./d
+ def augmented_compliance(self,d):
+ return d
+ def augmented_compliance_deriv(self,d):
+ return 1.
+
+class softening_function:
+ def val(self,d):
+ return
+ def deriv(self,d):
+ return
+
+class softening_function_linear_czm:
+ def __init__(self,k,tc,Gc):
+ wc = 2.*Gc/tc
+ self.a = k*wc/tc
+ def val(self,d):
+ return d/(d+self.a*(1.-d))
+ def deriv(self,d):
+ return self.a/(d+self.a*(1.-d))**2
+
+class damage_updater:
+ def __init__(self,Gc, sigc, k, g, h):
+ self.Gc = Gc
+ self.sigc = sigc
+ self.k = k
+ self.g = g
+ self.h = h
+ def update_d(self,d_previous,lda):
+ dz = d_previous.flatten()
+ shape_d = d_previous.shape
+ lda_lda_on_k = (lda[:,0]**2 + lda[:,1]**2)/self.k
+
+ def F(d):
+ return -0.5*np.dot(lda_lda_on_k,self.g.augmented_compliance(d))+self.Gc*np.sum(self.h.val(d))
+ def dF(d):
+ return -0.5*self.g.augmented_compliance_deriv(d)*lda_lda_on_k+self.Gc*self.h.deriv(d)
+
+ #dd = np.linspace(0.,1.,100)
+ #D = dz.copy()
+ #FF = []
+ #for i in range(len(dd)):
+ # D[0] = dd[i]
+ # FF.append(F(D))
+ #plt.plot(dd,FF)
+ #plt.show()
+ #print("dz = ",dz)
+ #print("lda_lda_on_k = ",lda_lda_on_k)
+ #print("F1 = ",-0.5*np.dot(lda_lda_on_k,self.g.augmented_compliance(dz)))
+ #print("F2 = ",self.Gc*np.sum(self.h.val(dz)))
+ #print("F(dz) = ",F(dz))
+ #print("self.k = ",self.k)
+
+ boundd = scipy.optimize.Bounds(dz, np.ones(len(dz)))
+ options = {'verbose':1}
+ resd = scipy.optimize.minimize(F, dz, bounds=boundd
+ , jac =dF
+ #, method = 'trust-constr'
+ #, options = options
+ #, tol = 1e-10
+ )
+ if not resd.success:
+ print("resd.message = ",resd.message)
+ #print("niterd = ",resd.nit)
+ return np.resize(resd.x,shape_d)
+
+def energy(model):
+ return model.getEnergy("potential") + model.getEnergy("dissipated")
+
+def set_dumpers(model):
+ model.setBaseName("cohesive")
+ model.addDumpFieldVector("displacement")
+ model.addDumpFieldVector("external_force")
+ model.addDumpField("strain")
+ model.addDumpField("stress")
+ model.addDumpField("blocked_dofs")
+
+ model.setBaseNameToDumper("cohesive elements", "cohesive")
+ model.addDumpFieldVectorToDumper("cohesive elements", "displacement")
+ model.addDumpFieldToDumper("cohesive elements", "damage")
+ model.addDumpFieldVectorToDumper("cohesive elements", "tractions")
+ model.addDumpFieldVectorToDumper("cohesive elements", "opening")
+
+
+def solve(material_file, mesh_file, increment):
+ aka.parseInput(material_file)
+ spatial_dimension = 2
+
+ # -------------------------------------------------------------------------
+ # Initialization
+ # -------------------------------------------------------------------------
+ mesh = aka.Mesh(spatial_dimension)
+ mesh.read(mesh_file)
+
+ model = aka.SolidMechanicsModelCohesive(mesh)
+ model.getElementInserter().setLimit(aka._x, -0.1, 0.1);
+
+ model.initFull(_analysis_method=aka._static, _is_extrinsic=False)
+
+ set_dumpers(model)
+
+ L = 0.4
+ E = model.getMaterial("steel").getReal("E")
+ Gc = model.getMaterial("cohesive").getReal("G_c")
+ tc = model.getMaterial("cohesive").getReal("sigma_c")
+ wc = 2.*Gc/tc
+ k = model.getMaterial("cohesive").getReal("k")
+ g = degradation_function_linear_czm()
+ h = softening_function_linear_czm(k,tc,Gc)
+ updater = damage_updater(Gc,tc,k,g,h)
+
+ imposed_disp = [0.]
+ reaction_force = [0.]
+ U = 0.
+ F = 0.
+
+ # -------------------------------------------------------------------------
+ # Boundary conditions
+ # -------------------------------------------------------------------------
+ model.applyBC(aka.FixedValue(0.0, aka._x), "left")
+ model.applyBC(aka.FixedValue(0.0, aka._y), "point")
+ model.applyBC(aka.IncrementValue(increment, aka._x), "right")
+ U+=increment
+ imposed_disp.append(U)
+
+ model.getExternalForce()[:] = 0
+
+ solver = model.getNonLinearSolver("static")
+ solver.set("max_iterations", 2)
+ solver.set("threshold", 1e-10)
+ solver.set("convergence_type", aka.SolveConvergenceCriteria.residual)
+
+ model.getNewSolver("linear_static", aka.TimeStepSolverType.static,
+ aka.NonLinearSolverType.linear)
+ model.setIntegrationScheme("linear_static", "displacement",
+ aka.IntegrationSchemeType.pseudo_time)
+ model.setIntegrationScheme("linear_static", "lambda",
+ aka.IntegrationSchemeType.pseudo_time)
+
+ d = model.getMaterial("cohesive").getInternalReal("czm_damage")(aka._cohesive_2d_4)
+ lda = model.getMaterial("cohesive").getInternalReal("lambda")(aka._cohesive_2d_4)
+ Fint = model.getInternalForce()
+
+ nodes_right = mesh.getElementGroup("right").getNodeGroup().getNodes()
+ model.solveStep("linear_static")
+ F = -np.sum(Fint[nodes_right,0])
+ reaction_force.append(F)
+
+ model.dump()
+ #model.dump("cohesive elements")
+
+ maxsteps = 8
+ maxiter = 10
+ tol = 1e-6
+
+ for i in range(0, maxsteps):
+ print("---- step {0}/{1}".format(i, maxsteps))
+ model.applyBC(aka.IncrementValue(increment, aka._x), "right")
+ U+=increment
+ imposed_disp.append(U)
+ it = 0
+ d_previous = model.getMaterial("cohesive").getInternalReal("czm_damage")(aka._cohesive_2d_4).copy()
+ W0 = energy(model)
+ while it < maxiter:
+ #aka.debug.setDebugLevel(aka.dblInfo);
+ model.solveStep("linear_static")
+ d[:] = updater.update_d(d_previous,lda)
+ W = energy(model)
+ err = abs(W-W0)/W
+ print("iter ",it,", err = ",err)
+ if err < tol:
+ break
+ else:
+ W0 = W
+ it+=1
+ F = -np.sum(Fint[nodes_right,0])
+ reaction_force.append(F)
+
+ model.dump()
+ #model.dump("cohesive elements")
+
+ plt.plot(imposed_disp,reaction_force,'-o')
+ plt.plot([0., (tc/E)*L, wc], [0., tc*L, 0.], 'k-o')
+ plt.show()
+ W = model.getEnergy("potential");
+ D = model.getEnergy("dissipated");
+ print("Elastic energy = ",W)
+ print("Dissipated energy = ",D)
+
+# -----------------------------------------------------------------------------
+# main
+# -----------------------------------------------------------------------------
+def main():
+ mesh_file = "triangle.msh"
+ material_file = "material.dat"
+ increment = 0.01e-3
+ solve(material_file, mesh_file, increment)
+
+
+# -----------------------------------------------------------------------------
+if __name__ == "__main__":
+ main()
diff --git a/packages/cohesive_element.cmake b/packages/cohesive_element.cmake
index 8914b039e..a735cb8ce 100644
--- a/packages/cohesive_element.cmake
+++ b/packages/cohesive_element.cmake
@@ -1,74 +1,79 @@
#===============================================================================
# Copyright (©) 2012-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This file is part of Akantu
#
# Akantu is free software: you can redistribute it and/or modify it under the
# terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with Akantu. If not, see .
#
#===============================================================================
package_declare(cohesive_element DEFAULT ON
DESCRIPTION "Use cohesive_element package of Akantu"
DEPENDS solid_mechanics)
package_declare_sources(cohesive_element
fe_engine/cohesive_element.hh
fe_engine/fe_engine_template_cohesive.cc
fe_engine/shape_cohesive.hh
fe_engine/shape_cohesive_inline_impl.hh
mesh_utils/cohesive_element_inserter.cc
mesh_utils/cohesive_element_inserter.hh
mesh_utils/cohesive_element_inserter_inline_impl.hh
mesh_utils/cohesive_element_inserter_parallel.cc
mesh_utils/cohesive_element_inserter_helper.cc
mesh_utils/cohesive_element_inserter_helper.hh
model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc
model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.hh
model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.cc
model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/cohesive_internal_field.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/cohesive_internal_field_tmpl.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_bilinear.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_fatigue.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_fatigue.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_uncoupled.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_uncoupled.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_inline_impl.hh
model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc
model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh
- model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_inline_impl.hh
-
+ model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.cc
+ model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.hh
+ model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic_inline_impl.hh
+ model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.cc
+ model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.hh
+ model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic_inline_impl.hh
+
model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc
model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.hh
model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_inline_impl.hh
model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_parallel.cc
)
diff --git a/python/py_material.cc b/python/py_material.cc
index 2b56be2d4..49e4946dd 100644
--- a/python/py_material.cc
+++ b/python/py_material.cc
@@ -1,254 +1,260 @@
/**
* Copyright (©) 2019-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "py_aka_array.hh"
#include "py_constitutive_law.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
#if defined(AKANTU_COHESIVE_ELEMENT)
#include
#include
#include
#endif
#include
/* -------------------------------------------------------------------------- */
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
/* -------------------------------------------------------------------------- */
namespace akantu {
namespace {
template class PyMaterial : public _Material {
public:
/* Inherit the constructors */
using _Material::_Material;
void initMaterial() override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE(void, _Material, initMaterial, );
};
void computeStress(ElementType el_type,
GhostType ghost_type = _not_ghost) override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE(void, _Material, computeStress, el_type, ghost_type);
}
void computeTangentModuli(ElementType el_type, Array & tangent_matrix,
GhostType ghost_type = _not_ghost) override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE(void, _Material, computeTangentModuli, el_type,
tangent_matrix, ghost_type);
}
void computePotentialEnergy(ElementType el_type) override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE(void, _Material, computePotentialEnergy, el_type);
}
[[nodiscard]] Real
getPushWaveSpeed(const Element & element) const override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE(Real, _Material, getPushWaveSpeed, element);
}
[[nodiscard]] Real
getShearWaveSpeed(const Element & element) const override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE(Real, _Material, getShearWaveSpeed, element);
}
};
/* ------------------------------------------------------------------------ */
template
void register_material_classes(py::module & mod, const std::string & name) {
py::class_<_Material, Material, PyMaterial<_Material>>(
mod, name.c_str(), py::multiple_inheritance())
.def(py::init());
}
#if defined(AKANTU_COHESIVE_ELEMENT)
// trampoline for the cohesive materials
template > * =
nullptr>
class PyMaterialCohesive : public PyMaterial<_Material> {
public:
using Parent = PyMaterial<_Material>;
/* Inherit the constructors */
using Parent::Parent;
void checkInsertion(bool check_only) override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE(void, _Material, checkInsertion, check_only);
}
void computeTraction(ElementType el_type,
GhostType ghost_type = _not_ghost) override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE(void, _Material, computeTraction, el_type, ghost_type);
}
void computeTangentTraction(ElementType el_type,
Array & tangent_matrix,
GhostType ghost_type = _not_ghost) override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE(void, _Material, computeTangentTraction, el_type,
tangent_matrix, ghost_type);
}
void computeStress(ElementType /*el_type*/,
GhostType /*ghost_type*/ = _not_ghost) final {}
void computeTangentModuli(ElementType /*el_type*/,
Array & /*tangent_matrix*/,
GhostType /*ghost_type*/ = _not_ghost) final {}
};
template
void register_material_cohesive_classes(py::module & mod,
const std::string & name) {
py::class_<_Material, MaterialCohesive, PyMaterialCohesive<_Material>>(
mod, name.c_str(), py::multiple_inheritance())
.def(py::init());
}
#endif
} // namespace
/* -------------------------------------------------------------------------- */
void register_material(py::module & mod) {
register_constitutive_law(mod);
py::class_,
ConstitutiveLaw>(mod, "Material",
py::multiple_inheritance())
.def(py::init())
.def(
"getGradU",
[](Material & self, ElementType el_type,
GhostType ghost_type = _not_ghost) -> decltype(auto) {
return self.getGradU(el_type, ghost_type);
},
py::arg("el_type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
.def(
"getStress",
[](Material & self, ElementType el_type,
GhostType ghost_type = _not_ghost) -> decltype(auto) {
return self.getStress(el_type, ghost_type);
},
py::arg("el_type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
.def(
"getPotentialEnergy",
[](Material & self, ElementType el_type) -> decltype(auto) {
return self.getPotentialEnergy(el_type);
},
py::arg("el_type"), py::return_value_policy::reference)
.def(
"getPotentialEnergy",
[](Material & self, Element element) -> Real {
return self.getPotentialEnergy(element);
},
py::arg("element"))
.def("getPotentialEnergy",
[](Material & self) -> Real { return self.getPotentialEnergy(); })
.def("initMaterial", &Material::initMaterial)
.def("getModel",
[](Material & self) -> SolidMechanicsModel & {
return self.getModel();
})
.def("getPushWaveSpeed", &Material::getPushWaveSpeed)
.def("getShearWaveSpeed", &Material::getShearWaveSpeed);
register_material_classes>(mod, "MaterialLinearElastic1D");
register_material_classes>(mod, "MaterialLinearElastic2D");
register_material_classes>(mod, "MaterialLinearElastic3D");
#if defined(AKANTU_COHESIVE_ELEMENT)
/* ------------------------------------------------------------------------ */
py::class_>(
mod, "MaterialCohesive", py::multiple_inheritance())
.def(py::init())
.def(
"getFacetFilter",
[](MaterialCohesive & self) -> decltype(auto) {
return self.getFacetFilter();
},
py::return_value_policy::reference)
.def(
"getFacetFilter",
[](MaterialCohesive & self, ElementType type,
GhostType ghost_type) -> decltype(auto) {
return self.getFacetFilter(type, ghost_type);
},
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
.def(
"getOpening",
[](MaterialCohesive & self, ElementType type, GhostType ghost_type)
-> decltype(auto) { return self.getOpening(type, ghost_type); },
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
.def(
"getTraction",
[](MaterialCohesive & self, ElementType type, GhostType ghost_type)
-> decltype(auto) { return self.getTraction(type, ghost_type); },
py::arg("type"), py::arg("ghost_type") = _not_ghost,
+ py::return_value_policy::reference)
+ .def(
+ "getDamage",
+ [](MaterialCohesive & self, ElementType type, GhostType ghost_type)
+ -> decltype(auto) { return self.getDamage(type, ghost_type); },
+ py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference);
register_material_cohesive_classes>(
mod, "MaterialCohesiveLinear2D");
register_material_cohesive_classes>(
mod, "MaterialCohesiveLinear3D");
register_material_cohesive_classes>(
mod, "MaterialCohesiveLinearFriction2D");
register_material_cohesive_classes>(
mod, "MaterialCohesiveLinearFriction3D");
#endif
py::class_(mod, "MaterialFactory")
.def_static(
"getInstance",
[]() -> MaterialFactory & { return Material::getFactory(); },
py::return_value_policy::reference)
.def("registerAllocator",
[](MaterialFactory & self, const std::string id, py::function func) {
self.registerAllocator(
id,
[func, id](Int dim, const ID & /*unused*/,
SolidMechanicsModel & model,
const ID & option) -> std::unique_ptr {
py::object obj = func(dim, id, model, option);
auto & ptr = py::cast(obj);
obj.release();
return std::unique_ptr(&ptr);
});
})
.def("getPossibleAllocators", &MaterialFactory::getPossibleAllocators);
}
} // namespace akantu
diff --git a/python/py_solid_mechanics_model_cohesive.cc b/python/py_solid_mechanics_model_cohesive.cc
index d4d8c8db9..2ff63c793 100644
--- a/python/py_solid_mechanics_model_cohesive.cc
+++ b/python/py_solid_mechanics_model_cohesive.cc
@@ -1,119 +1,120 @@
/**
* Copyright (©) 2020-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "py_aka_array.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
/* -------------------------------------------------------------------------- */
namespace akantu {
#define def_function_nocopy(func_name) \
def( \
#func_name, \
- [](SolidMechanicsModel & self) -> decltype(auto) { \
+ [](SolidMechanicsModelCohesive & self) -> decltype(auto) { \
return self.func_name(); \
}, \
py::return_value_policy::reference)
#define def_function(func_name) \
- def(#func_name, [](SolidMechanicsModel & self) -> decltype(auto) { \
+ def(#func_name, [](SolidMechanicsModelCohesive & self) -> decltype(auto) { \
return self.func_name(); \
})
void register_solid_mechanics_model_cohesive(py::module & mod) {
py::class_(mod, "CohesiveElementInserter")
.def("setLimit", &CohesiveElementInserter::setLimit)
.def(
"getCheckFacets",
[](CohesiveElementInserter & self) -> decltype(auto) {
return self.getCheckFacets();
},
py::return_value_policy::reference)
.def(
"getCheckFacets",
[](CohesiveElementInserter & self, ElementType type,
GhostType ghost_type) -> decltype(auto) {
return self.getCheckFacets(type, ghost_type);
},
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
.def(
"getInsertionFacets",
[](CohesiveElementInserter & self) -> decltype(auto) {
return self.getInsertionFacetsByElement();
},
py::return_value_policy::reference)
.def(
"getInsertionFacets",
[](CohesiveElementInserter & self, ElementType type,
GhostType ghost_type) -> decltype(auto) {
return self.getInsertionFacets(type, ghost_type);
},
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
.def("addPhysicalSurface", &CohesiveElementInserter::addPhysicalSurface)
.def("addPhysicalVolume", &CohesiveElementInserter::addPhysicalVolume);
py::class_(
mod, "SolidMechanicsModelCohesiveOptions")
.def(py::init(),
py::arg("analysis_method") = _explicit_lumped_mass,
py::arg("is_extrinsic") = false);
py::class_(
mod, "SolidMechanicsModelCohesive")
.def(py::init(), py::arg("mesh"),
py::arg("spatial_dimension") = _all_dimensions,
py::arg("id") = "solid_mechanics_model")
.def(
"initFull",
[](SolidMechanicsModel & self, const AnalysisMethod & analysis_method,
bool is_extrinsic) {
self.initFull(_analysis_method = analysis_method,
_is_extrinsic = is_extrinsic);
},
py::arg("_analysis_method") = _explicit_lumped_mass,
py::arg("_is_extrinsic") = false)
.def("checkCohesiveStress",
&SolidMechanicsModelCohesive::checkCohesiveStress)
.def("getElementInserter",
&SolidMechanicsModelCohesive::getElementInserter,
py::return_value_policy::reference)
.def("getStressOnFacets", &SolidMechanicsModelCohesive::getStressOnFacets,
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
.def("getTangents", &SolidMechanicsModelCohesive::getTangents,
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
+ .def_function_nocopy(getLambda)
.def("updateAutomaticInsertion",
&SolidMechanicsModelCohesive::updateAutomaticInsertion);
}
} // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh
index e93f23105..844b72b15 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh
@@ -1,244 +1,244 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "material.hh"
#include "solid_mechanics_model_cohesive.hh"
/* -------------------------------------------------------------------------- */
#include "cohesive_internal_field.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_MATERIAL_COHESIVE_HH_
#define AKANTU_MATERIAL_COHESIVE_HH_
/* -------------------------------------------------------------------------- */
namespace akantu {
class SolidMechanicsModelCohesive;
}
namespace akantu {
class MaterialCohesive : public Material {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
- using MyFEEngineCohesiveType =
- FEEngineTemplate;
+ using MyFEEngineCohesiveType =
+ FEEngineTemplate;
public:
MaterialCohesive(SolidMechanicsModel & model, const ID & id = "");
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// compute tractions (including normals and openings)
void computeTraction(GhostType ghost_type = _not_ghost);
/// assemble residual
void assembleInternalForces(GhostType ghost_type = _not_ghost) override;
/// check stress for cohesive elements' insertion, by default it
/// also updates the cohesive elements' data
virtual void checkInsertion(bool /*check_only*/ = false) {
AKANTU_TO_IMPLEMENT();
}
/// interpolate stress on given positions for each element (empty
/// implemantation to avoid the generic call to be done on cohesive elements)
virtual void interpolateStress(const ElementType /*type*/,
Array & /*result*/) {}
/// compute the stresses
void computeAllStresses(GhostType /*ghost_type*/ = _not_ghost) final{};
// add the facet to be handled by the material
Idx addFacet(const Element & element);
template
inline decltype(auto) getArguments(ElementType element_type,
GhostType ghost_type) {
return zip(
"normal"_n = make_view(this->normals(element_type, ghost_type)),
"opening"_n = make_view(this->opening(element_type, ghost_type)),
"traction"_n =
make_view(this->tractions(element_type, ghost_type)),
"contact_opening"_n =
make_view(this->contact_opening(element_type, ghost_type)),
"contact_traction"_n =
make_view(this->contact_tractions(element_type, ghost_type)),
"previous_opening"_n =
make_view(this->opening.previous(element_type, ghost_type)),
"previous_traction"_n =
make_view(this->tractions.previous(element_type, ghost_type)),
"delta_max"_n = this->delta_max(element_type, ghost_type),
"damage"_n = this->damage(element_type, ghost_type));
}
protected:
virtual void computeTangentTraction(ElementType /*el_type*/,
Array & /*tangent_matrix*/,
GhostType /*ghost_type*/ = _not_ghost) {
AKANTU_TO_IMPLEMENT();
}
/// compute the normal
void computeNormal(const Array & position, Array & normal,
ElementType type, GhostType ghost_type);
/// compute the opening
void computeOpening(const Array & displacement, Array & opening,
ElementType type, GhostType ghost_type);
template
void computeNormal(const Array & position, Array & normal,
GhostType ghost_type);
/// assemble stiffness
void assembleStiffnessMatrix(GhostType ghost_type) override;
/// constitutive law
virtual void computeTraction(ElementType /*el_type*/,
GhostType /*ghost_type*/ = _not_ghost) {
AKANTU_TO_IMPLEMENT();
}
/// parallelism functions
[[nodiscard]] inline Int
getNbData(const Array & elements,
const SynchronizationTag & tag) const override;
inline void packData(CommunicationBuffer & buffer,
const Array & elements,
const SynchronizationTag & tag) const override;
inline void unpackData(CommunicationBuffer & buffer,
const Array & elements,
const SynchronizationTag & tag) override;
protected:
void updateEnergies(ElementType el_type) override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get the opening
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Opening, opening, Real);
/// get the traction
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Traction, tractions, Real);
/// get damage
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real);
/// get facet filter
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FacetFilter, (*facet_filter), Idx);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(FacetFilter, (*facet_filter), Idx);
AKANTU_GET_MACRO_AUTO(FacetFilter, *facet_filter);
AKANTU_GET_MACRO_AUTO_NOT_CONST(FacetFilter, *facet_filter);
// AKANTU_GET_MACRO(ElementFilter, element_filter, const
// ElementTypeMapArray &);
/// compute reversible energy
auto getReversibleEnergy() -> Real;
/// compute dissipated energy
auto getDissipatedEnergy() -> Real;
/// compute contact energy
auto getContactEnergy() -> Real;
/// get energy
auto getEnergy(const std::string & type) -> Real override;
/// return the energy (identified by id) for the provided element
// Real getEnergy(const std::string & energy_id, const Element & element
// ) override {
// return Material::getEnergy(energy_id, element);
// }
[[nodiscard]] virtual bool needLambda() const { return false; }
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// list of facets assigned to this material
std::shared_ptr> facet_filter;
/// Link to the cohesive fem object in the model
FEEngine & fem_cohesive;
private:
/// reversible energy by quadrature point
CohesiveInternalField & reversible_energy;
/// total energy by quadrature point
CohesiveInternalField & total_energy;
protected:
/// opening in all elements and quadrature points
CohesiveInternalField & opening;
/// traction in all elements and quadrature points
CohesiveInternalField & tractions;
/// traction due to contact
CohesiveInternalField & contact_tractions;
/// normal openings for contact tractions
CohesiveInternalField & contact_opening;
/// maximum displacement
CohesiveInternalField & delta_max;
/// tell if the previous delta_max state is needed (in iterative schemes)
bool use_previous_delta_max;
/// tell if the previous opening state is needed (in iterative schemes)
bool use_previous_opening;
/// damage
CohesiveInternalField & damage;
/// critical stress
FacetRandomInternalField & sigma_c;
/// pointer to the solid mechanics model for cohesive elements
SolidMechanicsModelCohesive * model;
/// critical displacement
Real delta_c{0.};
/// array to temporarily store the normals
CohesiveInternalField & normals;
};
} // namespace akantu
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "cohesive_internal_field_tmpl.hh" // NOLINT(unused-includes)
#include "material_cohesive_inline_impl.hh"
#endif /* AKANTU_MATERIAL_COHESIVE_HH_ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc
index eb312510b..d5648e028 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc
@@ -1,496 +1,70 @@
/**
* Copyright (©) 2012-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "material_cohesive_damage.hh"
-#include "dof_synchronizer.hh"
-#include "fe_engine_template.hh"
-#include "integrator_gauss.hh"
-#include "shape_cohesive.hh"
-#include "solid_mechanics_model_cohesive.hh"
-#include "sparse_matrix.hh"
+//#include "solid_mechanics_model_cohesive.hh"
/* -------------------------------------------------------------------------- */
#include
#include
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
-template
-MaterialCohesiveDamage::MaterialCohesiveDamage(SolidMechanicsModel & model,
+MaterialCohesiveDamage::MaterialCohesiveDamage(SolidMechanicsModel & model,
const ID & id)
- : MaterialCohesive(model, id),
- lambda(registerInternal("lambda",
- spatial_dimension)),
- err_openings(registerInternal(
- "err_openings", spatial_dimension)),
- czm_damage(
- registerInternal("czm_damage", 1)) {
+ : MaterialCohesive(model, id)
+{
AKANTU_DEBUG_IN();
this->registerParam("k", k, Real(0.), _pat_parsable | _pat_readable,
- "Beta parameter");
+ "Cohesive stiffness parameter");
this->registerParam("G_c", G_c, Real(0.), _pat_parsable | _pat_readable,
"Mode I fracture energy");
+ a = 2.*k*G_c/(sigma_c*sigma_c);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
-template void MaterialCohesiveDamage::initMaterial() {
+void MaterialCohesiveDamage::initMaterial() {
AKANTU_DEBUG_IN();
MaterialCohesive::initMaterial();
AKANTU_DEBUG_OUT();
}
-
/* -------------------------------------------------------------------------- */
-template
-void MaterialCohesiveDamage::assembleInternalForces(GhostType ghost_type) {
- AKANTU_DEBUG_IN();
-
- auto & internal_force = const_cast &>(model->getInternalForce());
-
- for (auto type : getElementFilter().elementTypes(spatial_dimension,
- ghost_type, _ek_cohesive)) {
- auto & elem_filter = getElementFilter(type, ghost_type);
- auto nb_element = elem_filter.size();
- if (nb_element == 0) {
- continue;
- }
-
- const auto & shapes = fem_cohesive.getShapes(type, ghost_type);
- auto & traction = tractions(type, ghost_type);
- auto & err_opening = err_openings(type, ghost_type);
-
- auto size_of_shapes = shapes.getNbComponent();
- auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
- auto nb_quadrature_points =
- fem_cohesive.getNbIntegrationPoints(type, ghost_type);
-
- /// compute @f$t_i N_a@f$
-
- auto traction_cpy = std::make_shared>(
- nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes);
-
- auto err_opening_cpy = std::make_shared>(
- nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes);
-
- auto traction_it = traction.begin(spatial_dimension, 1);
- auto err_opening_it = err_opening.begin(spatial_dimension, 1);
-
- auto shapes_filtered_begin = shapes.begin(1, size_of_shapes);
- auto traction_cpy_it =
- traction_cpy->begin(spatial_dimension, size_of_shapes);
- auto err_opening_cpy_it =
- err_opening_cpy->begin(spatial_dimension, size_of_shapes);
-
- for (Int el = 0; el < nb_element; ++el) {
- auto current_quad = elem_filter(el) * nb_quadrature_points;
-
- for (Int q = 0; q < nb_quadrature_points; ++q, ++traction_it,
- ++err_opening_it, ++current_quad, ++traction_cpy_it,
- ++err_opening_cpy_it) {
-
- auto && shapes_filtered = shapes_filtered_begin[current_quad];
-
- *traction_cpy_it = (*traction_it) * shapes_filtered;
- *err_opening_cpy_it = (*err_opening_it) * shapes_filtered;
- }
- }
-
- /**
- * compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t
- * \mathbf{t}_q \overline w_q J_q@f$
- */
- auto partial_int_t_N = std::make_shared>(
- nb_element, spatial_dimension * size_of_shapes, "int_t_N");
-
- fem_cohesive.integrate(*traction_cpy, *partial_int_t_N,
- spatial_dimension * size_of_shapes, type, ghost_type,
- elem_filter);
-
- auto int_t_N = std::make_shared>(
- nb_element, 2 * spatial_dimension * size_of_shapes, "int_t_N");
-
- auto * int_t_N_val = int_t_N->data();
- auto * partial_int_t_N_val = partial_int_t_N->data();
- for (Int el = 0; el < nb_element; ++el) {
- std::copy_n(partial_int_t_N_val, size_of_shapes * spatial_dimension,
- int_t_N_val);
- std::copy_n(partial_int_t_N_val, size_of_shapes * spatial_dimension,
- int_t_N_val + size_of_shapes * spatial_dimension);
-
- for (Int n = 0; n < size_of_shapes * spatial_dimension; ++n) {
- int_t_N_val[n] *= -1.;
- }
-
- int_t_N_val += nb_nodes_per_element * spatial_dimension;
- partial_int_t_N_val += size_of_shapes * spatial_dimension;
- }
-
- /**
- * compute @f$\int err \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t
- * \mathbf{err}_q \overline w_q J_q@f$
- */
- auto int_err_N = std::make_shared>(
- nb_element, spatial_dimension * size_of_shapes, "int_err_N");
-
- fem_cohesive.integrate(*err_opening_cpy, *int_err_N,
- spatial_dimension * size_of_shapes, type, ghost_type,
- elem_filter);
-
- /// assemble
- model->getDOFManager().assembleElementalArrayLocalArray(
- "displacement", *int_t_N, internal_force, type, ghost_type, 1,
- elem_filter);
-
- // auto lambda_connectivity = lambda_connectivities(type, ghost_type);
- auto underlying_type = Mesh::getFacetType(type);
- model->getDOFManager().assembleElementalArrayToResidual(
- "lambda", *int_err_N, underlying_type, ghost_type, 1., elem_filter);
- }
-
- AKANTU_DEBUG_OUT();
-}
-
-/* -------------------------------------------------------------------------- */
-template
-void MaterialCohesiveDamage::assembleStiffnessMatrix(
- GhostType ghost_type) {
- AKANTU_DEBUG_IN();
-
- // auto & lambda_connectivities =
- // model->getMesh().getElementalData("lambda_connectivities");
-
- for (auto type : getElementFilter().elementTypes(spatial_dimension,
- ghost_type, _ek_cohesive)) {
- auto nb_quadrature_points =
- fem_cohesive.getNbIntegrationPoints(type, ghost_type);
- auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
-
- const auto & elem_filter = getElementFilter(type, ghost_type);
- auto nb_element = elem_filter.size();
-
- if (nb_element == 0U) {
- continue;
- }
-
- const auto & shapes = fem_cohesive.getShapes(type, ghost_type);
- auto size_of_shapes = shapes.getNbComponent();
-
- auto shapes_filtered = std::make_shared>(
- nb_element * nb_quadrature_points, size_of_shapes, "filtered shapes");
-
- for (auto && data :
- zip(filter(elem_filter,
- make_view(shapes, size_of_shapes, nb_quadrature_points)),
- make_view(*shapes_filtered, size_of_shapes,
- nb_quadrature_points))) {
- std::get<1>(data) = std::get<0>(data);
- }
-
- Matrix A(spatial_dimension * size_of_shapes,
- spatial_dimension * nb_nodes_per_element);
- A.zero();
- for (Int i = 0; i < spatial_dimension * size_of_shapes; ++i) {
- A(i, i) = 1;
- A(i, i + spatial_dimension * size_of_shapes) = -1;
- }
-
- /// get the tangent matrix @f$\frac{\partial{(t/\delta)}}{\partial{\delta}}
- /// @f$
- /// TODO : optimisation not to reassemble uu term, which does not change
- /// during computation
- auto tangent_stiffness_matrix_uu = std::make_unique>(
- nb_element * nb_quadrature_points,
- spatial_dimension * spatial_dimension, "tangent_stiffness_matrix_uu");
-
- auto tangent_stiffness_matrix_ll = std::make_unique>(
- nb_element * nb_quadrature_points,
- spatial_dimension * spatial_dimension, "tangent_stiffness_matrix_ll");
-
- computeNormal(model->getCurrentPosition(), normals(type, ghost_type), type,
- ghost_type);
-
- /// compute openings @f$\mathbf{\delta}@f$
- // computeOpening(model->getDisplacement(), opening(type, ghost_type), type,
- // ghost_type);
-
- tangent_stiffness_matrix_uu->zero();
- tangent_stiffness_matrix_ll->zero();
-
- computeTangentTraction(type, *tangent_stiffness_matrix_uu,
- *tangent_stiffness_matrix_ll, ghost_type);
-
- UInt size_at_nt_duu_n_a = spatial_dimension * nb_nodes_per_element *
- spatial_dimension * nb_nodes_per_element;
- auto at_nt_duu_n_a =
- std::make_unique>(nb_element * nb_quadrature_points,
- size_at_nt_duu_n_a, "A^t*N^t*Duu*N*A");
-
- UInt size_nt_dll_n =
- spatial_dimension * size_of_shapes * spatial_dimension * size_of_shapes;
- auto nt_dll_n = std::make_unique>(
- nb_element * nb_quadrature_points, size_nt_dll_n, "N^t*Dll*N");
-
- UInt size_at_nt_dul_n = spatial_dimension * nb_nodes_per_element *
- spatial_dimension * size_of_shapes;
- auto at_nt_dul_n = std::make_unique>(
- nb_element * nb_quadrature_points, size_at_nt_dul_n, "A^t*N^t*Dul*N");
-
- Matrix N(spatial_dimension, spatial_dimension * size_of_shapes);
-
- for (auto && [At_Nt_Duu_N_A, Duu, Nt_Dll_N, Dll, At_Nt_Dul_N, shapes] :
- zip(make_view(*at_nt_duu_n_a, spatial_dimension * nb_nodes_per_element,
- spatial_dimension * nb_nodes_per_element),
- make_view(*tangent_stiffness_matrix_uu, spatial_dimension,
- spatial_dimension),
- make_view(*nt_dll_n, spatial_dimension * size_of_shapes,
- spatial_dimension * size_of_shapes),
- make_view(*tangent_stiffness_matrix_ll, spatial_dimension,
- spatial_dimension),
- make_view(*at_nt_dul_n, spatial_dimension * nb_nodes_per_element,
- spatial_dimension * size_of_shapes),
- make_view(*shapes_filtered, size_of_shapes))) {
- N.zero();
- /**
- * store the shapes in voigt notations matrix @f$\mathbf{N} =
- * \begin{array}{cccccc} N_0(\xi) & 0 & N_1(\xi) &0 & N_2(\xi) & 0 \\
- * 0 & * N_0(\xi)& 0 &N_1(\xi)& 0 & N_2(\xi) \end{array} @f$
- **/
- for (Int i = 0; i < spatial_dimension; ++i) {
- for (Int n = 0; n < size_of_shapes; ++n) {
- N(i, i + spatial_dimension * n) = shapes(n);
- }
- }
-
- /**
- * compute stiffness matrix @f$ \mathbf{K} = \delta
- *\mathbf{U}^T \int_{\Gamma_c} {\mathbf{P}^t
- *\frac{\partial{\mathbf{t}}}
- *{\partial{\delta}}
- * \mathbf{P} d\Gamma \Delta \mathbf{U}} @f$
- **/
- auto && NA = N * A;
- At_Nt_Duu_N_A = (Duu * NA).transpose() * NA;
- Nt_Dll_N = (Dll * N).transpose() * N;
- At_Nt_Dul_N = NA.transpose() * N;
- }
-
- auto Kuu_e =
- std::make_unique>(nb_element, size_at_nt_duu_n_a, "Kuu_e");
-
- fem_cohesive.integrate(*at_nt_duu_n_a, *Kuu_e, size_at_nt_duu_n_a, type,
- ghost_type, elem_filter);
-
- auto Kll_e =
- std::make_unique>(nb_element, size_nt_dll_n, "Kll_e");
-
- fem_cohesive.integrate(*nt_dll_n, *Kll_e, size_nt_dll_n, type, ghost_type,
- elem_filter);
-
- auto Kul_e =
- std::make_unique>(nb_element, size_at_nt_dul_n, "Kul_e");
-
- fem_cohesive.integrate(*at_nt_dul_n, *Kul_e, size_at_nt_dul_n, type,
- ghost_type, elem_filter);
-
- model->getDOFManager().assembleElementalMatricesToMatrix(
- "K", "displacement", *Kuu_e, type, ghost_type, _symmetric, elem_filter);
-
- auto underlying_type = Mesh::getFacetType(type);
- model->getDOFManager().assembleElementalMatricesToMatrix(
- "K", "lambda", *Kll_e, underlying_type, ghost_type, _symmetric,
- elem_filter);
-
- auto connectivity = model->getMesh().getConnectivity(type, ghost_type);
- auto conn = make_view(connectivity, connectivity.getNbComponent()).begin();
-
- auto lambda_connectivity =
- model->getLambdaMesh().getConnectivity(underlying_type, ghost_type);
- auto lambda_conn =
- make_view(lambda_connectivity, lambda_connectivity.getNbComponent())
- .begin();
-
- /// Assemble Kll_e
- // TermsToAssemble term_ll("lambda", "lambda");
- // auto el_mat_it_ll = Kll_e->begin(spatial_dimension * size_of_shapes,
- // spatial_dimension * size_of_shapes);
-
- // auto compute_ll = [&](const auto & el) {
- // auto kll_e = *el_mat_it_ll;
- // auto && lda_conn_el = lambda_conn[el];
- // auto N = lda_conn_el.rows();
- // for (Int m = 0; m < N; ++m) {
- // auto ldai = lda_conn_el(m);
- // for (Int n = m; n < N; ++n) {
- // auto ldaj = lda_conn_el(n);
- // for (Int k = 0; k < spatial_dimension; ++k) {
- // for (Int l = 0; l < spatial_dimension; ++l) {
- // auto && term_ll_ij = term_ll(ldai * spatial_dimension + k,
- // ldaj * spatial_dimension + l);
- // term_ll_ij =
- // kll_e(m * spatial_dimension + k, n * spatial_dimension +
- // l);
- // }
- // }
- // }
- // }
- // ++el_mat_it_ll;
- // };
- // for_each_element(nb_element, elem_filter, compute_ll);
-
- // model->getDOFManager().assemblePreassembledMatrix("K", term_ll);
- // model->getDOFManager().getMatrix("K").saveMatrix("Kuu_terms.mtx");
-
- /// Assemble Klu_e
- TermsToAssemble term_ul("displacement", "lambda");
- auto el_mat_it_ul = Kul_e->begin(spatial_dimension * nb_nodes_per_element,
- spatial_dimension * size_of_shapes);
-
- auto compute_ul = [&](const auto & el) {
- auto kul_e = *el_mat_it_ul;
- auto && u_conn_el = conn[el];
- auto && lda_conn_el = lambda_conn[el];
- auto M = u_conn_el.rows();
- auto N = lda_conn_el.rows();
- for (Int m = 0; m < M; ++m) {
- for (Int n = 0; n < N; ++n) {
- auto u = u_conn_el(m);
- auto lda = lda_conn_el(n);
- for (Int k = 0; k < spatial_dimension; ++k) {
- for (Int l = 0; l < spatial_dimension; ++l) {
- auto && term_ul_ij = term_ul(u * spatial_dimension + k,
- lda * spatial_dimension + l);
- term_ul_ij =
- kul_e(m * spatial_dimension + k, n * spatial_dimension + l);
- }
- }
- }
- }
- ++el_mat_it_ul;
- };
- for_each_element(nb_element, elem_filter, compute_ul);
-
- model->getDOFManager().assemblePreassembledMatrix("K", term_ul);
- }
-
- AKANTU_DEBUG_OUT();
+Real MaterialCohesiveDamage::stiffness(Real d)
+{
+ return k*(1./d-1.);
}
-
-/* -------------------------------------------------------------------------- */
-template
-void MaterialCohesiveDamage::computeLambdaOnQuad(ElementType type,
- GhostType ghost_type) {
- auto & fem_lambda = this->model->getFEEngine("LambdaFEEngine");
- const auto & lambda = this->model->getLambda();
- auto & lambda_on_quad = this->lambda(type, ghost_type);
-
-// std::cout << "lambda in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl;
-// lambda.printself(std::cout << std::endl);
-// ArrayPrintHelper::print_content(lambda,std::cout,0);
-
-
- auto underlying_type = Mesh::getFacetType(type);
- fem_lambda.interpolateOnIntegrationPoints(
- lambda, lambda_on_quad, dim, underlying_type, ghost_type,
- this->getElementFilter(type, ghost_type));
-
-// std::cout << "lambda_on_quad in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl;
-// lambda_on_quad.printself(std::cout << std::endl);
-// ArrayPrintHelper::print_content(lambda_on_quad,std::cout,0);
-
-// std::cout << "shapes_lambda in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl;
-// auto shapes_lambda = fem_lambda.getShapes(underlying_type,ghost_type,0);
-// ArrayPrintHelper::print_content(shapes_lambda,std::cout,0);
-
-// std::cout << "shapes_cohesive in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl;
-// auto & fem_cohesive =
-// this->model->getFEEngineClass("CohesiveFEEngine");
-// auto shapes_cohesive = fem_cohesive.getShapes(type,ghost_type,0);
-// ArrayPrintHelper::print_content(shapes_cohesive,std::cout,0);
+Real MaterialCohesiveDamage::augmented_stiffness(Real d)
+{
+ return k*1./d;
}
-
-/* -------------------------------------------------------------------------- */
-template
-void MaterialCohesiveDamage::computeTraction(ElementType el_type,
- GhostType ghost_type) {
-
- for (const auto & type : getElementFilter().elementTypes(
- spatial_dimension, ghost_type, _ek_cohesive)) {
- auto & elem_filter = getElementFilter(type, ghost_type);
- auto nb_element = elem_filter.size();
- if (nb_element == 0) {
- continue;
- }
-
- /// compute normals @f$\mathbf{n}@f$
- computeNormal(model->getCurrentPosition(), normals(type, ghost_type), type,
- ghost_type);
-
- /// compute openings @f$\mathbf{\delta}@f$
- computeOpening(model->getDisplacement(), opening(type, ghost_type), type,
- ghost_type);
-
- computeLambdaOnQuad(el_type, ghost_type);
-
-// auto & traction = tractions(type, ghost_type);
-// std::cout << "traction in : MaterialCohesiveDamage::computeTraction " << std::endl;
-// ArrayPrintHelper::print_content(traction,std::cout,0);
-
-// auto & lambda_ = lambda(type, ghost_type);
-// std::cout << "lambda in : MaterialCohesiveDamage::computeTraction " << std::endl;
-// ArrayPrintHelper::print_content(lambda_,std::cout,0);
-
- for (auto && args : getArguments(el_type, ghost_type)) {
- this->computeTractionOnQuad(args);
- }
- }
+Real MaterialCohesiveDamage::augmented_compliance(Real d)
+{
+ return d/k;
}
-/* -------------------------------------------------------------------------- */
-template
-void MaterialCohesiveDamage::computeTangentTraction(
- ElementType el_type, Array & tangent_matrix_uu,
- Array & tangent_matrix_ll, GhostType ghost_type) {
- AKANTU_DEBUG_IN();
-
- for (auto && [args, tangent_uu, tangent_ll] :
- zip(getArguments(el_type, ghost_type),
- make_view(tangent_matrix_uu),
- make_view(tangent_matrix_ll))) {
- computeTangentTractionOnQuad(tangent_uu, tangent_ll, args);
- }
-
- AKANTU_DEBUG_OUT();
}
-
-/* -------------------------------------------------------------------------- */
-template class MaterialCohesiveDamage<1>;
-template class MaterialCohesiveDamage<2>;
-template class MaterialCohesiveDamage<3>;
-const bool material_is_alocated_cohesive_damage [[maybe_unused]] =
- instantiateMaterial("cohesive_damage");
-
-} // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh
index 965cd4e85..538ad747d 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh
@@ -1,126 +1,77 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "material_cohesive.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_
#define AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_
namespace akantu {
/**
- * Cohesive material linear damage for extrinsic case
+ * Cohesive material linear damage
*
* parameters in the material files :
* - sigma_c : critical stress sigma_c (default: 0)
- * - beta : weighting parameter for sliding and normal opening (default:
- * 0)
- * - G_cI : fracture energy for mode I (default: 0)
- * - G_cII : fracture energy for mode II (default: 0)
- * - penalty : stiffness in compression to prevent penetration
+ * - G_c : fracture energy for mode I (default: 0)
+ * - k : cohesive stiffness
*/
-template class MaterialCohesiveDamage : public MaterialCohesive {
+class MaterialCohesiveDamage : public MaterialCohesive {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MaterialCohesiveDamage(SolidMechanicsModel & model, const ID & id = "");
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// initialize the material parameters
void initMaterial() override;
- /// assemble stiffness
- void assembleStiffnessMatrix(GhostType ghost_type) override;
+ /// these functions should be put in a separate class
+ Real stiffness(Real d);
+ Real augmented_stiffness(Real d);
+ Real augmented_compliance(Real d);
- /// assemble residual
- void assembleInternalForces(GhostType ghost_type = _not_ghost) override;
-
-protected:
- void computeLambdaOnQuad(ElementType type, GhostType ghost_type);
-
- inline decltype(auto) getArguments(ElementType element_type,
- GhostType ghost_type) {
- using namespace tuple;
- return zip_append(
- MaterialCohesive::getArguments(element_type, ghost_type),
- "lambda"_n = make_view(this->lambda(element_type, ghost_type)),
- "err_opening"_n =
- make_view(this->err_openings(element_type, ghost_type)));
- }
-
- /// constitutive law
- void computeTraction(ElementType el_type,
- GhostType ghost_type = _not_ghost) override;
-
- /// compute tangent stiffness matrix
- /// WARNING : override was removed, not sure it should have
- void computeTangentTraction(ElementType el_type,
- Array & tangent_matrix_uu,
- Array & tangent_matrix_ll,
- GhostType ghost_type);
-
- /// compute the traction for a given quadrature point
- template inline void computeTractionOnQuad(Args && args);
-
- template
- inline void
- computeTangentTractionOnQuad(Eigen::MatrixBase & tangent_uu,
- Eigen::MatrixBase & tangent_ll,
- Args && args);
-
- [[nodiscard]] bool needLambda() const override { return true; }
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// cohesive stiffness
Real k;
/// mode I fracture energy
Real G_c;
- /// augmented lagrange multiplier
- CohesiveInternalField & lambda;
-
- /// target opening
- CohesiveInternalField & err_openings;
-
- /// cohesive damage
- CohesiveInternalField & czm_damage;
-
- Vector normal_opening;
- Real normal_opening_norm{0.};
+ /// parameter for softening function
+ Real a;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
} // namespace akantu
-#include "material_cohesive_damage_inline_impl.hh"
-
#endif /* AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.cc
new file mode 100644
index 000000000..abdf02aba
--- /dev/null
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.cc
@@ -0,0 +1,169 @@
+/**
+ * Copyright (©) 2012-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
+ * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ *
+ * This file is part of Akantu
+ *
+ * Akantu is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU Lesser General Public License as published by the Free
+ * Software Foundation, either version 3 of the License, or (at your option) any
+ * later version.
+ *
+ * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+ * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with Akantu. If not, see .
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "material_cohesive_damage_extrinsic.hh"
+#include "dof_synchronizer.hh"
+#include "fe_engine_template.hh"
+#include "integrator_gauss.hh"
+#include "shape_cohesive.hh"
+#include "solid_mechanics_model_cohesive.hh"
+#include "sparse_matrix.hh"
+/* -------------------------------------------------------------------------- */
+#include
+#include
+/* -------------------------------------------------------------------------- */
+
+namespace akantu {
+
+/* -------------------------------------------------------------------------- */
+template
+MaterialCohesiveDamageExtrinsic::MaterialCohesiveDamageExtrinsic(SolidMechanicsModel & model,
+ const ID & id)
+ : MaterialCohesiveDamage(model, id),
+ lambda(registerInternal("lambda",
+ spatial_dimension)),
+ err_openings(registerInternal(
+ "err_openings", spatial_dimension)),
+ czm_damage(
+ registerInternal("czm_damage", 1))
+{
+ if(not this->model->getIsExtrinsic())
+ {
+ AKANTU_EXCEPTION(
+ "MaterialCohesiveDamageExtrinsic can be used only with extrinsic cohesive elements");
+ }
+}
+
+/* -------------------------------------------------------------------------- */
+template // NOLINTNEXTLINE(readability-function-cognitive-complexity)
+void MaterialCohesiveDamageExtrinsic::checkInsertion(bool check_only) {
+ AKANTU_DEBUG_IN();
+
+// const auto & mesh_facets = model->getMeshFacets();
+// auto & inserter = model->getElementInserter();
+
+// for (const auto & type_facet : mesh_facets.elementTypes(dim - 1)) {
+// auto type_cohesive = FEEngine::getCohesiveElementType(type_facet);
+// auto & f_insertion = inserter.getInsertionFacets(type_facet);
+// auto & damage = czm_damage(type_cohesive);
+
+// const auto & facets_check = inserter.getCheckFacets(type_facet);
+// const auto & facet_filter_array = getFacetFilter(type_facet);
+
+// auto nb_quad_facet =
+// model->getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet);
+
+//#ifndef AKANTU_NDEBUG
+// auto nb_quad_cohesive = model->getFEEngine("CohesiveFEEngine")
+// .getNbIntegrationPoints(type_cohesive);
+
+// AKANTU_DEBUG_ASSERT(nb_quad_cohesive == nb_quad_facet,
+// "The cohesive element and the corresponding facet do "
+// "not have the same numbers of integration points");
+//#endif
+
+// // loop over each facet belonging to this material
+// for (auto && [facet] : facet_filter_array) {
+// // skip facets where check shouldn't be realized
+// if (!facets_check(facet)) {
+// continue;
+// }
+
+// Vector d_check(nb_quad_facet);
+// // compute the effective norm on each quadrature point of the facet
+// for (Int q = 0; q < nb_quad_facet; ++q) {
+// auto current_quad = facet * nb_quad_facet + q;
+// auto && normal = normal_begin[current_quad];
+// auto && tangent = tangent_begin[current_quad];
+// auto && facet_stress = facet_stress_begin[current_quad];
+
+// // compute average stress on the current quadrature point
+// auto && d_1 = facet_stress(0);
+// auto && d_2 = facet_stress(1);
+
+// auto && stress_avg = (stress_1 + stress_2) / 2.;
+
+// // compute normal and effective stress
+// stress_check(q) = computeEffectiveNorm(stress_avg, normal, tangent,
+// normal_traction(q));
+// }
+
+// // verify if the effective stress overcomes the threshold
+// auto final_stress = stress_check.mean();
+// if (max_quad_stress_insertion) {
+// final_stress =
+// *std::max_element(stress_check.begin(), stress_check.end());
+// }
+
+// if (final_stress > sigma_limit) {
+// f_insertion(facet) = true;
+
+// if (check_only) {
+// continue;
+// }
+// }
+// }
+// }
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+template
+void MaterialCohesiveDamageExtrinsic::computeLambdaOnQuad(ElementType type,
+ GhostType ghost_type) {
+ auto & fem_lambda = model->getFEEngine("LambdaFEEngine");
+ const auto & lambda = model->getLambda();
+ auto & lambda_on_quad = this->lambda(type, ghost_type);
+
+ auto underlying_type = Mesh::getFacetType(type);
+ fem_lambda.interpolateOnIntegrationPoints(
+ lambda, lambda_on_quad, dim, underlying_type, ghost_type,
+ this->getElementFilter(type, ghost_type));
+}
+
+/* -------------------------------------------------------------------------- */
+template
+void MaterialCohesiveDamageExtrinsic::computeTraction(ElementType el_type,
+ GhostType ghost_type) {
+
+ for (const auto & type : getElementFilter().elementTypes(
+ spatial_dimension, ghost_type, _ek_cohesive)) {
+ auto & elem_filter = getElementFilter(type, ghost_type);
+ auto nb_element = elem_filter.size();
+ if (nb_element == 0) {
+ continue;
+ }
+
+ for (auto && args : getArguments(el_type, ghost_type)) {
+ computeTractionOnQuad(args);
+ }
+ }
+}
+
+/* -------------------------------------------------------------------------- */
+template class MaterialCohesiveDamageExtrinsic<1>;
+template class MaterialCohesiveDamageExtrinsic<2>;
+template class MaterialCohesiveDamageExtrinsic<3>;
+const bool material_is_alocated_cohesive_damage [[maybe_unused]] =
+ instantiateMaterial("cohesive_damage_extrinsic");
+
+} // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.hh
similarity index 62%
copy from src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh
copy to src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.hh
index 965cd4e85..17d69e9da 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic.hh
@@ -1,126 +1,98 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
-#include "material_cohesive.hh"
+#include "material_cohesive_damage.hh"
/* -------------------------------------------------------------------------- */
-#ifndef AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_
-#define AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_
+#ifndef AKANTU_MATERIAL_COHESIVE_DAMAGE_EXTRINSIC_HH_
+#define AKANTU_MATERIAL_COHESIVE_DAMAGE_EXTRINSIC_HH_
namespace akantu {
/**
- * Cohesive material linear damage for extrinsic case
+ * Cohesive material linear damage for extrinsinc case
*
* parameters in the material files :
* - sigma_c : critical stress sigma_c (default: 0)
- * - beta : weighting parameter for sliding and normal opening (default:
- * 0)
- * - G_cI : fracture energy for mode I (default: 0)
- * - G_cII : fracture energy for mode II (default: 0)
- * - penalty : stiffness in compression to prevent penetration
+ * - G_c : fracture energy for mode I (default: 0)
+ * - k : cohesive stiffness
*/
-template class MaterialCohesiveDamage : public MaterialCohesive {
+template class MaterialCohesiveDamageExtrinsic : public MaterialCohesiveDamage {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
- MaterialCohesiveDamage(SolidMechanicsModel & model, const ID & id = "");
+ MaterialCohesiveDamageExtrinsic(SolidMechanicsModel & model, const ID & id = "");
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
-public:
- /// initialize the material parameters
- void initMaterial() override;
-
- /// assemble stiffness
- void assembleStiffnessMatrix(GhostType ghost_type) override;
- /// assemble residual
- void assembleInternalForces(GhostType ghost_type = _not_ghost) override;
+ // /// check stress for cohesive elements' insertion
+ void checkInsertion(bool check_only = false) override;
protected:
void computeLambdaOnQuad(ElementType type, GhostType ghost_type);
inline decltype(auto) getArguments(ElementType element_type,
GhostType ghost_type) {
using namespace tuple;
return zip_append(
MaterialCohesive::getArguments(element_type, ghost_type),
- "lambda"_n = make_view(this->lambda(element_type, ghost_type)),
- "err_opening"_n =
- make_view(this->err_openings(element_type, ghost_type)));
+ "czm_damage"_n = this->czm_damage(element_type, ghost_type));
}
/// constitutive law
void computeTraction(ElementType el_type,
GhostType ghost_type = _not_ghost) override;
- /// compute tangent stiffness matrix
- /// WARNING : override was removed, not sure it should have
- void computeTangentTraction(ElementType el_type,
- Array & tangent_matrix_uu,
- Array & tangent_matrix_ll,
- GhostType ghost_type);
-
/// compute the traction for a given quadrature point
template inline void computeTractionOnQuad(Args && args);
- template
- inline void
- computeTangentTractionOnQuad(Eigen::MatrixBase & tangent_uu,
- Eigen::MatrixBase & tangent_ll,
- Args && args);
+ [[nodiscard]] bool needLambda() const override { return false; }
- [[nodiscard]] bool needLambda() const override { return true; }
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
- /// cohesive stiffness
- Real k;
-
- /// mode I fracture energy
- Real G_c;
/// augmented lagrange multiplier
CohesiveInternalField & lambda;
/// target opening
CohesiveInternalField & err_openings;
/// cohesive damage
CohesiveInternalField & czm_damage;
Vector normal_opening;
Real normal_opening_norm{0.};
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
} // namespace akantu
-#include "material_cohesive_damage_inline_impl.hh"
+#include "material_cohesive_damage_extrinsic_inline_impl.hh"
-#endif /* AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_ */
+#endif /* AKANTU_MATERIAL_COHESIVE_DAMAGE_EXTRINSIC_HH_ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic_inline_impl.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic_inline_impl.hh
new file mode 100644
index 000000000..1914d8552
--- /dev/null
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_extrinsic_inline_impl.hh
@@ -0,0 +1,47 @@
+/**
+ * Copyright (©) 2015-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
+ * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ *
+ * This file is part of Akantu
+ *
+ * Akantu is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU Lesser General Public License as published by the Free
+ * Software Foundation, either version 3 of the License, or (at your option) any
+ * later version.
+ *
+ * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+ * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with Akantu. If not, see .
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "material_cohesive_damage_extrinsic.hh"
+
+/* -------------------------------------------------------------------------- */
+#ifndef AKANTU_MATERIAL_COHESIVE_DAMAGE_EXTRINSIC_INLINE_IMPL_HH_
+#define AKANTU_MATERIAL_COHESIVE_DAMAGE_EXTRINSIC_INLINE_IMPL_HH_
+
+/* -------------------------------------------------------------------------- */
+
+namespace akantu {
+
+/* -------------------------------------------------------------------------- */
+template
+template
+inline void MaterialCohesiveDamageExtrinsic::computeTractionOnQuad(Args && args) {
+ auto && d = args["czm_damage"_n];
+ auto && opening = args["opening"_n];
+ auto && traction = args["traction"_n];
+
+ Real A = this->augmented_stiffness(d);
+ traction = k*A*opening;
+}
+
+} // namespace akantu
+
+/* -------------------------------------------------------------------------- */
+#endif // AKANTU_MATERIAL_COHESIVE_DAMAGE_EXTRINSIC_INLINE_IMPL_HH_
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.cc
similarity index 89%
copy from src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc
copy to src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.cc
index eb312510b..8f823be64 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.cc
@@ -1,496 +1,484 @@
/**
* Copyright (©) 2012-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
-#include "material_cohesive_damage.hh"
+#include "material_cohesive_damage_intrinsic.hh"
#include "dof_synchronizer.hh"
#include "fe_engine_template.hh"
#include "integrator_gauss.hh"
#include "shape_cohesive.hh"
#include "solid_mechanics_model_cohesive.hh"
#include "sparse_matrix.hh"
/* -------------------------------------------------------------------------- */
#include
#include
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
template
-MaterialCohesiveDamage::MaterialCohesiveDamage(SolidMechanicsModel & model,
+MaterialCohesiveDamageIntrinsic::MaterialCohesiveDamageIntrinsic(SolidMechanicsModel & model,
const ID & id)
- : MaterialCohesive(model, id),
+ : MaterialCohesiveDamage(model, id),
lambda(registerInternal("lambda",
spatial_dimension)),
err_openings(registerInternal(
"err_openings", spatial_dimension)),
czm_damage(
registerInternal("czm_damage", 1)) {
- AKANTU_DEBUG_IN();
-
- this->registerParam("k", k, Real(0.), _pat_parsable | _pat_readable,
- "Beta parameter");
-
- this->registerParam("G_c", G_c, Real(0.), _pat_parsable | _pat_readable,
- "Mode I fracture energy");
-
- AKANTU_DEBUG_OUT();
-}
-
-/* -------------------------------------------------------------------------- */
-template void MaterialCohesiveDamage::initMaterial() {
- AKANTU_DEBUG_IN();
-
- MaterialCohesive::initMaterial();
-
- AKANTU_DEBUG_OUT();
+ if(this->model->getIsExtrinsic())
+ {
+ AKANTU_EXCEPTION(
+ "MaterialCohesiveDamageIntrinsic can be used only with intrinsic cohesive elements");
+ }
}
/* -------------------------------------------------------------------------- */
template
-void MaterialCohesiveDamage::assembleInternalForces(GhostType ghost_type) {
+void MaterialCohesiveDamageIntrinsic::assembleInternalForces(GhostType ghost_type) {
AKANTU_DEBUG_IN();
auto & internal_force = const_cast &>(model->getInternalForce());
for (auto type : getElementFilter().elementTypes(spatial_dimension,
ghost_type, _ek_cohesive)) {
auto & elem_filter = getElementFilter(type, ghost_type);
auto nb_element = elem_filter.size();
if (nb_element == 0) {
continue;
}
const auto & shapes = fem_cohesive.getShapes(type, ghost_type);
auto & traction = tractions(type, ghost_type);
auto & err_opening = err_openings(type, ghost_type);
auto size_of_shapes = shapes.getNbComponent();
auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
auto nb_quadrature_points =
fem_cohesive.getNbIntegrationPoints(type, ghost_type);
/// compute @f$t_i N_a@f$
auto traction_cpy = std::make_shared>(
nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes);
auto err_opening_cpy = std::make_shared>(
nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes);
auto traction_it = traction.begin(spatial_dimension, 1);
auto err_opening_it = err_opening.begin(spatial_dimension, 1);
auto shapes_filtered_begin = shapes.begin(1, size_of_shapes);
auto traction_cpy_it =
traction_cpy->begin(spatial_dimension, size_of_shapes);
auto err_opening_cpy_it =
err_opening_cpy->begin(spatial_dimension, size_of_shapes);
for (Int el = 0; el < nb_element; ++el) {
auto current_quad = elem_filter(el) * nb_quadrature_points;
for (Int q = 0; q < nb_quadrature_points; ++q, ++traction_it,
++err_opening_it, ++current_quad, ++traction_cpy_it,
++err_opening_cpy_it) {
auto && shapes_filtered = shapes_filtered_begin[current_quad];
*traction_cpy_it = (*traction_it) * shapes_filtered;
*err_opening_cpy_it = (*err_opening_it) * shapes_filtered;
}
}
/**
* compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t
* \mathbf{t}_q \overline w_q J_q@f$
*/
auto partial_int_t_N = std::make_shared>(
nb_element, spatial_dimension * size_of_shapes, "int_t_N");
fem_cohesive.integrate(*traction_cpy, *partial_int_t_N,
spatial_dimension * size_of_shapes, type, ghost_type,
elem_filter);
auto int_t_N = std::make_shared>(
nb_element, 2 * spatial_dimension * size_of_shapes, "int_t_N");
auto * int_t_N_val = int_t_N->data();
auto * partial_int_t_N_val = partial_int_t_N->data();
for (Int el = 0; el < nb_element; ++el) {
std::copy_n(partial_int_t_N_val, size_of_shapes * spatial_dimension,
int_t_N_val);
std::copy_n(partial_int_t_N_val, size_of_shapes * spatial_dimension,
int_t_N_val + size_of_shapes * spatial_dimension);
for (Int n = 0; n < size_of_shapes * spatial_dimension; ++n) {
int_t_N_val[n] *= -1.;
}
int_t_N_val += nb_nodes_per_element * spatial_dimension;
partial_int_t_N_val += size_of_shapes * spatial_dimension;
}
/**
* compute @f$\int err \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t
* \mathbf{err}_q \overline w_q J_q@f$
*/
auto int_err_N = std::make_shared>(
nb_element, spatial_dimension * size_of_shapes, "int_err_N");
fem_cohesive.integrate(*err_opening_cpy, *int_err_N,
spatial_dimension * size_of_shapes, type, ghost_type,
elem_filter);
/// assemble
model->getDOFManager().assembleElementalArrayLocalArray(
"displacement", *int_t_N, internal_force, type, ghost_type, 1,
elem_filter);
// auto lambda_connectivity = lambda_connectivities(type, ghost_type);
auto underlying_type = Mesh::getFacetType(type);
model->getDOFManager().assembleElementalArrayToResidual(
"lambda", *int_err_N, underlying_type, ghost_type, 1., elem_filter);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
-void MaterialCohesiveDamage::assembleStiffnessMatrix(
+void MaterialCohesiveDamageIntrinsic::assembleStiffnessMatrix(
GhostType ghost_type) {
AKANTU_DEBUG_IN();
// auto & lambda_connectivities =
// model->getMesh().getElementalData("lambda_connectivities");
for (auto type : getElementFilter().elementTypes(spatial_dimension,
ghost_type, _ek_cohesive)) {
auto nb_quadrature_points =
fem_cohesive.getNbIntegrationPoints(type, ghost_type);
auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const auto & elem_filter = getElementFilter(type, ghost_type);
auto nb_element = elem_filter.size();
if (nb_element == 0U) {
continue;
}
const auto & shapes = fem_cohesive.getShapes(type, ghost_type);
auto size_of_shapes = shapes.getNbComponent();
auto shapes_filtered = std::make_shared>(
nb_element * nb_quadrature_points, size_of_shapes, "filtered shapes");
for (auto && data :
zip(filter(elem_filter,
make_view(shapes, size_of_shapes, nb_quadrature_points)),
make_view(*shapes_filtered, size_of_shapes,
nb_quadrature_points))) {
std::get<1>(data) = std::get<0>(data);
}
Matrix A(spatial_dimension * size_of_shapes,
spatial_dimension * nb_nodes_per_element);
A.zero();
for (Int i = 0; i < spatial_dimension * size_of_shapes; ++i) {
A(i, i) = 1;
A(i, i + spatial_dimension * size_of_shapes) = -1;
}
/// get the tangent matrix @f$\frac{\partial{(t/\delta)}}{\partial{\delta}}
/// @f$
/// TODO : optimisation not to reassemble uu term, which does not change
/// during computation
auto tangent_stiffness_matrix_uu = std::make_unique>(
nb_element * nb_quadrature_points,
spatial_dimension * spatial_dimension, "tangent_stiffness_matrix_uu");
auto tangent_stiffness_matrix_ll = std::make_unique>(
nb_element * nb_quadrature_points,
spatial_dimension * spatial_dimension, "tangent_stiffness_matrix_ll");
computeNormal(model->getCurrentPosition(), normals(type, ghost_type), type,
ghost_type);
/// compute openings @f$\mathbf{\delta}@f$
// computeOpening(model->getDisplacement(), opening(type, ghost_type), type,
// ghost_type);
tangent_stiffness_matrix_uu->zero();
tangent_stiffness_matrix_ll->zero();
computeTangentTraction(type, *tangent_stiffness_matrix_uu,
*tangent_stiffness_matrix_ll, ghost_type);
UInt size_at_nt_duu_n_a = spatial_dimension * nb_nodes_per_element *
spatial_dimension * nb_nodes_per_element;
auto at_nt_duu_n_a =
std::make_unique>(nb_element * nb_quadrature_points,
size_at_nt_duu_n_a, "A^t*N^t*Duu*N*A");
UInt size_nt_dll_n =
spatial_dimension * size_of_shapes * spatial_dimension * size_of_shapes;
auto nt_dll_n = std::make_unique>(
nb_element * nb_quadrature_points, size_nt_dll_n, "N^t*Dll*N");
UInt size_at_nt_dul_n = spatial_dimension * nb_nodes_per_element *
spatial_dimension * size_of_shapes;
auto at_nt_dul_n = std::make_unique>(
nb_element * nb_quadrature_points, size_at_nt_dul_n, "A^t*N^t*Dul*N");
Matrix N(spatial_dimension, spatial_dimension * size_of_shapes);
for (auto && [At_Nt_Duu_N_A, Duu, Nt_Dll_N, Dll, At_Nt_Dul_N, shapes] :
zip(make_view(*at_nt_duu_n_a, spatial_dimension * nb_nodes_per_element,
spatial_dimension * nb_nodes_per_element),
make_view(*tangent_stiffness_matrix_uu, spatial_dimension,
spatial_dimension),
make_view(*nt_dll_n, spatial_dimension * size_of_shapes,
spatial_dimension * size_of_shapes),
make_view(*tangent_stiffness_matrix_ll, spatial_dimension,
spatial_dimension),
make_view(*at_nt_dul_n, spatial_dimension * nb_nodes_per_element,
spatial_dimension * size_of_shapes),
make_view(*shapes_filtered, size_of_shapes))) {
N.zero();
/**
* store the shapes in voigt notations matrix @f$\mathbf{N} =
* \begin{array}{cccccc} N_0(\xi) & 0 & N_1(\xi) &0 & N_2(\xi) & 0 \\
* 0 & * N_0(\xi)& 0 &N_1(\xi)& 0 & N_2(\xi) \end{array} @f$
**/
for (Int i = 0; i < spatial_dimension; ++i) {
for (Int n = 0; n < size_of_shapes; ++n) {
N(i, i + spatial_dimension * n) = shapes(n);
}
}
/**
* compute stiffness matrix @f$ \mathbf{K} = \delta
*\mathbf{U}^T \int_{\Gamma_c} {\mathbf{P}^t
*\frac{\partial{\mathbf{t}}}
*{\partial{\delta}}
* \mathbf{P} d\Gamma \Delta \mathbf{U}} @f$
**/
auto && NA = N * A;
At_Nt_Duu_N_A = (Duu * NA).transpose() * NA;
Nt_Dll_N = (Dll * N).transpose() * N;
At_Nt_Dul_N = NA.transpose() * N;
}
auto Kuu_e =
std::make_unique>(nb_element, size_at_nt_duu_n_a, "Kuu_e");
fem_cohesive.integrate(*at_nt_duu_n_a, *Kuu_e, size_at_nt_duu_n_a, type,
ghost_type, elem_filter);
auto Kll_e =
std::make_unique>(nb_element, size_nt_dll_n, "Kll_e");
fem_cohesive.integrate(*nt_dll_n, *Kll_e, size_nt_dll_n, type, ghost_type,
elem_filter);
auto Kul_e =
std::make_unique>(nb_element, size_at_nt_dul_n, "Kul_e");
fem_cohesive.integrate(*at_nt_dul_n, *Kul_e, size_at_nt_dul_n, type,
ghost_type, elem_filter);
model->getDOFManager().assembleElementalMatricesToMatrix(
"K", "displacement", *Kuu_e, type, ghost_type, _symmetric, elem_filter);
auto underlying_type = Mesh::getFacetType(type);
model->getDOFManager().assembleElementalMatricesToMatrix(
"K", "lambda", *Kll_e, underlying_type, ghost_type, _symmetric,
elem_filter);
auto connectivity = model->getMesh().getConnectivity(type, ghost_type);
auto conn = make_view(connectivity, connectivity.getNbComponent()).begin();
auto lambda_connectivity =
model->getLambdaMesh().getConnectivity(underlying_type, ghost_type);
auto lambda_conn =
make_view(lambda_connectivity, lambda_connectivity.getNbComponent())
.begin();
/// Assemble Kll_e
// TermsToAssemble term_ll("lambda", "lambda");
// auto el_mat_it_ll = Kll_e->begin(spatial_dimension * size_of_shapes,
// spatial_dimension * size_of_shapes);
// auto compute_ll = [&](const auto & el) {
// auto kll_e = *el_mat_it_ll;
// auto && lda_conn_el = lambda_conn[el];
// auto N = lda_conn_el.rows();
// for (Int m = 0; m < N; ++m) {
// auto ldai = lda_conn_el(m);
// for (Int n = m; n < N; ++n) {
// auto ldaj = lda_conn_el(n);
// for (Int k = 0; k < spatial_dimension; ++k) {
// for (Int l = 0; l < spatial_dimension; ++l) {
// auto && term_ll_ij = term_ll(ldai * spatial_dimension + k,
// ldaj * spatial_dimension + l);
// term_ll_ij =
// kll_e(m * spatial_dimension + k, n * spatial_dimension +
// l);
// }
// }
// }
// }
// ++el_mat_it_ll;
// };
// for_each_element(nb_element, elem_filter, compute_ll);
// model->getDOFManager().assemblePreassembledMatrix("K", term_ll);
// model->getDOFManager().getMatrix("K").saveMatrix("Kuu_terms.mtx");
/// Assemble Klu_e
TermsToAssemble term_ul("displacement", "lambda");
auto el_mat_it_ul = Kul_e->begin(spatial_dimension * nb_nodes_per_element,
spatial_dimension * size_of_shapes);
auto compute_ul = [&](const auto & el) {
auto kul_e = *el_mat_it_ul;
auto && u_conn_el = conn[el];
auto && lda_conn_el = lambda_conn[el];
auto M = u_conn_el.rows();
auto N = lda_conn_el.rows();
for (Int m = 0; m < M; ++m) {
for (Int n = 0; n < N; ++n) {
auto u = u_conn_el(m);
auto lda = lda_conn_el(n);
for (Int k = 0; k < spatial_dimension; ++k) {
for (Int l = 0; l < spatial_dimension; ++l) {
auto && term_ul_ij = term_ul(u * spatial_dimension + k,
lda * spatial_dimension + l);
term_ul_ij =
kul_e(m * spatial_dimension + k, n * spatial_dimension + l);
}
}
}
}
++el_mat_it_ul;
};
for_each_element(nb_element, elem_filter, compute_ul);
model->getDOFManager().assemblePreassembledMatrix("K", term_ul);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
-void MaterialCohesiveDamage::computeLambdaOnQuad(ElementType type,
+void MaterialCohesiveDamageIntrinsic::computeLambdaOnQuad(ElementType type,
GhostType ghost_type) {
- auto & fem_lambda = this->model->getFEEngine("LambdaFEEngine");
- const auto & lambda = this->model->getLambda();
+ auto & fem_lambda = model->getFEEngine("LambdaFEEngine");
+ const auto & lambda = model->getLambda();
auto & lambda_on_quad = this->lambda(type, ghost_type);
-// std::cout << "lambda in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl;
+// std::cout << "lambda in : MaterialCohesiveDamageIntrinsic::computeLambdaOnQuad " << std::endl;
// lambda.printself(std::cout << std::endl);
// ArrayPrintHelper::print_content(lambda,std::cout,0);
auto underlying_type = Mesh::getFacetType(type);
fem_lambda.interpolateOnIntegrationPoints(
lambda, lambda_on_quad, dim, underlying_type, ghost_type,
this->getElementFilter(type, ghost_type));
-// std::cout << "lambda_on_quad in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl;
+// std::cout << "lambda_on_quad in : MaterialCohesiveDamageIntrinsic::computeLambdaOnQuad " << std::endl;
// lambda_on_quad.printself(std::cout << std::endl);
// ArrayPrintHelper::print_content(lambda_on_quad,std::cout,0);
-// std::cout << "shapes_lambda in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl;
+// std::cout << "shapes_lambda in : MaterialCohesiveDamageIntrinsic::computeLambdaOnQuad " << std::endl;
// auto shapes_lambda = fem_lambda.getShapes(underlying_type,ghost_type,0);
// ArrayPrintHelper::print_content(shapes_lambda,std::cout,0);
-// std::cout << "shapes_cohesive in : MaterialCohesiveDamage::computeLambdaOnQuad " << std::endl;
+// std::cout << "shapes_cohesive in : MaterialCohesiveDamageIntrinsic::computeLambdaOnQuad " << std::endl;
// auto & fem_cohesive =
-// this->model->getFEEngineClass("CohesiveFEEngine");
+// model->getFEEngineClass("CohesiveFEEngine");
// auto shapes_cohesive = fem_cohesive.getShapes(type,ghost_type,0);
// ArrayPrintHelper::print_content(shapes_cohesive,std::cout,0);
}
/* -------------------------------------------------------------------------- */
template
-void MaterialCohesiveDamage::computeTraction(ElementType el_type,
+void MaterialCohesiveDamageIntrinsic::computeTraction(ElementType el_type,
GhostType ghost_type) {
for (const auto & type : getElementFilter().elementTypes(
spatial_dimension, ghost_type, _ek_cohesive)) {
auto & elem_filter = getElementFilter(type, ghost_type);
auto nb_element = elem_filter.size();
if (nb_element == 0) {
continue;
}
/// compute normals @f$\mathbf{n}@f$
computeNormal(model->getCurrentPosition(), normals(type, ghost_type), type,
ghost_type);
/// compute openings @f$\mathbf{\delta}@f$
computeOpening(model->getDisplacement(), opening(type, ghost_type), type,
ghost_type);
computeLambdaOnQuad(el_type, ghost_type);
// auto & traction = tractions(type, ghost_type);
-// std::cout << "traction in : MaterialCohesiveDamage::computeTraction " << std::endl;
+// std::cout << "traction in : MaterialCohesiveDamageIntrinsic::computeTraction " << std::endl;
// ArrayPrintHelper::print_content(traction,std::cout,0);
// auto & lambda_ = lambda(type, ghost_type);
-// std::cout << "lambda in : MaterialCohesiveDamage::computeTraction " << std::endl;
+// std::cout << "lambda in : MaterialCohesiveDamageIntrinsic::computeTraction " << std::endl;
// ArrayPrintHelper::print_content(lambda_,std::cout,0);
+// lambda.printself(std::cout);
for (auto && args : getArguments(el_type, ghost_type)) {
- this->computeTractionOnQuad(args);
+ computeTractionOnQuad(args);
}
}
}
/* -------------------------------------------------------------------------- */
template
-void MaterialCohesiveDamage::computeTangentTraction(
+void MaterialCohesiveDamageIntrinsic::computeTangentTraction(
ElementType el_type, Array & tangent_matrix_uu,
Array & tangent_matrix_ll, GhostType ghost_type) {
AKANTU_DEBUG_IN();
for (auto && [args, tangent_uu, tangent_ll] :
zip(getArguments(el_type, ghost_type),
make_view(tangent_matrix_uu),
make_view(tangent_matrix_ll))) {
computeTangentTractionOnQuad(tangent_uu, tangent_ll, args);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
-template class MaterialCohesiveDamage<1>;
-template class MaterialCohesiveDamage<2>;
-template class MaterialCohesiveDamage<3>;
+template class MaterialCohesiveDamageIntrinsic<1>;
+template class MaterialCohesiveDamageIntrinsic<2>;
+template class MaterialCohesiveDamageIntrinsic<3>;
const bool material_is_alocated_cohesive_damage [[maybe_unused]] =
- instantiateMaterial("cohesive_damage");
+ instantiateMaterial("cohesive_damage_intrinsic");
} // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.hh
similarity index 82%
copy from src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh
copy to src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.hh
index 965cd4e85..ce3112be6 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive_damage_intrinsic.hh
@@ -1,126 +1,117 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
-#include "material_cohesive.hh"
+#include "material_cohesive_damage.hh"
/* -------------------------------------------------------------------------- */
-#ifndef AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_
-#define AKANTU_MATERIAL_COHESIVE_DAMAGE_HH_
+#ifndef AKANTU_MATERIAL_COHESIVE_DAMAGE_INTRINSIC_HH_
+#define AKANTU_MATERIAL_COHESIVE_DAMAGE_INTRINSIC_HH_
namespace akantu {
/**
- * Cohesive material linear damage for extrinsic case
+ * Cohesive material linear damage for intrinsinc case
*
* parameters in the material files :
* - sigma_c : critical stress sigma_c (default: 0)
- * - beta : weighting parameter for sliding and normal opening (default:
- * 0)
- * - G_cI : fracture energy for mode I (default: 0)
- * - G_cII : fracture energy for mode II (default: 0)
- * - penalty : stiffness in compression to prevent penetration
+ * - G_c : fracture energy for mode I (default: 0)
+ * - k : cohesive stiffness
*/
-template class MaterialCohesiveDamage : public MaterialCohesive {
+template class MaterialCohesiveDamageIntrinsic : public MaterialCohesiveDamage {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
- MaterialCohesiveDamage(SolidMechanicsModel & model, const ID & id = "");
+ MaterialCohesiveDamageIntrinsic(SolidMechanicsModel & model, const ID & id = "");
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
- /// initialize the material parameters
- void initMaterial() override;
-
/// assemble stiffness
void assembleStiffnessMatrix(GhostType ghost_type) override;
/// assemble residual
void assembleInternalForces(GhostType ghost_type = _not_ghost) override;
protected:
void computeLambdaOnQuad(ElementType type, GhostType ghost_type);
inline decltype(auto) getArguments(ElementType element_type,
GhostType ghost_type) {
using namespace tuple;
return zip_append(
MaterialCohesive::getArguments(element_type, ghost_type),
"lambda"_n = make_view(this->lambda(element_type, ghost_type)),
+ "czm_damage"_n = this->czm_damage(element_type, ghost_type),
"err_opening"_n =
make_view