Skip to content

Commit

Permalink
[SolidMechanics.TensorMass] Implement buildStiffnessMatrix for Tetrah…
Browse files Browse the repository at this point in the history 
…edralTensorMassForceField (#4127)
  • Loading branch information
@alxbilger
alxbilger authored Sep 6, 2023
1 parent 0a89f43 commit daa20ca
Show file tree
Hide file tree
Showing 3 changed files with 81 additions and 21 deletions.
18 changes: 2 additions & 16 deletions ...TensorMass/src/sofa/component/solidmechanics/tensormass/TetrahedralTensorMassForceField.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,27 +58,12 @@ class TetrahedralTensorMassForceField : public core::behavior::ForceField<DataTy

using Index = sofa::Index;

class Mat3 : public type::fixed_array<Deriv,3>
{
public:
Deriv operator*(const Deriv& v) const
{
return Deriv((*this)[0]*v,(*this)[1]*v,(*this)[2]*v);
}
Deriv transposeMultiply(const Deriv& v) const
{
return Deriv(v[0]*((*this)[0])[0]+v[1]*((*this)[1])[0]+v[2]*((*this)[2][0]),
v[0]*((*this)[0][1])+v[1]*((*this)[1][1])+v[2]*((*this)[2][1]),
v[0]*((*this)[0][2])+v[1]*((*this)[1][2])+v[2]*((*this)[2][2]));
}
};

protected:

class EdgeRestInformation
{
public:
Mat3 DfDx; /// the edge stiffness matrix
sofa::type::Mat<3, 3, Real> DfDx; /// the edge stiffness matrix
float vertices[2]; // the vertices of this edge

EdgeRestInformation()
Expand Down Expand Up @@ -123,6 +108,7 @@ class TetrahedralTensorMassForceField : public core::behavior::ForceField<DataTy

void addForce(const core::MechanicalParams* mparams, DataVecDeriv& d_f, const DataVecCoord& d_x, const DataVecDeriv& d_v) override;
void addDForce(const core::MechanicalParams* mparams, DataVecDeriv& d_df, const DataVecDeriv& d_dx) override;
void buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) override;
void buildDampingMatrix(core::behavior::DampingMatrix* /*matrix*/) final;
SReal getPotentialEnergy(const core::MechanicalParams* /*mparams*/, const DataVecCoord& /* x */) const override
{
Expand Down
37 changes: 32 additions & 5 deletions ...nsorMass/src/sofa/component/solidmechanics/tensormass/TetrahedralTensorMassForceField.inl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@ void TetrahedralTensorMassForceField<DataTypes>::applyTetrahedronCreation(const
k = m_topology->getLocalEdgesInTetrahedron(j)[0];
l = m_topology->getLocalEdgesInTetrahedron(j)[1];

Mat3 &m=edgeData[te[j]].DfDx;
auto& m = edgeData[te[j]].DfDx;

val1= dot(shapeVector[k],shapeVector[l])*mustar;

Expand Down Expand Up @@ -185,7 +185,7 @@ void TetrahedralTensorMassForceField<DataTypes>::applyTetrahedronDestruction(con
k = m_topology->getLocalEdgesInTetrahedron(j)[0];
l = m_topology->getLocalEdgesInTetrahedron(j)[1];

Mat3 &m=edgeData[te[j]].DfDx;
auto& m = edgeData[te[j]].DfDx;

val1= dot(shapeVector[k],shapeVector[l])*mustar;
// print if obtuse tetrahedron along that edge
Expand Down Expand Up @@ -365,7 +365,7 @@ SReal TetrahedralTensorMassForceField<DataTypes>::getPotentialEnergy(const core
dp = dp1-dp0;
force=einfo->DfDx*dp;
energy+=dot(force,dp1);
force=einfo->DfDx.transposeMultiply(dp);
force=einfo->DfDx.multTranspose(dp);
energy-=dot(force,dp0);
}

Expand Down Expand Up @@ -404,7 +404,7 @@ void TetrahedralTensorMassForceField<DataTypes>::addForce(const core::Mechanical
dp = dp1-dp0;

f[v1]+=einfo->DfDx*dp;
f[v0]-=einfo->DfDx.transposeMultiply(dp);
f[v0]-=einfo->DfDx.multTranspose(dp);
}

edgeInfo.endEdit();
Expand Down Expand Up @@ -442,7 +442,7 @@ void TetrahedralTensorMassForceField<DataTypes>::addDForce(const core::Mechanica
dp = dp1-dp0;

df[v1]+= (einfo->DfDx*dp) * kFactor;
df[v0]-= (einfo->DfDx.transposeMultiply(dp)) * kFactor;
df[v0]-= (einfo->DfDx.multTranspose(dp)) * kFactor;
}
edgeInfo.endEdit();

Expand All @@ -451,6 +451,33 @@ void TetrahedralTensorMassForceField<DataTypes>::addDForce(const core::Mechanica
sofa::helper::AdvancedTimer::stepEnd("addDForceTetraTensorMass");
}

template <class DataTypes>
void TetrahedralTensorMassForceField<DataTypes>::buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix)
{
auto dfdx = matrix->getForceDerivativeIn(this->mstate)
.withRespectToPositionsIn(this->mstate);
const sofa::Size nbEdges = m_topology->getNbEdges();

const auto edgeInf = sofa::helper::getReadAccessor(edgeInfo);
const auto edges = m_topology->getEdges();

for (sofa::Size i = 0; i < nbEdges; ++i)
{
const sofa::topology::Edge& edge = edges[i];

const unsigned p0 = Deriv::total_size * edge[0];
const unsigned p1 = Deriv::total_size * edge[1];

const auto& localDfdx = edgeInf[i].DfDx;
const auto localDfdx_T = edgeInf[i].DfDx.transposed();

dfdx(p0, p0) += localDfdx_T;
dfdx(p0, p1) += -localDfdx_T;
dfdx(p1, p0) += -localDfdx;
dfdx(p1, p1) += localDfdx;
}
}

template <class DataTypes>
void TetrahedralTensorMassForceField<DataTypes>::buildDampingMatrix(core::behavior::DampingMatrix*)
{
Expand Down
47 changes: 47 additions & 0 deletions examples/Component/SolidMechanics/TensorMass/TetrahedralTensorMassForceField.scn
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@
<Node name="root">
<Node name="plugins">
<RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedConstraint] -->
<RequiredPlugin name="Sofa.Component.Engine.Select"/> <!-- Needed to use components [BoxROI] -->
<RequiredPlugin name="Sofa.Component.LinearSolver.Direct"/> <!-- Needed to use components [EigenSimplicialLDLT] -->
<RequiredPlugin name="Sofa.Component.Mapping.Linear"/> <!-- Needed to use components [IdentityMapping] -->
<RequiredPlugin name="Sofa.Component.Mass"/> <!-- Needed to use components [UniformMass] -->
<RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <!-- Needed to use components [EulerImplicitSolver] -->
<RequiredPlugin name="Sofa.Component.SolidMechanics.TensorMass"/> <!-- Needed to use components [TetrahedralTensorMassForceField] -->
<RequiredPlugin name="Sofa.Component.StateContainer"/> <!-- Needed to use components [MechanicalObject] -->
<RequiredPlugin name="Sofa.Component.Topology.Container.Dynamic"/> <!-- Needed to use components [QuadSetGeometryAlgorithms,QuadSetTopologyContainer,QuadSetTopologyModifier,TetrahedronSetGeometryAlgorithms,TetrahedronSetTopologyContainer,TetrahedronSetTopologyModifier] -->
<RequiredPlugin name="Sofa.Component.Topology.Container.Grid"/> <!-- Needed to use components [RegularGridTopology] -->
<RequiredPlugin name="Sofa.Component.Topology.Mapping"/> <!-- Needed to use components [Hexa2QuadTopologicalMapping,Hexa2TetraTopologicalMapping] -->
<RequiredPlugin name="Sofa.Component.Visual"/> <!-- Needed to use components [VisualStyle] -->
<RequiredPlugin name="Sofa.GL.Component.Rendering3D"/> <!-- Needed to use components [OglModel] -->
</Node>

<VisualStyle displayFlags="showBehaviorModels showForceFields" />

<DefaultAnimationLoop name="animationLoop"/>
<DefaultVisualManagerLoop name="visualLoop"/>

<EulerImplicitSolver name="odesolver" rayleighStiffness="0.1" rayleighMass="0.1" />
<EigenSimplicialLDLT template="CompressedRowSparseMatrixMat3x3"/>
<MechanicalObject name="DoFs" />
<UniformMass name="mass" totalMass="320" />
<RegularGridTopology name="grid" nx="4" ny="4" nz="20" xmin="-9" xmax="-6" ymin="0" ymax="3" zmin="0" zmax="19" />
<BoxROI name="box" box="-10 -1 -0.0001 -5 4 0.0001"/>
<FixedConstraint indices="@box.indices" />

<TetrahedronSetTopologyContainer name="Tetra_topo"/>
<TetrahedronSetTopologyModifier name="Modifier" />
<TetrahedronSetGeometryAlgorithms template="Vec3" name="GeomAlgo" />
<Hexa2TetraTopologicalMapping input="@grid" output="@Tetra_topo" />
<TetrahedralTensorMassForceField name="deformable" youngModulus="100000" poissonRatio="0.4" />

<Node name="quads">
<QuadSetTopologyContainer name="Container" />
<QuadSetTopologyModifier name="Modifier" />
<QuadSetGeometryAlgorithms name="GeomAlgo" template="Vec3" />
<Hexa2QuadTopologicalMapping input="@../grid" output="@Container" />
<Node name="Visu">
<OglModel name="Visual" color="yellow" />
<IdentityMapping input="@../../DoFs" output="@Visual" />
</Node>
</Node>
</Node>

0 comments on commit daa20ca

Please sign in to comment.