FemElastoDynamics.h

Go to the documentation of this file.

 
#ifndef PBAT_SIM_DYNAMICS_FEMELASTODYNAMICS_H
#define PBAT_SIM_DYNAMICS_FEMELASTODYNAMICS_H
 
#include "pbat/common/Concepts.h"
#include "pbat/fem/HyperElasticPotential.h"
#include "pbat/fem/Jacobian.h"
#include "pbat/fem/Mesh.h"
#include "pbat/fem/MeshQuadrature.h"
#include "pbat/fem/ShapeFunctions.h"
#include "pbat/io/Archive.h"
#include "pbat/physics/HyperElasticity.h"
#include "pbat/sim/integration/Bdf.h"
 
#include <algorithm>
#include <cassert>
#include <exception>
#include <type_traits>
 
namespace pbat::sim::dynamics {

template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar = Scalar,
    common::CIndex TIndex          = Index>
struct FemElastoDynamics
{
    using MeshType              = fem::Mesh<TElement, Dims, TScalar, TIndex>; 
    using ElementType           = TElement;                                   
    using ElasticEnergyType     = THyperElasticEnergy;               
    using BdfType               = integration::Bdf<TScalar, TIndex>; 
    using ScalarType            = TScalar;                           
    using IndexType             = TIndex;                            
    static auto constexpr kDims = Dims;                              
 
    MeshType mesh; 
    BdfType bdf;   
 
    Eigen::Matrix<ScalarType, kDims, Eigen::Dynamic>
        fext; 
    Eigen::Vector<ScalarType, Eigen::Dynamic> m; 
    Eigen::Matrix<ScalarType, kDims, Eigen::Dynamic>
        xtilde;                                         
    Eigen::Matrix<ScalarType, kDims, Eigen::Dynamic> x; 
    Eigen::Matrix<ScalarType, kDims, Eigen::Dynamic> v; 
 
    Eigen::Vector<IndexType, Eigen::Dynamic> egU; 
    Eigen::Vector<ScalarType, Eigen::Dynamic>
        wgU; 
    Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic>
        GNegU; 
    Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic>
        lamegU; 
    Eigen::Vector<ScalarType, Eigen::Dynamic>
        UgU; 
    Eigen::Matrix<ScalarType, kDims * ElementType::kNodes, Eigen::Dynamic>
        GgU; 
    Eigen::Matrix<ScalarType, kDims * ElementType::kNodes, Eigen::Dynamic>
        HgU; 
 
    IndexType ndbc; 
    Eigen::Vector<IndexType, Eigen::Dynamic>
        dbc; 
    Eigen::Vector<bool, Eigen::Dynamic> dmask; 

