NeoHookeanMaterial.h

#pragma once
 
#include <Caribou/Algebra/Tensor.h>
#include <SofaCaribou/Material/HyperelasticMaterial.h>
 
namespace SofaCaribou::material {
 
template<class DataTypes>
class NeoHookeanMaterial : 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(NeoHookeanMaterial, DataTypes), SOFA_TEMPLATE(HyperelasticMaterial, DataTypes));
 
    NeoHookeanMaterial()
        : 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));
    }
 
    Real
    strain_energy_density(const Real & J, const Eigen::Matrix<Real, Dimension, Dimension>  & C) const override {
        const auto lnJ = log(J);
        return mu/2.*(C.trace()-3) - mu*lnJ + l/2 *lnJ*lnJ;
    }
 
    virtual 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 Ci = C.inverse();
 
        return (Id - Ci)*mu + Ci*(l*log(J));
    }
 
    virtual Eigen::Matrix<Real, 6, 6>
    PK2_stress_jacobian(const Real & J, const Eigen::Matrix<Real, Dimension, Dimension> & C) const override {
        using caribou::algebra::symmetric_dyad_1;
        using caribou::algebra::symmetric_dyad_2;
 
        const auto Ci = C.inverse().eval();
 
        Eigen::Matrix<Real, 6, 6> D = l*symmetric_dyad_1(Ci) + 2*(mu - l*log(J))*symmetric_dyad_2(Ci);
        return D;
    }
 
private:
    // Private members
    Real mu; // Lame's mu parameter
    Real l;  // Lame's lambda parameter
 
    // Data members
    sofa::core::objectmodel::Data<Real> d_young_modulus;
    sofa::core::objectmodel::Data<Real> d_poisson_ratio;
};
 
} // namespace SofaCaribou::material