# Writing a New Constitutive Law¶

The abstraction for a constitutive law in µSpectre** is the Material, and new such materials must inherit from the class muSpectre::MaterialBase. Most often, however, it will be most convenient to inherit from the derived class muSpectre::MaterialMuSpectre, as it implements a lot of the machinery that is commonly used by constitutive laws. This section describes how to implement a new constitutive law with internal variables (sometimes also called state variables). The example material implemented here is MaterialTutorial, an objective linear elastic law with a distribution of eigenstrains as internal variables. The constitutive law is defined by the relationship between material (or second Piola-Kirchhoff) stress $$\mathbf{S}$$ and Green-Lagrange strain $$\mathbf{E}$$

\begin{align} \mathbf{S} &= \mathbb C:\left(\mathbf{E}-\mathbf{e}\right), \\ S_{ij} &= C_{ijkl}\left(E_{kl}-e_{kl}\right), \\ \end{align}

where $$\mathbb C$$ is the elastic stiffness tensor and $$\mathbf e$$ is the local eigenstrain. Note that the implementation sketched out here is the most convenient to quickly get started with using µSpectre**, but not the most efficient one. For a most efficient implementation, refer to the implementation of muSpectre::MaterialLinearElastic2.

## The muSpectre::MaterialMuSpectre class¶

The class muSpectre::MaterialMuSpectre is defined in material_muSpectre_base.hh and takes three template parameters;

1. class Material is a CRTP and names the material inheriting from it. The reason for this construction is that we want to avoid virtual method calls from muSpectre::MaterialMuSpectre to its derived classes. Rather, we want muSpectre::MaterialMuSpectre to be able to call methods of the inheriting class directly without runtime overhead.

2. Dim_t DimS defines the number of spatial dimensions of the problem, i.e., whether we are dealing with a two- or three-dimensional grid of pixels/voxels.

3. Dim_t DimM defines the number of dimensions of our material description. This value will typically be the same as DimS, but in cases like generalised plane strain, we can for instance have a three three-dimensional material response in a two-dimensional pixel grid.

The main job of muSpectre::MaterialMuSpectre is to

1. loop over all pixels to which this material has been assigned, transform the global gradient $$\mathbf{F}$$ (or small strain tensor $$\boldsymbol\varepsilon$$) into the new material’s required strain measure (e.g., the Green-Lagrange strain tensor $$\mathbf{E}$$),

2. for each pixel evaluate the constitutive law by calling its evaluate_stress (computes the stress response) or evaluate_stress_tangent (computes both stress and consistent tangent) method with the local strain and internal variables, and finally

3. transform the stress (and possibly tangent) response from the material’s stress measure into first Piola-Kirchhoff stress $$\mathbf{P}$$ (or Cauchy stress $$\boldsymbol\sigma$$ in small strain).

In order to use these facilities, the new material needs to inherit from muSpectre::MaterialMuSpectre (where we calculation of the response) and specialise the type muSpectre::MaterialMuSpectre_traits (where we tell muSpectre::MaterialMuSpectre how to use the new material). These two steps are described here for our example material.

Specialising the muSpectre::MaterialMuSpectre_traits structure *********************************************************************** This structure is templated by the new material (in this case MaterialTutorial) and needs to specify

1. the types used to communicate per-pixel strains, stresses and stiffness tensors to the material (i.e., whether you want to get maps to Eigen matrices or raw pointers, or …). Here we will use the convenient muSpectre::MatrixFieldMap for strains and stresses, and muSpectre::T4MatrixFieldMap for the stiffness. Look through the classes deriving from muSpectre::FieldMap for all available options.

2. the strain measure that is expected (e.g., gradient, Green-Lagrange strain, left Cauchy-Green strain, etc.). Here we will use Green-Lagrange strain. The options are defined by the enum muSpectre::StrainMeasure.

3. the stress measure that is computed by the law (e.g., Cauchy, first Piola-Kirchhoff, etc,). Here, it will be first Piola-Kirchhoff stress. The available options are defined by the enum muSpectre::StressMeasure.