    FemElastoDynamics() = default;
    FemElastoDynamics(
        Eigen::Ref<Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic> const> const& V,
        Eigen::Ref<Eigen::Matrix<IndexType, Eigen::Dynamic, Eigen::Dynamic> const> const& C);
    void Construct(
        Eigen::Ref<Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic> const> const& V,
        Eigen::Ref<Eigen::Matrix<IndexType, Eigen::Dynamic, Eigen::Dynamic> const> const& C);
    template <class TDerivedX0, class TDerivedV0>
    void SetInitialConditions(
        Eigen::DenseBase<TDerivedX0> const& x0,
        Eigen::DenseBase<TDerivedV0> const& v0);
    void SetMassMatrix(ScalarType rho);
    void SetElasticEnergy(ScalarType mu, ScalarType lambda);
    void SetExternalLoad(Eigen::Vector<ScalarType, kDims> const& b);
    void SetTimeIntegrationScheme(ScalarType dt = ScalarType(1e-2), int s = 1);
    template <typename TDerivedDirichletMask>
    void Constrain(Eigen::DenseBase<TDerivedDirichletMask> const& D);
    template <class TDerivedEg, class TDerivedWg, class TDerivedXg, class TDerivedRhog>
    void SetMassMatrix(
        Eigen::DenseBase<TDerivedEg> const& eg,
        Eigen::MatrixBase<TDerivedWg> const& wg,
        Eigen::MatrixBase<TDerivedXg> const& Xg,
        Eigen::MatrixBase<TDerivedRhog> const& rhog);
    template <
        class TDerivedEg,
        class TDerivedWg,
        class TDerivedXg,
        class TDerivedMug,
        class TDerivedLambdag>
    void SetElasticEnergy(
        Eigen::DenseBase<TDerivedEg> const& eg,
        Eigen::DenseBase<TDerivedWg> const& wg,
        Eigen::MatrixBase<TDerivedXg> const& Xg,
        Eigen::DenseBase<TDerivedMug> const& mug,
        Eigen::DenseBase<TDerivedLambdag> const& lambdag);
    template <class TDerivedEg, class TDerivedWg, class TDerivedXg, class TDerivedBg>
    void SetExternalLoad(
        Eigen::DenseBase<TDerivedEg> const& eg,
        Eigen::MatrixBase<TDerivedWg> const& wg,
        Eigen::MatrixBase<TDerivedXg> const& Xg,
        Eigen::MatrixBase<TDerivedBg> const& bg);
    void SetupTimeIntegrationOptimization();
    void ComputeElasticEnergy(
        int eElasticComputationFlags,
        fem::EHyperElasticSpdCorrection eSpdCorrectionFlags);
    auto M() const { return m.replicate(1, kDims).transpose().reshaped(); };
    auto aext() const { return (fext * m.cwiseInverse().asDiagonal()); }
    bool IsDirichletNode(IndexType node) const { return dmask(node); }
    bool IsDirichletDof(IndexType i) const { return dmask(i / kDims); }
    auto DirichletNodes() const { return dbc.tail(ndbc); }
    auto DirichletCoordinates() const { return x(Eigen::placeholders::all, DirichletNodes()); }
    auto DirichletVelocities() const { return v(Eigen::placeholders::all, DirichletNodes()); }
    auto DirichletCoordinates() { return x(Eigen::placeholders::all, DirichletNodes()); }
    auto DirichletVelocities() { return v(Eigen::placeholders::all, DirichletNodes()); }
    auto FreeNodes() const { return dbc.head(dbc.size() - ndbc); }
    auto FreeCoordinates() const { return x(Eigen::placeholders::all, FreeNodes()); }
    auto FreeVelocities() const { return v(Eigen::placeholders::all, FreeNodes()); }
    auto FreeCoordinates() { return x(Eigen::placeholders::all, FreeNodes()); }
    auto FreeVelocities() { return v(Eigen::placeholders::all, FreeNodes()); }
    void Serialize(io::Archive& archive) const;
    void Deserialize(io::Archive const& archive);
};
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
inline FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::FemElastoDynamics(
    Eigen::Ref<Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic> const> const& V,
    Eigen::Ref<Eigen::Matrix<IndexType, Eigen::Dynamic, Eigen::Dynamic> const> const& C)
{
    Construct(V, C);
}
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
inline void FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::Construct(
    Eigen::Ref<Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic> const> const& V,
    Eigen::Ref<Eigen::Matrix<IndexType, Eigen::Dynamic, Eigen::Dynamic> const> const& C)
{
    mesh.Construct(V, C);
    auto const nNodes = mesh.X.cols();
    SetInitialConditions(
        mesh.X,
        Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic>::Zero(kDims, nNodes));
    // Mass
    ScalarType constexpr rho{1e3};
    SetMassMatrix(rho);
    // Elasticity
    ScalarType constexpr Y  = 1e6;
    ScalarType constexpr nu = 0.45;
    auto const [mu, lambda] = physics::LameCoefficients(Y, nu);
    SetElasticEnergy(mu, lambda);
    // External load
    Eigen::Vector<ScalarType, kDims> load = Eigen::Vector<ScalarType, kDims>::Zero();
    load(load.rows() - 1)                 = rho * ScalarType(-9.81);
    SetExternalLoad(load);
    // Time integration scheme
    ScalarType constexpr dt{1e-2};
    int constexpr bdfstep = 1;
    SetTimeIntegrationScheme(dt, bdfstep);
    // Unconstrained
    Constrain(Eigen::Vector<bool, Eigen::Dynamic>::Constant(nNodes, false));
}
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
template <class TDerivedX0, class TDerivedV0>
inline void
FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::SetInitialConditions(
    Eigen::DenseBase<TDerivedX0> const& x0,
    Eigen::DenseBase<TDerivedV0> const& v0)
{
    x = x0.reshaped(kDims, x0.size() / kDims);
    v = v0.reshaped(kDims, v0.size() / kDims);
    bdf.SetOrder(2);
    bdf.SetInitialConditions(x0.reshaped(), v0.reshaped());
}
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
inline void FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::SetMassMatrix(
    ScalarType rho)
{
    auto constexpr kOrder     = 2 * TElement::kOrder;
    IndexType const nElements = static_cast<IndexType>(mesh.E.cols());
    Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic> const XgM =
        mesh.template QuadraturePoints<kOrder>();
    auto const nQuadPtsPerElementM = XgM.cols() / nElements;
    // Mass
    SetMassMatrix(
        Eigen::Vector<IndexType, Eigen::Dynamic>::LinSpaced(nElements, IndexType(0), nElements - 1)
            .replicate(1, nQuadPtsPerElementM)
            .transpose()
            .reshaped() /*eg*/,
        fem::MeshQuadratureWeights<kOrder>(mesh).reshaped() /*wg*/,
        XgM /*Xg*/,
        Eigen::Vector<ScalarType, Eigen::Dynamic>::Constant(XgM.cols(), rho) /*rhog*/
    );
}
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
inline void
FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::SetElasticEnergy(
    ScalarType mu,
    ScalarType lambda)
{
    auto constexpr kOrder     = TElement::kOrder;
    IndexType const nElements = static_cast<IndexType>(mesh.E.cols());
    // Compute mesh quadrature points
    Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic> const XgU =
        mesh.template QuadraturePoints<kOrder>();
    auto const nQuadPtsPerElementU = XgU.cols() / nElements;
    // Elasticity
    SetElasticEnergy(
        Eigen::Vector<IndexType, Eigen::Dynamic>::LinSpaced(nElements, IndexType(0), nElements - 1)
            .replicate(1, nQuadPtsPerElementU)
            .transpose()
            .reshaped() /*eg*/,
        fem::MeshQuadratureWeights<kOrder>(mesh).reshaped() /*wg*/,
        XgU /*Xg*/,
        Eigen::Vector<ScalarType, Eigen::Dynamic>::Constant(XgU.cols(), mu) /*mug*/,
        Eigen::Vector<ScalarType, Eigen::Dynamic>::Constant(XgU.cols(), lambda) /*lambdag*/);
}
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
inline void
FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::SetExternalLoad(
    Eigen::Vector<ScalarType, kDims> const& b)
{
    auto constexpr kOrder     = TElement::kOrder;
    IndexType const nElements = static_cast<IndexType>(mesh.E.cols());
    Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic> const XgB =
        mesh.template QuadraturePoints<kOrder>();
    auto const nQuadPtsPerElementB = XgB.cols() / nElements;
    // External load
    SetExternalLoad(
        Eigen::Vector<IndexType, Eigen::Dynamic>::LinSpaced(nElements, IndexType(0), nElements - 1)
            .replicate(1, nQuadPtsPerElementB)
            .transpose()
            .reshaped() /*eg*/,
        fem::MeshQuadratureWeights<kOrder>(mesh).reshaped() /*wg*/,
        XgB /*Xg*/,
        b.replicate(1, XgB.cols()) /*bg*/
    );
}
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
inline void
FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::SetTimeIntegrationScheme(
    ScalarType dt,
    int s)
{
    bdf.SetStep(s);
    bdf.SetTimeStep(dt);
}
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
template <typename TDerivedDirichletMask>
inline void FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::Constrain(
    Eigen::DenseBase<TDerivedDirichletMask> const& D)
{
    static_assert(
        std::is_same_v<typename TDerivedDirichletMask::Scalar, bool>,
        "Dirichlet mask must be of type bool");
    IndexType const nNodes = static_cast<IndexType>(mesh.X.cols());
    assert(D.size() == nNodes);
    dmask = D.template cast<bool>();
    dbc.setLinSpaced(nNodes, IndexType(0), nNodes - 1);
    auto it = std::stable_partition(dbc.begin(), dbc.end(), [&D](IndexType i) { return not D[i]; });
    ndbc    = nNodes - std::distance(dbc.begin(), it);
}
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
template <class TDerivedEg, class TDerivedWg, class TDerivedXg, class TDerivedRhog>
inline void FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::SetMassMatrix(
    Eigen::DenseBase<TDerivedEg> const& eg,
    Eigen::MatrixBase<TDerivedWg> const& wg,
    Eigen::MatrixBase<TDerivedXg> const& Xg,
    Eigen::MatrixBase<TDerivedRhog> const& rhog)
{
    auto N = fem::ShapeFunctionMatrixAt(mesh, eg.derived(), Xg.derived());
    Eigen::SparseMatrix<ScalarType, Eigen::RowMajor, IndexType> rhogwgN =
        rhog.cwiseProduct(wg).asDiagonal() * N;
    Eigen::SparseMatrix<ScalarType, Eigen::ColMajor, IndexType> M = N.transpose() * rhogwgN;
    m.resize(M.cols());
    for (auto j = 0; j < M.cols(); ++j)
        m(j) = M.col(j).sum();
}
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
template <
    class TDerivedEg,
    class TDerivedWg,
    class TDerivedXg,
    class TDerivedMug,
    class TDerivedLambdag>
inline void
FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::SetElasticEnergy(
    Eigen::DenseBase<TDerivedEg> const& eg,
    Eigen::DenseBase<TDerivedWg> const& wg,
    Eigen::MatrixBase<TDerivedXg> const& Xg,
    Eigen::DenseBase<TDerivedMug> const& mug,
    Eigen::DenseBase<TDerivedLambdag> const& lambdag)
{
    egU   = eg;
    wgU   = wg;
    GNegU = fem::ShapeFunctionGradientsAt(mesh, eg, Xg);
    lamegU.resize(2, eg.size());
    lamegU.row(0) = mug;
    lamegU.row(1) = lambdag;
    UgU.resize(eg.size());
    GgU.resize(kDims * ElementType::kNodes, eg.size());
    HgU.resize(kDims * ElementType::kNodes, kDims * ElementType::kNodes * eg.size());
}
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
template <class TDerivedEg, class TDerivedWg, class TDerivedXg, class TDerivedBg>
inline void
FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::SetExternalLoad(
    Eigen::DenseBase<TDerivedEg> const& eg,
    Eigen::MatrixBase<TDerivedWg> const& wg,
    Eigen::MatrixBase<TDerivedXg> const& Xg,
    Eigen::MatrixBase<TDerivedBg> const& bg)
{
    Eigen::SparseMatrix<ScalarType, Eigen::RowMajor, IndexType> N =
        fem::ShapeFunctionMatrixAt(mesh, eg.derived(), Xg.derived());
    fext = bg * wg.asDiagonal() * N;
}
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
inline void FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::
    SetupTimeIntegrationOptimization()
{
    auto xtildeBdf = bdf.Inertia(0);
    auto vtildeBdf = bdf.Inertia(1);
    auto betaTilde = bdf.BetaTilde();
    xtilde.resize(kDims, xtildeBdf.size() / kDims);
    xtilde.reshaped() =
        -(xtildeBdf + betaTilde * vtildeBdf) + (betaTilde * betaTilde) * (aext().reshaped());
}
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
inline void
FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::ComputeElasticEnergy(
    int eElasticComputationFlags,
    fem::EHyperElasticSpdCorrection eSpdCorrectionFlags)
{
    fem::ToElementElasticity<ElasticEnergyType>(
        mesh,
        egU,
        wgU,
        GNegU,
        lamegU.row(0),
        lamegU.row(1),
        x.reshaped(),
        UgU,
        GgU,
        HgU,
        eElasticComputationFlags,
        eSpdCorrectionFlags);
}
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
inline void FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::Serialize(
    io::Archive& archive) const
{
    io::Archive femElastoDynamicsArchive = archive["pbat.sim.dynamics.FemElastoDynamics"];
    mesh.Serialize(femElastoDynamicsArchive);
    bdf.Serialize(femElastoDynamicsArchive);
    femElastoDynamicsArchive.WriteData("fext", fext);
    femElastoDynamicsArchive.WriteData("m", m);
    femElastoDynamicsArchive.WriteData("xtilde", xtilde);
    femElastoDynamicsArchive.WriteData("x", x);
    femElastoDynamicsArchive.WriteData("v", v);
    femElastoDynamicsArchive.WriteData("egU", egU);
    femElastoDynamicsArchive.WriteData("wgU", wgU);
    femElastoDynamicsArchive.WriteData("GNegU", GNegU);
    femElastoDynamicsArchive.WriteData("lamegU", lamegU);
    femElastoDynamicsArchive.WriteData("UgU", UgU);
    femElastoDynamicsArchive.WriteData("GgU", GgU);
    femElastoDynamicsArchive.WriteData("HgU", HgU);
    femElastoDynamicsArchive.WriteMetaData("ndbc", ndbc);
    femElastoDynamicsArchive.WriteData("dbc", dbc);
    femElastoDynamicsArchive.WriteData("dmask", dmask);
}
 
template <
    fem::CElement TElement,
    int Dims,
    physics::CHyperElasticEnergy THyperElasticEnergy,
    common::CFloatingPoint TScalar,
    common::CIndex TIndex>
inline void FemElastoDynamics<TElement, Dims, THyperElasticEnergy, TScalar, TIndex>::Deserialize(
    io::Archive const& archive)
{
    io::Archive const femElastoDynamicsArchive = archive["pbat.sim.dynamics.FemElastoDynamics"];
    mesh.Deserialize(femElastoDynamicsArchive);
    bdf.Deserialize(femElastoDynamicsArchive);
    if (bdf.Order() != 2)
    {
        throw std::runtime_error("FemElastoDynamics only supports BDF of order 2");
    }
    fext = femElastoDynamicsArchive
               .ReadData<Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic>>("fext");
    m      = femElastoDynamicsArchive.ReadData<Eigen::Vector<ScalarType, Eigen::Dynamic>>("m");
    xtilde = femElastoDynamicsArchive
                 .ReadData<Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic>>("xtilde");
    x = femElastoDynamicsArchive
            .ReadData<Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic>>("x");
    v = femElastoDynamicsArchive
            .ReadData<Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic>>("v");
    egU   = femElastoDynamicsArchive.ReadData<Eigen::Vector<IndexType, Eigen::Dynamic>>("egU");
    wgU   = femElastoDynamicsArchive.ReadData<Eigen::Vector<ScalarType, Eigen::Dynamic>>("wgU");
    GNegU = femElastoDynamicsArchive
                .ReadData<Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic>>("GNegU");
    lamegU = femElastoDynamicsArchive
                 .ReadData<Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic>>("lamegU");
    UgU = femElastoDynamicsArchive.ReadData<Eigen::Vector<ScalarType, Eigen::Dynamic>>("UgU");
    GgU = femElastoDynamicsArchive
              .ReadData<Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic>>("GgU");
    HgU = femElastoDynamicsArchive
              .ReadData<Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic>>("HgU");
    ndbc  = femElastoDynamicsArchive.ReadMetaData<IndexType>("ndbc");
    dbc   = femElastoDynamicsArchive.ReadData<Eigen::Vector<IndexType, Eigen::Dynamic>>("dbc");
    dmask = femElastoDynamicsArchive.ReadData<Eigen::Vector<bool, Eigen::Dynamic>>("dmask");
}
 
} // namespace pbat::sim::dynamics
 
#endif // PBAT_SIM_DYNAMICS_FEMELASTODYNAMICS_H