SaintVenantKirchhoffMaterial.h

#pragma once
 
#include <SofaCaribou/Material/HyperelasticMaterial.h>
 
namespace SofaCaribou::material {
 
template<class DataTypes>
class SaintVenantKirchhoffMaterial : public HyperelasticMaterial<DataTypes> {
    static constexpr auto Dimension = DataTypes::spatial_dimensions;
    using Coord = typename DataTypes::Coord;
    using Real  = typename Coord::value_type;
public:
    SOFA_CLASS(SOFA_TEMPLATE(SaintVenantKirchhoffMaterial, DataTypes), SOFA_TEMPLATE(HyperelasticMaterial, DataTypes));
 
    SaintVenantKirchhoffMaterial()
        : d_young_modulus(initData(&d_young_modulus,
            Real(1000), "young_modulus",
            "Young's modulus of the material",
            true /*displayed_in_GUI*/, false /*read_only_in_GUI*/))
        , d_poisson_ratio(initData(&d_poisson_ratio,
            Real(0.3),  "poisson_ratio",
           "Poisson's ratio of the material",
           true /*displayed_in_GUI*/, false /*read_only_in_GUI*/))
    {
    }
 
    void before_update() override {
        const Real young_modulus = d_young_modulus.getValue();
        const Real poisson_ratio = d_poisson_ratio.getValue();
        mu = young_modulus / (2.0 * (1.0 + poisson_ratio));
        l = young_modulus * poisson_ratio / ((1.0 + poisson_ratio) * (1.0 - 2.0 * poisson_ratio));
 
        Eigen::Matrix<Real, 2*Dimension, 2*Dimension> D;
        D <<
              l + 2*mu,    l,          l,       0,  0,  0,
                l,       l + 2*mu,     l,       0,  0,  0,
                l,         l,        l + 2*mu,  0,  0,  0,
                0,         0,          0,      mu,  0,  0,
                0,         0,          0,       0, mu,  0,
                0,         0,          0,       0,  0, mu;
        C = D;
    }
 
    Real
    strain_energy_density(const Real & /*J*/, const Eigen::Matrix<Real, Dimension, Dimension>  & C) const override {
        static const auto Id = Eigen::Matrix<Real, Dimension, Dimension, Eigen::RowMajor>::Identity();
        const auto E = (1/2. * (C - Id)).eval();
        const auto trE  = E.trace();
        const auto trEE = (E*E).trace();
        return l/2.*(trE*trE) + mu*trEE;
    }
 
    Eigen::Matrix<Real, Dimension, Dimension>
    PK2_stress(const Real & /*J*/, const Eigen::Matrix<Real, Dimension, Dimension>  & C) const override {
        static const auto Id = Eigen::Matrix<Real, Dimension, Dimension, Eigen::RowMajor>::Identity();
        const auto E = (1/2. * (C - Id)).eval();
        return l*E.trace()*Id + 2*mu*E;
    }
 
    Eigen::Matrix<Real, 6, 6>
    PK2_stress_jacobian(const Real & /*J*/, const Eigen::Matrix<Real, Dimension, Dimension> & /*C*/) const override {
        return C;
    }
 
private:
    // Private members
    Real mu; // Lame's mu parameter
    Real l;  // Lame's lambda parameter
    Eigen::Matrix<Real, 2*Dimension, 2*Dimension> C; // Constant elasticity matrix (jacobian of the stress tensor)
 
    // Data members
    sofa::core::objectmodel::Data<Real> d_young_modulus;
    sofa::core::objectmodel::Data<Real> d_poisson_ratio;
};
 
} // namespace SofaCaribou::material