Our traits look like this (assuming we are in the namespace muSpectre:

template <Dim_t DimS, Dim_t DimM>
struct MaterialMuSpectre_traits<MaterialTutorial<DimS, DimM>>
{
//! global field collection
using GFieldCollection_t = typename
GlobalFieldCollection<DimS, DimM>;

//! expected map type for strain fields
using StrainMap_t = MatrixFieldMap<GFieldCollection_t, Real, DimM, DimM, true>;
//! expected map type for stress fields
using StressMap_t = MatrixFieldMap<GFieldCollection_t, Real, DimM, DimM>;
//! expected map type for tangent stiffness fields
using TangentMap_t = T4MatrixFieldMap<GFieldCollection_t, Real, DimM>;

//! declare what type of strain measure your law takes as input
constexpr static auto strain_measure{StrainMeasure::GreenLagrange};
//! declare what type of stress measure your law yields as output
constexpr static auto stress_measure{StressMeasure::PK2};

//! local field_collections used for internals
using LFieldColl_t = LocalFieldCollection<DimS, DimM>;
//! local strain type
using LStrainMap_t = MatrixFieldMap<LFieldColl_t, Real, DimM, DimM, true>;
//! elasticity with eigenstrain
using InternalVariables = std::tuple<LStrainMap_t>;

};


## Implementing the new material¶

The new law needs to implement the methods add_pixel, get_internals, evaluate_stress, and evaluate_stress_tangent. Below is a commented example header:

template <Dim_t DimS, Dim_t DimM>
class MaterialTutorial:
public MaterialMuSpectre<MaterialTutorial<DimS, DimM>, DimS, DimM>
{
public:
//! traits of this material
using traits = MaterialMuSpectre_traits<MaterialTutorial>;

//! Type of container used for storing eigenstrain
using InternalVariables = typename traits::InternalVariables;

//! Construct by name, Young's modulus and Poisson's ratio
MaterialTutorial(std::string name, Real young, Real poisson);

/**
* evaluates second Piola-Kirchhoff stress given the Green-Lagrange
* strain (or Cauchy stress if called with a small strain tensor)
*/
template <class s_t, class eigen_s_t>
inline decltype(auto) evaluate_stress(s_t && E, eigen_s_t && E_eig);

/**
* evaluates both second Piola-Kirchhoff stress and stiffness given
* the Green-Lagrange strain (or Cauchy stress and stiffness if
* called with a small strain tensor)
*/
template <class s_t, class eigen_s_t>
inline decltype(auto)
evaluate_stress_tangent(s_t &&  E, eigen_s_t && E_eig);

/**
* return the internals tuple (needed by muSpectre::MaterialMuSpectre)
*/
InternalVariables & get_internals() {
return this->internal_variables;};

/**
*/
const Eigen::Matrix<Real, DimM, DimM> & E_eig);

protected:
//! stiffness tensor
T4Mat<Real, DimM> C;
//! storage for eigenstrain
using Field_t =
TensorField<LocalFieldCollection<DimS,DimM>, Real, secondOrder, DimM>;
Field_t & eigen_field; //!< field of eigenstrains
//! tuple for iterable eigen_field
InternalVariables internal_variables;
private:
};


A possible implementation for the constructor would be:

template <Dim_t DimS, Dim_t DimM>
MaterialTutorial<DimS, DimM>::MaterialTutorial(std::string name,
Real young,
Real poisson)
:MaterialMuSpectre<MaterialTutorial, DimS, DimM>(name) {

// Lamé parameters
Real lambda{young*poisson/((1+poisson)*(1-2*poisson))};
Real mu{young/(2*(1+poisson))};

// Kronecker delta
Eigen::Matrix<Real, DimM, DimM> del{Eigen::Matrix<Real, DimM, DimM>::Identity()};

// fill the stiffness tensor
this->C.setZero();
for (Dim_t i = 0; i < DimM; ++i) {
for (Dim_t j = 0; j < DimM; ++j) {
for (Dim_t k = 0; k < DimM; ++k) {
for (Dim_t l = 0; l < DimM; ++l) {
get(this->C, i, j, k, l) += (lambda * del(i,j)*del(k,l) +
mu * (del(i,k)*del(j,l) + del(i,l)*del(j,k)));
}
}
}
}
}


as an exercise, you could check how muSpectre::MaterialLinearElastic1 uses µSpectre**’s materials toolbox (in namespace MatTB) to compute $$\mathbb C$$ in a much more convenient fashion. The evaluation of the stress could be (here, we make use of the Matrices namespace that defines common tensor algebra operations):

template <Dim_t DimS, Dim_t DimM>
template <class s_t, class eigen_s_t>
decltype(auto)
MaterialTutorial<DimS, DimM>::
evaluate_stress(s_t && E, eigen_s_t && E_eig) {
return Matrices::tens_mult(this->C, E-E_eig);
}


The remaining two methods are straight-forward:

template <Dim_t DimS, Dim_t DimM>
template <class s_t, class eigen_s_t>
decltype(auto)
MaterialTutorial<DimS, DimM>::
evaluate_stress_tangent(s_t && E, eigen_s_t && E_eig) {
return return std::make_tuple
(evaluate_stress(E, E_eig),
this->C);
}

template <Dim_t DimS, Dim_t DimM>
InternalVariables &
MaterialTutorial<DimS, DimM>::get_internals() {
return this->internal_variables;
}


Note that the methods evaluate_stress and evaluate_stress_tangent need to be in the header, as both their input parameter types and output type depend on the compile-time context.