diff --git a/examples/python/cohesive/custom_material.py b/examples/python/cohesive/custom_material.py
index 0732dee81..04b943051 100755
--- a/examples/python/cohesive/custom_material.py
+++ b/examples/python/cohesive/custom_material.py
@@ -1,194 +1,194 @@
#!/usr/bin/env python3
# pylint: disable=missing-module-docstring
# pylint: disable=missing-function-docstring
__copyright__ = (
"Copyright (©) 2022-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)"
"Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)"
)
__license__ = "LGPLv3"
import numpy as np
import akantu as aka
spatial_dimension = 2
class LinearCohesive(aka.MaterialCohesive):
"""Material Cohesive Linear."""
def __init__(self, model, _id):
"""Construct the material and register the parameters to the parser."""
super().__init__(model, _id)
super().registerParamReal(
"G_c", aka._pat_readable | aka._pat_parsable, "Fracture energy"
)
- super().registerParamReal("beta",
- aka._pat_readable | aka._pat_parsable)
+ super().registerParamReal(
+ "beta", aka._pat_readable | aka._pat_parsable, ""
+ )
self.registerInternalReal("delta_max", 1)
self.beta = 0.
self.sigma_c = 0.
self.delta_c = 0.
self.G_c = 0.
def initMaterial(self):
"""Initialize the parameters for the material."""
super().initMaterial()
self.sigma_c = self.getReal("sigma_c")
self.G_c = self.getReal("G_c")
self.beta = self.getReal("beta")
self.delta_c = 2 * self.G_c / self.sigma_c
- def checkInsertion(self, check_only):
+ def checkInsertion(self, _check_only):
"""Check if need to insert a cohesive element."""
model = self.getModel()
facets = self.getFacetFilter()
inserter = model.getElementInserter()
for type_facet in facets.elementTypes(dim=(spatial_dimension - 1)):
facet_filter = facets(type_facet)
nb_facet = facet_filter.shape[0]
if nb_facet == 0:
continue
fe_engine = model.getFEEngine("FacetsFEEngine")
facets_check = inserter.getCheckFacets(type_facet)
insertion = inserter.getInsertionFacets(type_facet)
nb_quad_facet = fe_engine.getNbIntegrationPoints(type_facet)
normals = fe_engine.getNormalsOnIntegrationPoints(type_facet)
facets_stresses = model.getStressOnFacets(type_facet).reshape(
normals.shape[0] // nb_quad_facet,
nb_quad_facet,
2,
spatial_dimension,
spatial_dimension,
)
tangents = model.getTangents(type_facet)
for facet, facet_stresses in zip(facet_filter, facets_stresses):
if not facets_check[facet]:
continue
ref_stress = 0
for q, quad_stresses in enumerate(facet_stresses):
current_quad = facet * nb_quad_facet + q
normal = normals[current_quad, :].ravel()
tangent = tangents[current_quad, :].ravel()
stress_1 = quad_stresses.T[0]
stress_2 = quad_stresses.T[1]
avg_stress = stress_1 + stress_2 / 2.0
traction = avg_stress.dot(normal)
n = traction.dot(normal)
n = max(0, n)
t = traction.dot(tangent)
ref_stress = max(
ref_stress, np.sqrt(n * n + t * t / (self.beta**2))
)
if ref_stress > self.sigma_c:
insertion[facet] = True
# constitutive law
- def computeTraction(self, normals, el_type, ghost_type):
+ def computeTraction(self, el_type, ghost_type):
"""Compute the traction for a given opening."""
+ normals = self.getNormals(el_type, ghost_type)
openings = self.getOpening(el_type, ghost_type)
tractions = self.getTraction(el_type, ghost_type)
delta_max = self.getInternalReal("delta_max")(el_type)
for el in range(normals.shape[0]):
normal = normals[el].ravel()
opening = openings[el].ravel()
delta_n = opening.dot(normal) * normal
delta_s = opening - delta_n
- delta = (
- self.beta * np.linalg.norm(delta_s) ** 2 +
+ delta = self.beta * np.linalg.norm(delta_s) ** 2 + \
np.linalg.norm(delta_n) ** 2
- )
- delta_max[el] = max(delta_max[el], delta)
+ delta_max[el] = max(delta, delta_max[el])
tractions[el, :] = (
(delta * delta_s + delta_n)
* self.sigma_c
/ delta
* (1 - delta / self.delta_c)
)
-def allocator(_dim, unused, model, _id):
+def allocator(_dim, _unused, model, _id): # NOQA
"""Register the material to the material factory."""
return LinearCohesive(model, _id)
mat_factory = aka.MaterialFactory.getInstance()
mat_factory.registerAllocator("local_cohesive", allocator)
# -------------------------------------------------------------------------
# Initialization
# -------------------------------------------------------------------------
aka.parseInput("local_material.dat")
mesh = aka.Mesh(spatial_dimension)
mesh.read("plate.msh")
model = aka.SolidMechanicsModelCohesive(mesh)
model.initFull(_analysis_method=aka._static, _is_extrinsic=True)
model.initNewSolver(aka._explicit_lumped_mass)
model.setBaseName("plate")
model.addDumpFieldVector("displacement")
model.addDumpFieldVector("external_force")
model.addDumpField("strain")
model.addDumpField("stress")
model.addDumpField("blocked_dofs")
model.setBaseNameToDumper("cohesive elements", "cohesive")
model.addDumpFieldVectorToDumper("cohesive elements", "displacement")
model.addDumpFieldToDumper("cohesive elements", "damage")
model.addDumpFieldVectorToDumper("cohesive elements", "tractions")
model.addDumpFieldVectorToDumper("cohesive elements", "opening")
# -------------------------------------------------------------------------
# Boundary conditions
# -------------------------------------------------------------------------
model.applyBC(aka.FixedValue(0.0, aka._x), "XBlocked")
model.applyBC(aka.FixedValue(0.0, aka._y), "YBlocked")
trac = np.zeros(spatial_dimension)
traction = 0.095
trac[int(aka._y)] = traction
model.getExternalForce()[:] = 0
model.applyBC(aka.FromTraction(trac), "Traction")
print("Solve for traction ", traction)
solver = model.getNonLinearSolver("static")
solver.set("max_iterations", 100)
solver.set("threshold", 1e-10)
solver.set("convergence_type", aka.SolveConvergenceCriteria.residual)
model.solveStep("static")
model.dump()
model.dump("cohesive elements")
model.setTimeStep(model.getStableTimeStep() * 0.1)
maxsteps = 100
for i in range(0, maxsteps):
print("{0}/{1}".format(i, maxsteps))
model.checkCohesiveStress()
model.solveStep("explicit_lumped")
if i % 10 == 0:
model.dump()
model.dump("cohesive elements")
diff --git a/python/py_material.cc b/python/py_material.cc
index ae02c9273..4376a2c13 100644
--- a/python/py_material.cc
+++ b/python/py_material.cc
@@ -1,366 +1,404 @@
/**
* 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 <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#include "py_aka_array.hh"
#include "py_akantu_pybind11_compatibility.hh"
/* -------------------------------------------------------------------------- */
#include <material_selector.hh>
#include <solid_mechanics_model.hh>
#if defined(AKANTU_COHESIVE_ELEMENT)
#include <material_cohesive_linear.hh>
+#include <material_cohesive_linear_friction.hh>
#include <solid_mechanics_model_cohesive.hh>
#endif
#include <material_elastic.hh>
/* -------------------------------------------------------------------------- */
#include <pybind11/operators.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
/* -------------------------------------------------------------------------- */
namespace akantu {
namespace {
template <typename _Material> class PyMaterial : public _Material {
public:
/* Inherit the constructors */
using _Material::_Material;
~PyMaterial() override = default;
void initMaterial() override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE(void, _Material, initMaterial, );
};
void computeStress(ElementType el_type,
GhostType ghost_type = _not_ghost) override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE_PURE(void, _Material, computeStress, el_type,
ghost_type);
}
void computeTangentModuli(ElementType el_type, Array<Real> & tangent_matrix,
GhostType ghost_type = _not_ghost) override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE(void, _Material, computeTangentModuli, el_type,
tangent_matrix, ghost_type);
}
void computePotentialEnergy(ElementType el_type) override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE(void, _Material, computePotentialEnergy, el_type);
}
Real getPushWaveSpeed(const Element & element) const override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE(Real, _Material, getPushWaveSpeed, element);
}
Real getShearWaveSpeed(const Element & element) const override {
// NOLINTNEXTLINE
PYBIND11_OVERRIDE(Real, _Material, getShearWaveSpeed, element);
}
template <typename T>
void registerInternal(const std::string & name, Int nb_component) {
auto && internal = std::make_shared<InternalField<T>>(name, *this);
AKANTU_DEBUG_INFO("alloc internal " << name << " "
<< &this->internals[name]);
internal->initialize(nb_component);
this->internals[name] = internal;
}
protected:
std::map<std::string, std::shared_ptr<ElementTypeMapBase>> internals;
};
/* ------------------------------------------------------------------------ */
template <typename _Material>
void register_material_classes(py::module & mod, const std::string & name) {
py::class_<_Material, Material, Parsable, PyMaterial<_Material>>(
mod, name.c_str(), py::multiple_inheritance())
- .def(py::init<SolidMechanicsModel &, const ID &>());
+ .def(py::init<SolidMechanicsModel &, const ID &>())
+ .def("registerInternalReal",
+ [](Material & self, const std::string & name, Int nb_component) {
+ return dynamic_cast<PyMaterial<_Material> &>(self)
+ .template registerInternal<Real>(name, nb_component);
+ })
+ .def("registerInternalUInt",
+ [](Material & self, const std::string & name, Int nb_component) {
+ return dynamic_cast<PyMaterial<_Material> &>(self)
+ .template registerInternal<UInt>(name, nb_component);
+ });
}
#if defined(AKANTU_COHESIVE_ELEMENT)
- template <typename _Material> class PyMaterialCohesive : public _Material {
- public:
- /* Inherit the constructors */
- using _Material::_Material;
-
- ~PyMaterialCohesive() override = default;
- void initMaterial() override {
- // NOLINTNEXTLINE
- PYBIND11_OVERRIDE(void, _Material, initMaterial, );
- };
+ // trampoline for the cohesive materials
+ template <typename _Material>
+ class PyMaterialCohesive : public PyMaterial<_Material> {
+ using Parent = PyMaterial<_Material>;
- void computeStress(ElementType /*el_type*/,
- GhostType /*ghost_type*/ = _not_ghost) final {}
+ public:
+ 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_PURE(void, _Material, computeTraction, el_type,
ghost_type);
}
- void computeTangentModuli(ElementType /*el_type*/,
- Array<Real> & /*tangent_matrix*/,
- GhostType /*ghost_type*/ = _not_ghost) final {}
-
void computeTangentTraction(ElementType el_type,
Array<Real> & 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<Real> & /*tangent_matrix*/,
+ GhostType /*ghost_type*/ = _not_ghost) final {}
+
template <typename T>
void registerInternal(const std::string & name, Int nb_component) {
auto && internal =
std::make_shared<CohesiveInternalField<T>>(name, *this);
AKANTU_DEBUG_INFO("alloc internal " << name << " "
<< &this->internals[name]);
internal->initialize(nb_component);
this->internals[name] = internal;
}
+ };
- protected:
- std::map<std::string, std::shared_ptr<ElementTypeMapBase>> internals;
+ // trampoline for the cohesive material inheritance where computeTraction is
+ // not pure virtual
+ template <typename _Material>
+ class PyMaterialCohesiveDaughters : public PyMaterialCohesive<_Material> {
+ using Parent = PyMaterialCohesive<_Material>;
+
+ public:
+ using Parent::Parent;
+
+ void computeTraction(ElementType el_type,
+ GhostType ghost_type = _not_ghost) override {
+ // NOLINTNEXTLINE
+ PYBIND11_OVERRIDE(void, _Material, computeTraction, el_type, ghost_type);
+ }
};
template <typename _Material>
void register_material_cohesive_classes(py::module & mod,
const std::string & name) {
- py::class_<_Material, Material, Parsable, PyMaterialCohesive<_Material>>(
+ py::class_<_Material, MaterialCohesive,
+ PyMaterialCohesiveDaughters<_Material>>(
mod, name.c_str(), py::multiple_inheritance())
- .def(py::init<SolidMechanicsModelCohesive &, const ID &>());
+ .def(py::init<SolidMechanicsModelCohesive &, const ID &>())
+ .def("registerInternalReal",
+ [](MaterialCohesive & self, const std::string & name,
+ UInt nb_component) {
+ auto & ref = dynamic_cast<PyMaterialCohesive<_Material> &>(self);
+ return ref.template registerInternal<Real>(name, nb_component);
+ })
+ .def("registerInternalUInt",
+ [](MaterialCohesive & self, const std::string & name,
+ UInt nb_component) {
+ return dynamic_cast<PyMaterialCohesive<_Material> &>(self)
+ .template registerInternal<UInt>(name, nb_component);
+ });
}
#endif
/* ------------------------------------------------------------------------ */
template <typename T>
void register_internal_field(py::module & mod, const std::string & name) {
py::class_<InternalField<T>, ElementTypeMapArray<T>,
std::shared_ptr<InternalField<T>>>(
mod, ("InternalField" + name).c_str());
}
} // namespace
/* -------------------------------------------------------------------------- */
void register_material(py::module & mod) {
py::class_<MaterialFactory>(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<Material> {
py::object obj = func(dim, id, model, option);
auto & ptr = py::cast<Material &>(obj);
obj.release();
return std::unique_ptr<Material>(&ptr);
});
})
.def("getPossibleAllocators", &MaterialFactory::getPossibleAllocators);
register_internal_field<Real>(mod, "Real");
register_internal_field<UInt>(mod, "UInt");
register_internal_field<Int>(mod, "Int");
py::class_<Material, Parsable, PyMaterial<Material>>(
mod, "Material", py::multiple_inheritance())
.def(py::init<SolidMechanicsModel &, const ID &>())
.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::getModel)
.def("registerInternalReal",
[](Material & self, const std::string & name, Int nb_component) {
return dynamic_cast<PyMaterial<Material> &>(self)
.registerInternal<Real>(name, nb_component);
})
.def("registerInternalUInt",
[](Material & self, const std::string & name, Int nb_component) {
return dynamic_cast<PyMaterial<Material> &>(self)
.registerInternal<UInt>(name, nb_component);
})
.def(
"getInternalReal",
[](Material & self, const ID & id) -> decltype(auto) {
return self.getInternal<Real>(id);
},
py::arg("id"), py::return_value_policy::reference)
.def(
"getInternalUInt",
[](Material & self, const ID & id) -> decltype(auto) {
return self.getInternal<UInt>(id);
},
py::arg("id"), py::return_value_policy::reference)
.def(
"getElementFilter",
[](Material & self) -> decltype(auto) {
return self.getElementFilter();
},
py::return_value_policy::reference)
/*
* These functions override the `Parsable` interface.
* This ensure that the `updateInternalParameters()` function is called.
*/
.def(
"setReal",
[](Material & self, const ID & name, const Real value) -> void {
self.setParam(name, value);
return;
},
py::arg("name"), py::arg("value"))
.def(
"setBool",
[](Material & self, const ID & name, const bool value) -> void {
self.setParam(name, value);
return;
},
py::arg("name"), py::arg("value"))
.def(
"setString",
[](Material & self, const ID & name,
const std::string & value) -> void {
self.setParam(name, value);
return;
},
py::arg("name"), py::arg("value"))
.def(
"setInt",
[](Material & self, const ID & name, const int value) -> void {
self.setParam(name, value);
return;
},
py::arg("name"), py::arg("value"))
.def("getPushWaveSpeed", &Material::getPushWaveSpeed)
.def("getShearWaveSpeed", &Material::getShearWaveSpeed)
.def("__repr__", [](Material & self) {
std::stringstream sstr;
sstr << self;
return sstr.str();
});
register_material_classes<MaterialElastic<2>>(mod, "MaterialElastic2D");
register_material_classes<MaterialElastic<3>>(mod, "MaterialElastic3D");
#if defined(AKANTU_COHESIVE_ELEMENT)
/* ------------------------------------------------------------------------ */
- py::class_<MaterialCohesive, Material, Parsable,
- PyMaterialCohesive<MaterialCohesive>>(mod, "MaterialCohesive",
- py::multiple_inheritance())
+ py::class_<MaterialCohesive, Material, PyMaterialCohesive<MaterialCohesive>>(
+ mod, "MaterialCohesive", py::multiple_inheritance())
.def(py::init<SolidMechanicsModelCohesive &, const ID &>())
.def("registerInternalReal",
[](MaterialCohesive & self, const std::string & name,
UInt nb_component) {
- return dynamic_cast<PyMaterialCohesive<MaterialCohesive> &>(self)
- .registerInternal<Real>(name, nb_component);
+ auto & ref =
+ dynamic_cast<PyMaterialCohesive<MaterialCohesive> &>(self);
+ return ref.registerInternal<Real>(name, nb_component);
})
.def("registerInternalUInt",
[](MaterialCohesive & self, const std::string & name,
UInt nb_component) {
return dynamic_cast<PyMaterialCohesive<MaterialCohesive> &>(self)
.registerInternal<UInt>(name, nb_component);
})
.def(
"getFacetFilter",
[](MaterialCohesive & self) -> decltype(auto) {
return self.getFacetFilter();
},
py::return_value_policy::reference)
.def(
"getFacetFilter",
[](MaterialCohesive & self, ElementType type,
GhostType ghost_type) -> decltype(auto) {
return self.getFacetFilter(type, ghost_type);
},
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
.def(
"getOpening",
[](MaterialCohesive & self, ElementType type, GhostType ghost_type)
-> decltype(auto) { return self.getOpening(type, ghost_type); },
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference)
.def(
"getTraction",
[](MaterialCohesive & self, ElementType type, GhostType ghost_type)
-> decltype(auto) { return self.getTraction(type, ghost_type); },
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::return_value_policy::reference);
register_material_cohesive_classes<MaterialCohesiveLinear<2>>(
mod, "MaterialCohesiveLinear2D");
register_material_cohesive_classes<MaterialCohesiveLinear<3>>(
mod, "MaterialCohesiveLinear3D");
+ register_material_cohesive_classes<MaterialCohesiveLinearFriction<2>>(
+ mod, "MaterialCohesiveLinearFriction2D");
+ register_material_cohesive_classes<MaterialCohesiveLinearFriction<3>>(
+ mod, "MaterialCohesiveLinearFriction3D");
#endif
}
} // namespace akantu
diff --git a/src/model/phase_field/phase_field_model_inline_impl.hh b/src/model/phase_field/phase_field_model_inline_impl.hh
index afc009270..0c7f2aaa6 100644
--- a/src/model/phase_field/phase_field_model_inline_impl.hh
+++ b/src/model/phase_field/phase_field_model_inline_impl.hh
@@ -1,90 +1,90 @@
/**
* 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 <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#include "aka_named_argument.hh"
#include "phasefield_selector.hh"
#include "phasefield_selector_tmpl.hh"
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_PHASE_FIELD_MODEL_INLINE_IMPL_CC__
#define __AKANTU_PHASE_FIELD_MODEL_INLINE_IMPL_CC__
namespace akantu {
/* -------------------------------------------------------------------------- */
inline decltype(auto) PhaseFieldModel::getPhaseFields() {
return make_dereference_adaptor(phasefields);
}
/* -------------------------------------------------------------------------- */
inline decltype(auto) PhaseFieldModel::getPhaseFields() const {
return make_dereference_adaptor(phasefields);
}
/* -------------------------------------------------------------------------- */
inline PhaseField & PhaseFieldModel::getPhaseField(Idx mat_index) {
- AKANTU_DEBUG_ASSERT(mat_index < phasefields.size(),
+ AKANTU_DEBUG_ASSERT(mat_index < Idx(phasefields.size()),
"The model " << id << " has no phasefield no "
<< mat_index);
return *phasefields[mat_index];
}
/* -------------------------------------------------------------------------- */
inline const PhaseField & PhaseFieldModel::getPhaseField(Idx mat_index) const {
- AKANTU_DEBUG_ASSERT(mat_index < phasefields.size(),
+ AKANTU_DEBUG_ASSERT(mat_index < Idx(phasefields.size()),
"The model " << id << " has no phasefield no "
<< mat_index);
return *phasefields[mat_index];
}
/* -------------------------------------------------------------------------- */
inline PhaseField & PhaseFieldModel::getPhaseField(const std::string & name) {
auto it = phasefields_names_to_id.find(name);
AKANTU_DEBUG_ASSERT(it != phasefields_names_to_id.end(),
"The model " << id << " has no phasefield named "
<< name);
return *phasefields[it->second];
}
/* -------------------------------------------------------------------------- */
inline Int PhaseFieldModel::getPhaseFieldIndex(const std::string & name) const {
auto it = phasefields_names_to_id.find(name);
AKANTU_DEBUG_ASSERT(it != phasefields_names_to_id.end(),
"The model " << id << " has no phasefield named "
<< name);
return it->second;
}
/* -------------------------------------------------------------------------- */
inline const PhaseField &
PhaseFieldModel::getPhaseField(const std::string & name) const {
auto it = phasefields_names_to_id.find(name);
AKANTU_DEBUG_ASSERT(it != phasefields_names_to_id.end(),
"The model " << id << " has no phasefield named "
<< name);
return *phasefields[it->second];
}
/* -------------------------------------------------------------------------- */
} // namespace akantu
#endif /* __AKANTU_PHASE_FIELD_MODEL_INLINE_IMPL_CC__ */
diff --git a/src/model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh b/src/model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh
index a5efd9ab0..e56aa8ced 100644
--- a/src/model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh
+++ b/src/model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh
@@ -1,367 +1,366 @@
/**
* 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 <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#include "aka_iterators.hh"
#include "material_anisotropic_damage.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_MATERIAL_ANISOTROPIC_DAMAGE_TMPL_HH_
#define AKANTU_MATERIAL_ANISOTROPIC_DAMAGE_TMPL_HH_
namespace akantu {
struct EmptyIteratorContainer {
using size_type = std::size_t;
struct iterator {
auto & operator++() { return *this; }
Real operator*() { return 0; }
bool operator!=(const iterator & /*unused*/) const { return true; }
bool operator==(const iterator & /*unused*/) const { return false; }
};
auto begin() const // NOLINT(readability-convert-member-functions-to-static)
{
return iterator();
}
auto end() const // NOLINT(readability-convert-member-functions-to-static)
{
return iterator();
}
};
} // namespace akantu
namespace std {
template <> struct iterator_traits<::akantu::EmptyIteratorContainer::iterator> {
using iterator_category = forward_iterator_tag;
using value_type = akantu::Real;
using difference_type = std::ptrdiff_t;
using pointer = akantu::Real *;
using reference = akantu::Real &;
};
} // namespace std
namespace akantu {
namespace {
template <Int dim, class D, class Op>
void tensorPlus_(const Eigen::MatrixBase<D> & A, Op && oper) {
Vector<Real, dim> A_eigs;
A.eig(A_eigs);
for (auto & ap : A_eigs) {
oper(ap);
}
}
template <Int dim, class D> auto tensorPlus2(const Eigen::MatrixBase<D> & A) {
Real square = 0;
tensorPlus_<dim>(A, [&](Real eig) {
eig = std::max(eig, 0.);
square += eig * eig;
});
return square;
}
template <Int dim, class D>
auto tensorPlusTrace(const Eigen::MatrixBase<D> & A) {
Real trace_plus = 0.;
Real trace_minus = 0.;
tensorPlus_<dim>(A, [&](Real eig) {
trace_plus += std::max(eig, 0.);
trace_minus += std::min(eig, 0.);
});
- return std::make_pair(trace_plus, trace_minus);
+ return std::pair(trace_plus, trace_minus);
}
template <Int dim, class Op, class D1, class D2>
auto tensorPlusOp(const Eigen::MatrixBase<D1> & A,
Eigen::MatrixBase<D2> & A_directions, Op && oper) {
Vector<Real, dim> A_eigs;
Matrix<Real, dim, dim> A_diag;
A.eig(A_eigs, A_directions);
for (auto && [i, eig] : enumerate(A_eigs)) {
A_diag(i, i) = oper(std::max(eig, 0.), i);
}
Matrix<Real, dim, dim> res =
A_directions * A_diag * A_directions.transpose();
return res;
}
template <Int dim, class D1, class D2, class Op>
auto tensorPlus(const Eigen::MatrixBase<D1> & A,
Eigen::MatrixBase<D2> & A_directions) {
return tensorPlusOp<dim>(A, A_directions, [](Real x, Real) { return x; });
}
template <Int dim, class D, class Op>
auto tensorPlusOp(const Eigen::MatrixBase<D> & A, Op && oper) {
Matrix<Real, dim, dim> A_directions;
A_directions.zero();
return tensorPlusOp<dim>(A, A_directions, std::forward<Op>(oper));
}
template <Int dim, class D> auto tensorPlus(const Eigen::MatrixBase<D> & A) {
return tensorPlusOp<dim>(A, [](Real x, Real) { return x; });
}
template <Int dim, class D> auto tensorSqrt(const Eigen::MatrixBase<D> & A) {
return tensorPlusOp<dim>(A, [](Real x, Int) { return std::sqrt(x); });
}
} // namespace
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template <Int dim, template <Int> class EquivalentStrain,
template <Int> class DamageThreshold, template <Int> class Parent>
MaterialAnisotropicDamage<dim, EquivalentStrain, DamageThreshold, Parent>::
MaterialAnisotropicDamage(SolidMechanicsModel & model, const ID & id)
: Parent<dim>(model, id), damage("damage_tensor", *this),
elastic_stress("elastic_stress", *this),
equivalent_strain("equivalent_strain", *this),
trace_damage("trace_damage", *this), equivalent_strain_function(*this),
damage_threshold_function(*this) {
this->registerParam("Dc", Dc, _pat_parsable, "Critical damage");
this->damage.initialize(dim * dim);
this->elastic_stress.initialize(dim * dim);
this->equivalent_strain.initialize(1);
this->trace_damage.initialize(1);
this->trace_damage.initializeHistory();
}
/* -------------------------------------------------------------------------- */
template <Int dim, template <Int> class EquivalentStrain,
template <Int> class DamageThreshold, template <Int> class Parent>
template <class D1, class D2, class D3>
void MaterialAnisotropicDamage<dim, EquivalentStrain, DamageThreshold, Parent>::
damageStress(Eigen::MatrixBase<D1> & sigma,
const Eigen::MatrixBase<D2> & sigma_el,
const Eigen::MatrixBase<D3> & D, Real TrD) {
// σ_(n + 1) = (1 − D_(n + 1))^(1/2) σ~_(n + 1) (1 − D_(n + 1))^(1 / 2)
// - ((1 − D_(n + 1)) : σ~_(n + 1))/ (3 - Tr(D_(n+1))) (1 − D_(n + 1))
// + 1/3 (1 - Tr(D_(n+1)) <Tr(σ~_(n + 1))>_+ + <Tr(σ~_(n + 1))>_-) I
Matrix<Real, dim, dim> one_D = Matrix<Real, dim, dim>::Identity() - D;
auto sqrt_one_D = tensorSqrt<dim>(one_D);
Real Tr_sigma_plus;
Real Tr_sigma_minus;
std::tie(Tr_sigma_plus, Tr_sigma_minus) = tensorPlusTrace<dim>(sigma_el);
auto I = Matrix<Real, dim, dim>::Identity();
sigma = sqrt_one_D * sigma_el * sqrt_one_D -
(one_D.doubleDot(sigma_el) / (dim - TrD) * one_D) +
1. / dim * ((1 - TrD) * Tr_sigma_plus - Tr_sigma_minus) * I;
}
/* -------------------------------------------------------------------------- */
template <Int dim, template <Int> class EquivalentStrain,
template <Int> class DamageThreshold, template <Int> class Parent>
template <class Args>
void MaterialAnisotropicDamage<dim, EquivalentStrain, DamageThreshold,
Parent>::computeStressOnQuad(Args && args) {
auto & sigma = args["sigma"_n];
auto & grad_u = args["grad_u"_n];
auto & sigma_el = args["sigma_el"_n];
auto & epsilon_hat = args["epsilon_hat"_n];
auto & D = args["damage"_n];
auto & TrD_n_1 = args["TrD_n_1"_n];
auto & TrD = args["TrD"_n];
auto & equivalent_strain_data = args["equivalent_strain_data"_n];
auto & damage_threshold_data = args["damage_threshold_data"_n];
Matrix<Real, dim, dim> Dtmp;
Real TrD_n_1_tmp;
Matrix<Real, dim, dim> epsilon;
// yes you read properly this is a label for a goto
auto computeDamage = [&]() {
MaterialElastic<dim>::computeStressOnQuad(args);
epsilon = this->template gradUToEpsilon<dim>(grad_u);
// evaluate the damage criteria
epsilon_hat = equivalent_strain_function(epsilon, equivalent_strain_data);
// evolve the damage if needed
auto K_TrD = damage_threshold_function.K(TrD, damage_threshold_data);
auto f = epsilon_hat - K_TrD;
// if test function > 0 evolve the damage
if (f > 0) {
TrD_n_1_tmp =
damage_threshold_function.K_inv(epsilon_hat, damage_threshold_data);
auto epsilon_plus = tensorPlus<dim>(epsilon);
auto delta_lambda = (TrD_n_1_tmp - TrD) / (epsilon_hat * epsilon_hat);
Dtmp = D + delta_lambda * epsilon_plus;
return true;
}
return false;
};
// compute a temporary version of the new damage
auto is_damage_updated = computeDamage();
if (is_damage_updated) {
/// Check and correct for broken case
if (Dtmp.trace() > Dc) {
if (epsilon.trace() > 0) { // tensile loading
auto kpa = this->kpa;
auto lambda = this->lambda;
// change kappa to Kappa_broken = (1-Dc) Kappa
kpa = (1 - Dc) * kpa;
this->E = 9 * kpa * (kpa - lambda) / (3 * kpa - lambda);
this->nu = lambda / (3 * kpa - lambda);
this->updateInternalParameters();
computeDamage();
} else if (std::abs(epsilon.trace()) < 1e-10) { // deviatoric case
Matrix<Real, dim, dim> n;
std::vector<Int> ns;
tensorPlusOp<dim>(Dtmp, n, [&](Real x, Int i) {
if (x > this->Dc) {
ns.push_back(i);
return this->Dc;
}
return x;
});
}
}
TrD_n_1 = TrD_n_1_tmp;
D = Dtmp;
} else {
TrD_n_1 = TrD;
}
// apply the damage to the stress
damageStress(sigma, sigma_el, D, TrD_n_1);
}
/* -------------------------------------------------------------------------- */
template <Int dim, template <Int> class EquivalentStrain,
template <Int> class DamageThreshold, template <Int> class Parent>
void MaterialAnisotropicDamage<dim, EquivalentStrain, DamageThreshold,
Parent>::computeStress(ElementType type,
GhostType ghost_type) {
for (auto && args : getArguments(type, ghost_type)) {
computeStressOnQuad(args);
}
}
/* -------------------------------------------------------------------------- */
/* EquivalentStrain functions */
/* -------------------------------------------------------------------------- */
template <Int dim>
class EquivalentStrainMazars : public EmptyIteratorContainer {
public:
EquivalentStrainMazars(Material & /*mat*/) {}
template <class D, class... Other>
Real operator()(const Eigen::MatrixBase<D> & epsilon, Other &&... /*other*/) {
- Real epsilon_hat = 0.;
- std::tie(epsilon_hat, std::ignore) = tensorPlusTrace<dim>(epsilon);
+ auto && [epsilon_hat, _] = tensorPlusTrace<dim>(epsilon);
return std::sqrt(epsilon_hat);
}
};
template <Int dim>
class EquivalentStrainMazarsDruckerPrager : public EquivalentStrainMazars<dim> {
public:
EquivalentStrainMazarsDruckerPrager(Material & mat)
: EquivalentStrainMazars<dim>(mat) {
mat.registerParam("k", k, _pat_parsable, "k");
}
template <class D, class... Other>
Real operator()(const Eigen::MatrixBase<D> & epsilon, Real) {
Real epsilon_hat = EquivalentStrainMazars<dim>::operator()(epsilon);
epsilon_hat += k * epsilon.trace();
return epsilon_hat;
}
protected:
Real k;
};
/* -------------------------------------------------------------------------- */
/* DamageThreshold functions */
/* -------------------------------------------------------------------------- */
template <Int dim> class DamageThresholdLinear : public EmptyIteratorContainer {
public:
DamageThresholdLinear(Material & mat) : mat(mat) {
mat.registerParam("A", A, _pat_parsable, "A");
mat.registerParam("K0", K0, _pat_parsable, "K0");
}
template <class... Other> Real K(Real x, Other &&... /*other*/) {
return 1. / A * x + K0;
}
template <class... Other> Real K_inv(Real x, Other &&... /*other*/) {
return A * (x - K0);
}
private:
Material & mat;
Real A;
Real K0;
};
template <Int dim> class DamageThresholdTan : public EmptyIteratorContainer {
public:
DamageThresholdTan(Material & mat) : mat(mat) {
mat.registerParam("a", a, _pat_parsable, "a");
mat.registerParam("A", A, _pat_parsable, "A");
mat.registerParam("K0", K0, _pat_parsable, "K0");
}
template <class... Other> Real K(Real x, Other &&... /*other*/) {
return a * std::tan(std::atan2(x, a) - std::atan2(K0, a));
}
template <class... Other> Real K_inv(Real x, Other &&... /*other*/) {
return a * A * (std::atan2(x, a) - std::atan2(K0, a));
}
private:
Material & mat;
Real a{2.93e-4};
Real A{5e3};
Real K0{5e-5};
};
} // namespace akantu
#endif /* AKANTU_MATERIAL_ANISOTROPIC_DAMAGE_TMPL_HH_ */