Skip to content

Commit daa20ca

Browse files
alxbilgeralxbilger
authored
[SolidMechanics.TensorMass] Implement buildStiffnessMatrix for TetrahedralTensorMassForceField (#4127)
1 parent 0a89f43 commit daa20ca

File tree

3 files changed

+81
-21
lines changed

3 files changed

+81
-21
lines changed

‎Sofa/Component/SolidMechanics/TensorMass/src/sofa/component/solidmechanics/tensormass/TetrahedralTensorMassForceField.h‎

Lines changed: 2 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -58,27 +58,12 @@ class TetrahedralTensorMassForceField : public core::behavior::ForceField<DataTy
5858

5959
using Index = sofa::Index;
6060

61-
class Mat3 : public type::fixed_array<Deriv,3>
62-
{
63-
public:
64-
Deriv operator*(const Deriv& v) const
65-
{
66-
return Deriv((*this)[0]*v,(*this)[1]*v,(*this)[2]*v);
67-
}
68-
Deriv transposeMultiply(const Deriv& v) const
69-
{
70-
return Deriv(v[0]*((*this)[0])[0]+v[1]*((*this)[1])[0]+v[2]*((*this)[2][0]),
71-
v[0]*((*this)[0][1])+v[1]*((*this)[1][1])+v[2]*((*this)[2][1]),
72-
v[0]*((*this)[0][2])+v[1]*((*this)[1][2])+v[2]*((*this)[2][2]));
73-
}
74-
};
75-
7661
protected:
7762

7863
class EdgeRestInformation
7964
{
8065
public:
81-
Mat3 DfDx; /// the edge stiffness matrix
66+
sofa::type::Mat<3, 3, Real> DfDx; /// the edge stiffness matrix
8267
float vertices[2]; // the vertices of this edge
8368

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

124109
void addForce(const core::MechanicalParams* mparams, DataVecDeriv& d_f, const DataVecCoord& d_x, const DataVecDeriv& d_v) override;
125110
void addDForce(const core::MechanicalParams* mparams, DataVecDeriv& d_df, const DataVecDeriv& d_dx) override;
111+
void buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix) override;
126112
void buildDampingMatrix(core::behavior::DampingMatrix* /*matrix*/) final;
127113
SReal getPotentialEnergy(const core::MechanicalParams* /*mparams*/, const DataVecCoord& /* x */) const override
128114
{

‎Sofa/Component/SolidMechanics/TensorMass/src/sofa/component/solidmechanics/tensormass/TetrahedralTensorMassForceField.inl‎

Lines changed: 32 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -104,7 +104,7 @@ void TetrahedralTensorMassForceField<DataTypes>::applyTetrahedronCreation(const
104104
k = m_topology->getLocalEdgesInTetrahedron(j)[0];
105105
l = m_topology->getLocalEdgesInTetrahedron(j)[1];
106106

107-
Mat3 &m=edgeData[te[j]].DfDx;
107+
auto& m = edgeData[te[j]].DfDx;
108108

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

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

188-
Mat3 &m=edgeData[te[j]].DfDx;
188+
auto& m = edgeData[te[j]].DfDx;
189189

190190
val1= dot(shapeVector[k],shapeVector[l])*mustar;
191191
// print if obtuse tetrahedron along that edge
@@ -365,7 +365,7 @@ SReal TetrahedralTensorMassForceField<DataTypes>::getPotentialEnergy(const core
365365
dp = dp1-dp0;
366366
force=einfo->DfDx*dp;
367367
energy+=dot(force,dp1);
368-
force=einfo->DfDx.transposeMultiply(dp);
368+
force=einfo->DfDx.multTranspose(dp);
369369
energy-=dot(force,dp0);
370370
}
371371

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

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

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

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

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

454+
template <class DataTypes>
455+
void TetrahedralTensorMassForceField<DataTypes>::buildStiffnessMatrix(sofa::core::behavior::StiffnessMatrix* matrix)
456+
{
457+
auto dfdx = matrix->getForceDerivativeIn(this->mstate)
458+
.withRespectToPositionsIn(this->mstate);
459+
const sofa::Size nbEdges = m_topology->getNbEdges();
460+
461+
const auto edgeInf = sofa::helper::getReadAccessor(edgeInfo);
462+
const auto edges = m_topology->getEdges();
463+
464+
for (sofa::Size i = 0; i < nbEdges; ++i)
465+
{
466+
const sofa::topology::Edge& edge = edges[i];
467+
468+
const unsigned p0 = Deriv::total_size * edge[0];
469+
const unsigned p1 = Deriv::total_size * edge[1];
470+
471+
const auto& localDfdx = edgeInf[i].DfDx;
472+
const auto localDfdx_T = edgeInf[i].DfDx.transposed();
473+
474+
dfdx(p0, p0) += localDfdx_T;
475+
dfdx(p0, p1) += -localDfdx_T;
476+
dfdx(p1, p0) += -localDfdx;
477+
dfdx(p1, p1) += localDfdx;
478+
}
479+
}
480+
454481
template <class DataTypes>
455482
void TetrahedralTensorMassForceField<DataTypes>::buildDampingMatrix(core::behavior::DampingMatrix*)
456483
{

‎examples/Component/SolidMechanics/TensorMass/TetrahedralTensorMassForceField.scn‎

Lines changed: 47 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,47 @@
1+
<Node name="root">
2+
<Node name="plugins">
3+
<RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedConstraint] -->
4+
<RequiredPlugin name="Sofa.Component.Engine.Select"/> <!-- Needed to use components [BoxROI] -->
5+
<RequiredPlugin name="Sofa.Component.LinearSolver.Direct"/> <!-- Needed to use components [EigenSimplicialLDLT] -->
6+
<RequiredPlugin name="Sofa.Component.Mapping.Linear"/> <!-- Needed to use components [IdentityMapping] -->
7+
<RequiredPlugin name="Sofa.Component.Mass"/> <!-- Needed to use components [UniformMass] -->
8+
<RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <!-- Needed to use components [EulerImplicitSolver] -->
9+
<RequiredPlugin name="Sofa.Component.SolidMechanics.TensorMass"/> <!-- Needed to use components [TetrahedralTensorMassForceField] -->
10+
<RequiredPlugin name="Sofa.Component.StateContainer"/> <!-- Needed to use components [MechanicalObject] -->
11+
<RequiredPlugin name="Sofa.Component.Topology.Container.Dynamic"/> <!-- Needed to use components [QuadSetGeometryAlgorithms,QuadSetTopologyContainer,QuadSetTopologyModifier,TetrahedronSetGeometryAlgorithms,TetrahedronSetTopologyContainer,TetrahedronSetTopologyModifier] -->
12+
<RequiredPlugin name="Sofa.Component.Topology.Container.Grid"/> <!-- Needed to use components [RegularGridTopology] -->
13+
<RequiredPlugin name="Sofa.Component.Topology.Mapping"/> <!-- Needed to use components [Hexa2QuadTopologicalMapping,Hexa2TetraTopologicalMapping] -->
14+
<RequiredPlugin name="Sofa.Component.Visual"/> <!-- Needed to use components [VisualStyle] -->
15+
<RequiredPlugin name="Sofa.GL.Component.Rendering3D"/> <!-- Needed to use components [OglModel] -->
16+
</Node>
17+
18+
<VisualStyle displayFlags="showBehaviorModels showForceFields" />
19+
20+
<DefaultAnimationLoop name="animationLoop"/>
21+
<DefaultVisualManagerLoop name="visualLoop"/>
22+
23+
<EulerImplicitSolver name="odesolver" rayleighStiffness="0.1" rayleighMass="0.1" />
24+
<EigenSimplicialLDLT template="CompressedRowSparseMatrixMat3x3"/>
25+
<MechanicalObject name="DoFs" />
26+
<UniformMass name="mass" totalMass="320" />
27+
<RegularGridTopology name="grid" nx="4" ny="4" nz="20" xmin="-9" xmax="-6" ymin="0" ymax="3" zmin="0" zmax="19" />
28+
<BoxROI name="box" box="-10 -1 -0.0001 -5 4 0.0001"/>
29+
<FixedConstraint indices="@box.indices" />
30+
31+
<TetrahedronSetTopologyContainer name="Tetra_topo"/>
32+
<TetrahedronSetTopologyModifier name="Modifier" />
33+
<TetrahedronSetGeometryAlgorithms template="Vec3" name="GeomAlgo" />
34+
<Hexa2TetraTopologicalMapping input="@grid" output="@Tetra_topo" />
35+
<TetrahedralTensorMassForceField name="deformable" youngModulus="100000" poissonRatio="0.4" />
36+
37+
<Node name="quads">
38+
<QuadSetTopologyContainer name="Container" />
39+
<QuadSetTopologyModifier name="Modifier" />
40+
<QuadSetGeometryAlgorithms name="GeomAlgo" template="Vec3" />
41+
<Hexa2QuadTopologicalMapping input="@../grid" output="@Container" />
42+
<Node name="Visu">
43+
<OglModel name="Visual" color="yellow" />
44+
<IdentityMapping input="@../../DoFs" output="@Visual" />
45+
</Node>
46+
</Node>
47+
</Node>

0 commit comments

Comments
 (0)