diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh
index 71494dbf3..02a71b592 100644
--- a/src/model/solid_mechanics/material.hh
+++ b/src/model/solid_mechanics/material.hh
@@ -1,803 +1,804 @@
/**
* 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 "aka_factory.hh"
#include "aka_voigthelper.hh"
#include "data_accessor.hh"
#include "integration_point.hh"
#include "parsable.hh"
#include "parser.hh"
/* -------------------------------------------------------------------------- */
#include "internal_field.hh"
#include "random_internal_field.hh"
/* -------------------------------------------------------------------------- */
#include "mesh_events.hh"
#include "solid_mechanics_model_event_handler.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_MATERIAL_HH_
#define AKANTU_MATERIAL_HH_
/* -------------------------------------------------------------------------- */
namespace akantu {
class Model;
class SolidMechanicsModel;
class Material;
} // namespace akantu
namespace akantu {
using MaterialFactory =
Factory;
/**
* Interface of all materials
* Prerequisites for a new material
* - inherit from this class
* - implement the following methods:
* \code
* virtual Real getStableTimeStep(Real h, const Element & element =
* ElementNull);
*
* virtual void computeStress(ElementType el_type,
* GhostType ghost_type = _not_ghost);
*
* virtual void computeTangentStiffness(ElementType el_type,
* Array & tangent_matrix,
* GhostType ghost_type = _not_ghost);
* \endcode
*
*/
class Material : public DataAccessor,
public Parsable,
public MeshEventHandler,
protected SolidMechanicsModelEventHandler {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
Material(const Material & mat) = delete;
Material & operator=(const Material & mat) = delete;
/// Initialize material with defaults
Material(SolidMechanicsModel & model, const ID & id = "");
/// Initialize material with custom mesh & fe_engine
Material(SolidMechanicsModel & model, Int dim, const Mesh & mesh,
FEEngine & fe_engine, const ID & id = "");
/// Destructor
~Material() override;
protected:
void initialize();
/* ------------------------------------------------------------------------ */
/* Function that materials can/should reimplement */
/* ------------------------------------------------------------------------ */
protected:
/// constitutive law
virtual void computeStress(ElementType /* el_type */,
GhostType /* ghost_type */ = _not_ghost) {
AKANTU_TO_IMPLEMENT();
}
/// compute the tangent stiffness matrix
virtual void computeTangentModuli(ElementType /*el_type*/,
Array & /*tangent_matrix*/,
GhostType /*ghost_type*/ = _not_ghost) {
AKANTU_TO_IMPLEMENT();
}
/// compute the potential energy
virtual void computePotentialEnergy(ElementType el_type);
/// compute the potential energy for an element
[[deprecated("Use the interface with an Element")]] virtual void
- computePotentialEnergyByElement(ElementType /*type*/, Int /*index*/,
- Vector & /*epot_on_quad_points*/) {
- AKANTU_TO_IMPLEMENT();
+ computePotentialEnergyByElement(ElementType type, Int index,
+ Vector & epot_on_quad_points) {
+ computePotentialEnergyByElement(Element{type, index, _not_ghost},
+ epot_on_quad_points);
}
virtual void
computePotentialEnergyByElement(const Element & /*element*/,
Vector & /*epot_on_quad_points*/) {
AKANTU_TO_IMPLEMENT();
}
virtual void updateEnergies(ElementType /*el_type*/) {}
virtual void updateEnergiesAfterDamage(ElementType /*el_type*/) {}
/// set the material to steady state (to be implemented for materials that
/// need it)
virtual void setToSteadyState(ElementType /*el_type*/,
GhostType /*ghost_type*/ = _not_ghost) {}
/// function called to update the internal parameters when the modifiable
/// parameters are modified
virtual void updateInternalParameters() {}
public:
/// extrapolate internal values
virtual void extrapolateInternal(const ID & id, const Element & element,
const Matrix & points,
Matrix & extrapolated);
/// compute the p-wave speed in the material
virtual Real getPushWaveSpeed(const Element & /*element*/) const {
AKANTU_TO_IMPLEMENT();
}
/// compute the s-wave speed in the material
virtual Real getShearWaveSpeed(const Element & /*element*/) const {
AKANTU_TO_IMPLEMENT();
}
/// get a material celerity to compute the stable time step (default: is the
/// push wave speed)
virtual Real getCelerity(const Element & element) const {
return getPushWaveSpeed(element);
}
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
template void registerInternal(InternalField & /*vect*/) {
AKANTU_TO_IMPLEMENT();
}
template void unregisterInternal(InternalField & /*vect*/) {
AKANTU_TO_IMPLEMENT();
}
/// initialize the material computed parameter
virtual void initMaterial();
/// compute the residual for this material
// virtual void updateResidual(GhostType ghost_type = _not_ghost);
/// assemble the residual for this material
virtual void assembleInternalForces(GhostType ghost_type);
/// save the internals in the previous_stress if needed
virtual void savePreviousState();
/// restore the internals from previous_stress if needed
virtual void restorePreviousState();
/// compute the stresses for this material
virtual void computeAllStresses(GhostType ghost_type = _not_ghost);
// virtual void
// computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost);
virtual void computeAllCauchyStresses(GhostType ghost_type = _not_ghost);
/// set material to steady state
void setToSteadyState(GhostType ghost_type = _not_ghost);
/// compute the stiffness matrix
virtual void assembleStiffnessMatrix(GhostType ghost_type);
/// add an element to the local mesh filter
inline auto addElement(ElementType type, Int element, GhostType ghost_type);
inline auto addElement(const Element & element);
/// add many elements at once
void addElements(const Array & elements_to_add);
/// remove many element at once
void removeElements(const Array & elements_to_remove);
/// function to print the contain of the class
void printself(std::ostream & stream, int indent = 0) const override;
/**
* interpolate stress on given positions for each element by means
* of a geometrical interpolation on quadrature points
*/
void interpolateStress(ElementTypeMapArray & result,
GhostType ghost_type = _not_ghost);
/**
* interpolate stress on given positions for each element by means
* of a geometrical interpolation on quadrature points and store the
* results per facet
*/
void interpolateStressOnFacets(ElementTypeMapArray & result,
ElementTypeMapArray & by_elem_result,
GhostType ghost_type = _not_ghost);
/**
* function to initialize the elemental field interpolation
* function by inverting the quadrature points' coordinates
*/
void initElementalFieldInterpolation(
const ElementTypeMapArray & interpolation_points_coordinates);
/* ------------------------------------------------------------------------ */
/* Common part */
/* ------------------------------------------------------------------------ */
protected:
/* ------------------------------------------------------------------------ */
constexpr static inline Int getTangentStiffnessVoigtSize(Int dim) {
return (dim * (dim - 1) / 2 + dim);
}
template
constexpr static inline Int getTangentStiffnessVoigtSize() {
return getTangentStiffnessVoigtSize(dim);
}
/// compute the potential energy by element
void computePotentialEnergyByElements();
/// resize the intenals arrays
virtual void resizeInternals();
template
decltype(auto) getArguments(ElementType el_type, GhostType ghost_type) {
using namespace tuple;
auto && args =
zip("grad_u"_n = make_view(this->gradu(el_type, ghost_type)),
"previous_sigma"_n =
make_view(this->stress.previous(el_type, ghost_type)),
"previous_grad_u"_n =
make_view(this->gradu.previous(el_type, ghost_type)));
if (not finite_deformation) {
return zip_append(
std::forward(args),
"sigma"_n = make_view(this->stress(el_type, ghost_type)));
}
return zip_append(std::forward(args),
"sigma"_n = make_view(
this->piola_kirchhoff_2(el_type, ghost_type)));
}
template
decltype(auto) getArgumentsTangent(Array & tangent_matrix,
ElementType el_type,
GhostType ghost_type) {
using namespace tuple;
constexpr auto tangent_size = Material::getTangentStiffnessVoigtSize(dim);
return zip("tangent_moduli"_n =
make_view(tangent_matrix),
"grad_u"_n =
make_view(this->gradu(el_type, ghost_type)));
}
/* ------------------------------------------------------------------------ */
/* Finite deformation functions */
/* This functions area implementing what is described in the paper of Bathe */
/* et al, in IJNME, Finite Element Formulations For Large Deformation */
/* Dynamic Analysis, Vol 9, 353-386, 1975 */
/* ------------------------------------------------------------------------ */
protected:
/// assemble the internal forces
template
void assembleInternalForces(ElementType type, GhostType ghost_type);
/// assemble the internal forces
template
void assembleInternalForces(GhostType ghost_type);
/// assemble the internal forces in the case of finite deformation
template
void assembleInternalForcesFiniteDeformation(ElementType type,
GhostType ghost_type);
/// assemble the internal forces in the case of finite deformation
template
void assembleInternalForcesFiniteDeformation(GhostType ghost_type);
template
void computeAllStressesFromTangentModuli(ElementType type,
GhostType ghost_type);
template
void assembleStiffnessMatrix(ElementType type, GhostType ghost_type);
/// assembling in finite deformation
template
void assembleStiffnessMatrixFiniteDeformation(ElementType type,
GhostType ghost_type);
template
void assembleStiffnessMatrix(GhostType ghost_type);
/// assembling in finite deformation
template
void assembleStiffnessMatrixNL(GhostType ghost_type);
template
void assembleStiffnessMatrixL2(GhostType ghost_type);
/* ------------------------------------------------------------------------ */
/* Conversion functions */
/* ------------------------------------------------------------------------ */
public:
/// Size of the Stress matrix for the case of finite deformation see: Bathe et
/// al, IJNME, Vol 9, 353-386, 1975
static constexpr inline Int getCauchyStressMatrixSize(Int dim) {
return (dim * dim);
}
/// Sets the stress matrix according to Bathe et al, IJNME, Vol 9, 353-386,
/// 1975
template
static constexpr inline void
setCauchyStressMatrix(const Eigen::MatrixBase & S_t,
Eigen::MatrixBase & sigma);
/// write the stress tensor in the Voigt notation.
template
static constexpr inline decltype(auto)
stressToVoigt(const Eigen::MatrixBase & stress) {
return VoigtHelper::matrixToVoigt(stress);
}
/// write the strain tensor in the Voigt notation.
template
static constexpr inline decltype(auto)
strainToVoigt(const Eigen::MatrixBase & strain) {
return VoigtHelper::matrixToVoigtWithFactors(strain);
}
/// write a voigt vector to stress
template
static constexpr inline void
voigtToStress(const Eigen::MatrixBase & voigt,
Eigen::MatrixBase & stress) {
VoigtHelper::voigtToMatrix(voigt, stress);
}
/// Computation of Cauchy stress tensor in the case of finite deformation from
/// the 2nd Piola-Kirchhoff for a given element type
template
void StoCauchy(ElementType el_type, GhostType ghost_type = _not_ghost);
/// Computation the Cauchy stress the 2nd Piola-Kirchhoff and the deformation
/// gradient
template
static constexpr inline void
StoCauchy(const Eigen::MatrixBase & F, const Eigen::MatrixBase & S,
Eigen::MatrixBase & sigma, const Real & C33 = 1.0);
template
static constexpr inline decltype(auto)
StoCauchy(const Eigen::MatrixBase & F, const Eigen::MatrixBase & S,
const Real & C33 = 1.0);
template
static constexpr inline void gradUToF(const Eigen::MatrixBase & grad_u,
Eigen::MatrixBase & F);
template
static constexpr inline decltype(auto)
gradUToF(const Eigen::MatrixBase & grad_u);
template
static constexpr inline void rightCauchy(const Eigen::MatrixBase & F,
Eigen::MatrixBase & C);
template
static constexpr inline decltype(auto)
rightCauchy(const Eigen::MatrixBase & F);
template
static constexpr inline void leftCauchy(const Eigen::MatrixBase & F,
Eigen::MatrixBase & B);
template
static constexpr inline decltype(auto)
leftCauchy(const Eigen::MatrixBase & F);
template
static constexpr inline void
gradUToEpsilon(const Eigen::MatrixBase & grad_u,
Eigen::MatrixBase & epsilon);
template
static constexpr inline decltype(auto)
gradUToEpsilon(const Eigen::MatrixBase & grad_u);
template
static constexpr inline void gradUToE(const Eigen::MatrixBase & grad_u,
Eigen::MatrixBase & E);
template
static constexpr inline decltype(auto)
gradUToE(const Eigen::MatrixBase & grad_u);
template
static constexpr inline void
computeDeviatoric(const Eigen::MatrixBase & sigma,
Eigen::MatrixBase & sigma_dev) {
sigma_dev =
sigma - Matrix::Identity() * sigma.trace() / dim;
}
template
static constexpr inline decltype(auto)
computeDeviatoric(const Eigen::MatrixBase & sigma) {
Matrix sigma_dev;
Material::computeDeviatoric(sigma, sigma_dev);
return sigma_dev;
}
template
static inline Real stressToVonMises(const Eigen::MatrixBase & stress);
protected:
/// converts global element to local element
inline Element convertToLocalElement(const Element & global_element) const;
/// converts local element to global element
inline Element convertToGlobalElement(const Element & local_element) const;
/// converts global quadrature point to local quadrature point
inline IntegrationPoint
convertToLocalPoint(const IntegrationPoint & global_point) const;
/// converts local quadrature point to global quadrature point
inline IntegrationPoint
convertToGlobalPoint(const IntegrationPoint & local_point) const;
/* ------------------------------------------------------------------------ */
/* DataAccessor inherited members */
/* ------------------------------------------------------------------------ */
public:
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;
template
inline void packElementDataHelper(const ElementTypeMapArray & data_to_pack,
CommunicationBuffer & buffer,
const Array & elements,
const ID & fem_id = ID()) const;
template
inline void unpackElementDataHelper(ElementTypeMapArray & data_to_unpack,
CommunicationBuffer & buffer,
const Array & elements,
const ID & fem_id = ID());
/* ------------------------------------------------------------------------ */
/* MeshEventHandler inherited members */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
void onNodesAdded(const Array &, const NewNodesEvent &) override{};
void onNodesRemoved(const Array &, const Array &,
const RemovedNodesEvent &) override{};
void onElementsAdded(const Array & element_list,
const NewElementsEvent & event) override;
void onElementsRemoved(const Array & element_list,
const ElementTypeMapArray & new_numbering,
const RemovedElementsEvent & event) override;
void onElementsChanged(const Array &, const Array &,
const ElementTypeMapArray &,
const ChangedElementsEvent &) override{};
/* ------------------------------------------------------------------------ */
/* SolidMechanicsModelEventHandler inherited members */
/* ------------------------------------------------------------------------ */
public:
virtual void beforeSolveStep();
virtual void afterSolveStep(bool converged = true);
void onDamageIteration() override;
void onDamageUpdate() override;
void onDump() override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(Name, name, const std::string &);
AKANTU_SET_MACRO(Name, name, const std::string &);
AKANTU_GET_MACRO(Model, model, const SolidMechanicsModel &)
AKANTU_GET_MACRO(ID, id, const ID &);
AKANTU_GET_MACRO(Rho, rho, Real);
AKANTU_SET_MACRO(Rho, rho, Real);
AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, Int);
/// return the potential energy for the subset of elements contained by the
/// material
Real getPotentialEnergy();
/// return the potential energy for the provided element
Real getPotentialEnergy(const Element & element);
[[deprecated("Use the interface with an Element")]] Real
getPotentialEnergy(ElementType type, Int index);
/// return the energy (identified by id) for the subset of elements contained
/// by the material
virtual Real getEnergy(const std::string & type);
/// return the energy (identified by id) for the provided element
virtual Real getEnergy(const std::string & energy_id,
const Element & element);
[[deprecated("Use the interface with an Element")]] virtual Real
getEnergy(const std::string & energy_id, ElementType type, Idx index) final {
return getEnergy(energy_id, {type, index, _not_ghost});
}
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, Idx);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GradU, gradu, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy,
Real);
AKANTU_GET_MACRO(GradU, gradu, const ElementTypeMapArray &);
AKANTU_GET_MACRO(Stress, stress, const ElementTypeMapArray &);
AKANTU_GET_MACRO(ElementFilter, element_filter,
const ElementTypeMapArray &);
AKANTU_GET_MACRO(FEEngine, fem, FEEngine &);
bool isNonLocal() const { return is_non_local; }
template
const Array & getArray(const ID & id, ElementType type,
GhostType ghost_type = _not_ghost) const;
template
Array & getArray(const ID & id, ElementType type,
GhostType ghost_type = _not_ghost);
template
const InternalField & getInternal(const ID & id) const;
template InternalField & getInternal(const ID & id);
template
inline bool isInternal(const ID & id, ElementKind element_kind) const;
template
ElementTypeMap getInternalDataPerElem(const ID & id,
ElementKind element_kind) const;
bool isFiniteDeformation() const { return finite_deformation; }
bool isInelasticDeformation() const { return inelastic_deformation; }
template inline void setParam(const ID & param, T value);
inline const Parameter & getParam(const ID & param) const;
template
void flattenInternal(const std::string & field_id,
ElementTypeMapArray & internal_flat,
GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined) const;
template
void inflateInternal(const std::string & field_id,
const ElementTypeMapArray & field,
GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined);
/// apply a constant eigengrad_u everywhere in the material
virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u,
GhostType /*ghost_type*/ = _not_ghost);
bool hasMatrixChanged(const ID & id) {
if (id == "K") {
return hasStiffnessMatrixChanged() or finite_deformation;
}
return true;
}
MatrixType getMatrixType(const ID & id) {
if (id == "K") {
return getTangentType();
}
if (id == "M") {
return _symmetric;
}
return _mt_not_defined;
}
/// specify if the matrix need to be recomputed for this material
virtual bool hasStiffnessMatrixChanged() { return true; }
/// specify the type of matrix, if not overloaded the material is not valid
/// for static or implicit computations
virtual MatrixType getTangentType() { return _mt_not_defined; }
/// static method to reteive the material factory
static MaterialFactory & getFactory();
protected:
bool isInit() const { return is_init; }
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// boolean to know if the material has been initialized
bool is_init{false};
std::map *> internal_vectors_real;
std::map *> internal_vectors_int;
std::map *> internal_vectors_bool;
protected:
ID id;
/// Link to the fem object in the model
FEEngine & fem;
/// Finite deformation
bool finite_deformation{false};
/// Finite deformation
bool inelastic_deformation{false};
/// material name
std::string name;
/// The model to witch the material belong
SolidMechanicsModel & model;
/// density : rho
Real rho{0.};
/// spatial dimension
Int spatial_dimension;
/// list of element handled by the material
ElementTypeMapArray element_filter;
/// stresses arrays ordered by element types
InternalField stress;
/// eigengrad_u arrays ordered by element types
InternalField eigengradu;
/// grad_u arrays ordered by element types
InternalField gradu;
/// Green Lagrange strain (Finite deformation)
InternalField green_strain;
/// Second Piola-Kirchhoff stress tensor arrays ordered by element types
/// (Finite deformation)
InternalField piola_kirchhoff_2;
/// potential energy by element
InternalField potential_energy;
/// tell if using in non local mode or not
bool is_non_local{false};
/// tell if the material need the previous stress state
bool use_previous_stress{false};
/// tell if the material need the previous strain state
bool use_previous_gradu{false};
/// elemental field interpolation coordinates
InternalField interpolation_inverse_coordinates;
/// elemental field interpolation points
InternalField interpolation_points_matrices;
/// vector that contains the names of all the internals that need to
/// be transferred when material interfaces move
std::vector internals_to_transfer;
private:
/// eigen_grad_u for the parser
Matrix eigen_grad_u;
};
/// standard output stream operator
inline std::ostream & operator<<(std::ostream & stream,
const Material & _this) {
_this.printself(stream);
return stream;
}
} // namespace akantu
#include "material_inline_impl.hh"
#include "internal_field_tmpl.hh"
#include "random_internal_field_tmpl.hh"
/* -------------------------------------------------------------------------- */
/* Auto loop */
/* -------------------------------------------------------------------------- */
/// This can be used to automatically write the loop on quadrature points in
/// functions such as computeStress. This macro in addition to write the loop
/// provides two tensors (matrices) sigma and grad_u
#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \
auto && grad_u_view = \
make_view(this->gradu(el_type, ghost_type), this->spatial_dimension, \
this->spatial_dimension); \
\
auto stress_view = \
make_view(this->stress(el_type, ghost_type), this->spatial_dimension, \
this->spatial_dimension); \
\
if (this->isFiniteDeformation()) { \
stress_view = make_view(this->piola_kirchhoff_2(el_type, ghost_type), \
this->spatial_dimension, this->spatial_dimension); \
} \
\
for (auto && data : zip(grad_u_view, stress_view)) { \
[[gnu::unused]] auto && grad_u = std::get<0>(data); \
[[gnu::unused]] auto && sigma = std::get<1>(data)
#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END }
/// This can be used to automatically write the loop on quadrature points in
/// functions such as computeTangentModuli. This macro in addition to write the
/// loop provides two tensors (matrices) sigma_tensor, grad_u, and a matrix
/// where the elemental tangent moduli should be stored in Voigt Notation
#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \
auto && grad_u_view = \
make_view(this->gradu(el_type, ghost_type), this->spatial_dimension, \
this->spatial_dimension); \
\
auto && stress_view = \
make_view(this->stress(el_type, ghost_type), this->spatial_dimension, \
this->spatial_dimension); \
\
auto tangent_size = \
Material::getTangentStiffnessVoigtSize(this->spatial_dimension); \
\
auto && tangent_view = make_view(tangent_mat, tangent_size, tangent_size); \
\
for (auto && data : zip(grad_u_view, stress_view, tangent_view)) { \
[[gnu::unused]] auto && grad_u = std::get<0>(data); \
[[gnu::unused]] auto && sigma = std::get<1>(data); \
auto && tangent = std::get<2>(data);
#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END }
/* -------------------------------------------------------------------------- */
namespace akantu {
namespace {
template class Mat> bool instantiateMaterial(const ID & id) {
return MaterialFactory::getInstance().registerAllocator(
id,
[](Int dim, const ID &, SolidMechanicsModel & model, const ID & id) {
return tuple_dispatch(
[&](auto && _) -> std::unique_ptr {
constexpr auto && dim_ = aka::decay_v;
return std::make_unique>(model, id);
},
dim);
});
}
} // namespace
} // namespace akantu
#endif /* AKANTU_MATERIAL_HH_ */
diff --git a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
index 158e5c66d..3166db6bd 100644
--- a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
+++ b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
@@ -1,248 +1,250 @@
/**
* Copyright (©) 2013-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_elastic_linear_anisotropic.hh"
#include "solid_mechanics_model.hh"
#include
#include
namespace akantu {
/* -------------------------------------------------------------------------- */
template
MaterialElasticLinearAnisotropic::MaterialElasticLinearAnisotropic(
SolidMechanicsModel & model, const ID & id, bool symmetric)
: Material(model, id), rot_mat(dim, dim), Cprime(dim * dim, dim * dim),
C(voigt_h::size, voigt_h::size), eigC(voigt_h::size),
symmetric(symmetric), was_stiffness_assembled(false) {
AKANTU_DEBUG_IN();
for (int i : arange(dim)) {
this->dir_vecs.emplace_back(std::make_unique>());
auto && n = *this->dir_vecs.back();
n.zero();
n[i] = 1.;
this->registerParam("n" + std::to_string(i + 1), *(this->dir_vecs.back()),
_pat_parsmod, "Direction of main material axis");
}
for (auto i : arange(voigt_h::size)) {
decltype(i) start = 0;
if (this->symmetric) {
start = i;
}
for (auto j : arange(start, voigt_h::size)) {
auto param = "C" + std::to_string(i + 1) + std::to_string(j + 1);
this->registerParam(param, this->Cprime(i, j), Real(0.), _pat_parsmod,
"Coefficient " + param);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template void MaterialElasticLinearAnisotropic::initMaterial() {
AKANTU_DEBUG_IN();
Material::initMaterial();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MaterialElasticLinearAnisotropic::updateInternalParameters() {
Material::updateInternalParameters();
if (this->symmetric) {
for (auto i : arange(voigt_h::size)) {
for (auto j : arange(i + 1, voigt_h::size)) {
this->Cprime(j, i) = this->Cprime(i, j);
}
}
}
this->rotateCprime();
this->C.eig(this->eigC);
this->was_stiffness_assembled = false;
}
/* -------------------------------------------------------------------------- */
template void MaterialElasticLinearAnisotropic::rotateCprime() {
// start by filling the empty parts fo Cprime
auto diff = Dim * Dim - voigt_h::size;
for (auto i : arange(voigt_h::size, Dim * Dim)) {
for (auto j : arange(Dim * Dim)) {
this->Cprime(i, j) = this->Cprime(i - diff, j);
}
}
for (auto i : arange(Dim * Dim)) {
for (auto j : arange(voigt_h::size, Dim * Dim)) {
this->Cprime(i, j) = this->Cprime(i, j - diff);
}
}
// construction of rotator tensor
// normalise rotation matrix
for (auto j : arange(Dim)) {
auto && rot_vec = this->rot_mat(j);
rot_vec = *this->dir_vecs[j];
rot_vec.normalize();
}
// make sure the vectors form a right-handed base
Real test_axis = 0.;
if (Dim == 2) {
Vector v1 = Vector::Zero();
Vector v2 = Vector::Zero();
v1.block(0, 0) = this->rot_mat(0);
v2.block(0, 0) = this->rot_mat(1);
Vector v3 = v1.cross(v2);
if (v3.norm() < 8 * std::numeric_limits::epsilon()) {
AKANTU_ERROR("The axis vectors parallel.");
}
v3.normalize();
test_axis = (v1.cross(v2) - v3).norm();
} else if (Dim == 3) {
Vector v1 = this->rot_mat(0);
Vector v2 = this->rot_mat(1);
Vector v3 = this->rot_mat(2);
test_axis = (v1.cross(v2) - v3).norm();
}
if (test_axis > 8 * std::numeric_limits::epsilon()) {
AKANTU_ERROR("The axis vectors do not form a right-handed coordinate "
<< "system. I. e., ||n1 x n2 - n3|| should be zero, but "
<< "it is " << test_axis << ".");
}
// create the rotator and the reverse rotator
Matrix rotator;
Matrix revrotor;
for (auto i : arange(Dim)) {
for (auto j : arange(Dim)) {
for (auto k : arange(Dim)) {
for (auto l : arange(Dim)) {
auto I = voigt_h::mat[i][j];
auto J = voigt_h::mat[k][l];
rotator(I, J) = this->rot_mat(k, i) * this->rot_mat(l, j);
revrotor(I, J) = this->rot_mat(i, k) * this->rot_mat(j, l);
}
}
}
}
// create the full rotated matrix
Matrix Cfull(Dim * Dim, Dim * Dim);
Cfull = rotator * Cprime * revrotor;
for (auto i : arange(voigt_h::size)) {
for (auto j : arange(voigt_h::size)) {
this->C(i, j) = Cfull(i, j);
}
}
}
/* -------------------------------------------------------------------------- */
template
void MaterialElasticLinearAnisotropic::computeStress(
ElementType el_type, GhostType ghost_type) {
auto && arguments = getArguments(el_type, ghost_type);
for (auto && data : arguments) {
this->computeStressOnQuad(data);
}
}
/* -------------------------------------------------------------------------- */
template
void MaterialElasticLinearAnisotropic::computeTangentModuli(
ElementType el_type, Array & tangent_matrix, GhostType ghost_type) {
auto && arguments =
Material::getArgumentsTangent(tangent_matrix, el_type, ghost_type);
for (auto && args : arguments) {
this->computeTangentModuliOnQuad(args);
}
this->was_stiffness_assembled = true;
}
/* -------------------------------------------------------------------------- */
template
void MaterialElasticLinearAnisotropic::computePotentialEnergy(
ElementType el_type) {
AKANTU_DEBUG_ASSERT(!this->finite_deformation,
"finite deformation not possible in material anisotropic "
"(TO BE IMPLEMENTED)");
auto && arguments = Material::getArguments(el_type, _not_ghost);
for (auto && args :
zip(arguments, this->potential_energy(el_type, _not_ghost))) {
this->computePotentialEnergyOnQuad(std::get<0>(args), std::get<1>(args));
}
}
/* -------------------------------------------------------------------------- */
template
void MaterialElasticLinearAnisotropic::computePotentialEnergyByElement(
- ElementType type, Int index, Vector & epot_on_quad_points) {
+ const Element & element, Vector & epot_on_quad_points) {
+ auto type = element.type;
+ auto index = element.element;
auto gradu_view = make_view(this->gradu(type));
auto stress_view = make_view(this->stress(type));
auto nb_quadrature_points = this->fem.getNbIntegrationPoints(type);
auto gradu_it = gradu_view.begin() + index * nb_quadrature_points;
auto gradu_end = gradu_it + nb_quadrature_points;
auto stress_it = stress_view.begin() + index * nb_quadrature_points;
auto stress_end = stress_it + nb_quadrature_points;
auto epot_quad = epot_on_quad_points.begin();
Matrix grad_u(dim, dim);
for (auto data : zip("grad_u"_n = range(gradu_it, gradu_end),
"sigma"_n = range(stress_it, stress_end),
"Epot"_n = epot_on_quad_points)) {
this->computePotentialEnergyOnQuad(data, data["Epot"_n]);
}
}
/* -------------------------------------------------------------------------- */
template
Real MaterialElasticLinearAnisotropic::getCelerity(
const Element & /*element*/) const {
return std::sqrt(this->eigC(0) / rho);
}
/* -------------------------------------------------------------------------- */
template class MaterialElasticLinearAnisotropic<1>;
template class MaterialElasticLinearAnisotropic<2>;
template class MaterialElasticLinearAnisotropic<3>;
static bool material_is_alocated_elastic =
instantiateMaterial(
"elastic_anisotropic");
} // namespace akantu
diff --git a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh
index e19541ca9..03383fbec 100644
--- a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh
+++ b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh
@@ -1,133 +1,133 @@
/**
* 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 "aka_common.hh"
#include "material.hh"
#include "material_elastic.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_MATERIAL_ELASTIC_LINEAR_ANISOTROPIC_HH_
#define AKANTU_MATERIAL_ELASTIC_LINEAR_ANISOTROPIC_HH_
namespace akantu {
/**
* General linear anisotropic elastic material
* The only constraint on the elastic tensor is that it can be represented
* as a symmetric 6x6 matrix (3D) or 3x3 (2D).
*
* parameters in the material files :
* - rho : density (default: 0)
* - C_ij : entry on the stiffness
*/
template class MaterialElasticLinearAnisotropic : public Material {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MaterialElasticLinearAnisotropic(SolidMechanicsModel & model,
const ID & id = "", bool symmetric = true);
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
void initMaterial() override;
/// constitutive law for all element of a type
void computeStress(ElementType el_type,
GhostType ghost_type = _not_ghost) override;
/// compute the tangent stiffness matrix for an element type
void computeTangentModuli(ElementType el_type, Array & tangent_matrix,
GhostType ghost_type = _not_ghost) override;
/// compute the elastic potential energy
void computePotentialEnergy(ElementType el_type) override;
void updateInternalParameters() override;
bool hasStiffnessMatrixChanged() override {
return (not was_stiffness_assembled);
}
MatrixType getTangentType() override { return _symmetric; }
protected:
// compute C from Cprime
void rotateCprime();
/// constitutive law for a given quadrature point
template inline void computeStressOnQuad(Args && args) const;
/// tangent matrix for a given quadrature point
template
inline void computeTangentModuliOnQuad(Args && args) const;
template
inline void computePotentialEnergyOnQuad(Args && args, Real & epot);
void
- computePotentialEnergyByElement(ElementType type, Int index,
+ computePotentialEnergyByElement(const Element & element,
Vector & epot_on_quad_points) override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
decltype(auto) getArguments(ElementType el_type, GhostType ghost_type) {
return Material::template getArguments(el_type, ghost_type);
}
/// compute max wave celerity
Real getCelerity(const Element & element) const override;
AKANTU_GET_MACRO(VoigtStiffness, C, Matrix);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
using voigt_h = VoigtHelper;
/// direction matrix and vectors
std::vector>> dir_vecs;
Matrix rot_mat;
/// Elastic stiffness tensor in material frame and full vectorised notation
Matrix Cprime;
/// Elastic stiffness tensor in voigt notation
Matrix C;
/// eigenvalues of stiffness tensor
Vector eigC;
bool symmetric;
/// defines if the stiffness was computed
bool was_stiffness_assembled;
};
} // namespace akantu
#include "material_elastic_linear_anisotropic_inline_impl.hh"
#endif /* AKANTU_MATERIAL_ELASTIC_LINEAR_ANISOTROPIC_HH_ */