Skip to content

Commit

Permalink
[SofaGeneralSimpleFem] Add comments and tests for TriangularFEMForceF…
Browse files Browse the repository at this point in the history
…ieldOptim (#2284)

* [SofaMiscFEM] Add test on TriangularFEMForceFieldOptim

* [SofaGeneralSimpleFem] Fix init of TriangularFEMForceFieldOptim

* [SofaGeneralSimpleFem] Fix AddForce method to TriangularFEMForceFieldOptim, the strainDisplacement was not computed on the current local position but always on the local init values therefore the strain and force were slightly wrong.

* [SofaGeneralSimpleFem] Fix AddDForce method of TriangularFEMForceFieldOptim, using saved StrainDisplacement value from AddForce. To be tested

* [SofaGeneralSimpleFem] Update AddDForce method of TriangularFEMForceFieldOptim. Should be better if I well understand what has been done....

* restore previous version

* Add scene to tests the different triangular forcefield

* [SofaGeneralSimpleFem] Add getter to element info to match TriangleFEM api. Those method should not be used internally

* [SofaMiscFem_test] Rename variable rotInit into rotatedInitPos which is more accurate

* [SofaMiscFem_test] Fix checkInit values for TriangularFEMForceFieldOptim

* [SofaMiscFem_test] Fix checkValues for TriangularFEMForceFieldOptim. Still some little diffs

* Update examples/Benchmark/Accuracy/TriangleFEMForceField_compare.scn

Co-authored-by: Frederick Roy <fredroy@users.noreply.github.com>

* Update examples/Benchmark/Accuracy/TriangleFEMForceField_compare.scn

Co-authored-by: Hugo <hugo.talbot@sofa-framework.org>

* revert some changes

Co-authored-by: Frederick Roy <fredroy@users.noreply.github.com>
Co-authored-by: Hugo <hugo.talbot@sofa-framework.org>
  • Loading branch information
3 people authored Dec 9, 2021
1 parent 6865a7e commit ba9b76b
Show file tree
Hide file tree
Showing 5 changed files with 494 additions and 69 deletions.
102 changes: 102 additions & 0 deletions examples/Benchmark/Accuracy/TriangleFEMForceField_compare.scn
Original file line number Diff line number Diff line change
@@ -0,0 +1,102 @@
<?xml version="1.0" ?>
<Node name="root" dt="0.01" gravity="0 10 0">
<RequiredPlugin name="SofaBoundaryCondition"/> <!-- Needed to use components [FixedConstraint] -->
<RequiredPlugin name="SofaGeneralSimpleFem"/> <!-- Needed to use components [TriangularFEMForceFieldOptim] -->
<RequiredPlugin name="SofaImplicitOdeSolver"/> <!-- Needed to use components [EulerImplicitSolver] -->
<RequiredPlugin name="SofaMiscFem"/> <!-- Needed to use components [TriangleFEMForceField, TriangularFEMForceField] -->
<RequiredPlugin name="SofaOpenglVisual"/> <!-- Needed to use components [OglModel] -->
<RequiredPlugin name="SofaTopologyMapping"/> <!-- Needed to use components [Quad2TriangleTopologicalMapping] -->

<VisualStyle displayFlags="showBehaviorModels showVisual" />
<DefaultPipeline verbose="0" />


<RegularGridTopology name="grid" n="40 1 40" min="0 0 0" max="10 0 10" />

<Node name="TriangularFEMForceField">
<EulerImplicitSolver name="cg_odesolver" printLog="false" rayleighStiffness="0.1" rayleighMass="0.1" />
<CGLinearSolver iterations="25" name="linear solver" tolerance="1.0e-9" threshold="1.0e-9" />

<QuadSetTopologyContainer name="Quad_topo" src="@../grid" />
<QuadSetTopologyModifier name="Modifier" />
<QuadSetGeometryAlgorithms name="GeomAlgo" template="Vec3d" />

<MechanicalObject name="TriangularFEMForceField_dof" />
<DiagonalMass totalMass="1.0" />
<FixedConstraint indices="0 39" />

<Node name="T">
<TriangleSetTopologyContainer name="Triangle_topo" />
<TriangleSetTopologyModifier name="Modifier" />
<TriangleSetGeometryAlgorithms name="GeomAlgo" template="Vec3d" />

<Quad2TriangleTopologicalMapping name="TriangularFEMForceField_mapTopo" input="@../Quad_topo" output="@Triangle_topo" />
<TriangularFEMForceField name="FEM" youngModulus="1000" poissonRatio="0.3" method="large" />

<Node name="Visu">
<OglModel name="TriangularFEMForceField_visu" color="red" />
<IdentityMapping name="TriangularFEMForceField_mapping" input="@../.." output="@." />
</Node>
</Node>

</Node>


<Node name="TriangularFEMForceFieldOptim">
<EulerImplicitSolver name="cg_odesolver" printLog="false" rayleighStiffness="0.1" rayleighMass="0.1" />
<CGLinearSolver iterations="25" name="linear solver" tolerance="1.0e-9" threshold="1.0e-9" />

<QuadSetTopologyContainer name="Quad_topo" src="@../grid" />
<QuadSetTopologyModifier name="Modifier" />
<QuadSetGeometryAlgorithms name="GeomAlgo" template="Vec3d" />

<MechanicalObject name="TriangularFEMForceFieldOptim_dof" />
<DiagonalMass totalMass="1.0" />
<FixedConstraint indices="0 39" />

<Node name="T">
<TriangleSetTopologyContainer name="Triangle_topo" />
<TriangleSetTopologyModifier name="Modifier" />
<TriangleSetGeometryAlgorithms name="GeomAlgo" template="Vec3d" />

<Quad2TriangleTopologicalMapping name="TriangularFEMForceFieldOptim_mapTopo" input="@../Quad_topo" output="@Triangle_topo" />
<TriangularFEMForceFieldOptim name="FEM" youngModulus="1000" poissonRatio="0.3" method="large" />

<Node name="Visu">
<OglModel name="TriangularFEMForceFieldOptim_visu" color="blue" />
<IdentityMapping name="TriangularFEMForceFieldOptim_mapping" input="@../.." output="@." />
</Node>
</Node>

</Node>


<Node name="TriangleFEMForceField">
<EulerImplicitSolver name="cg_odesolver" printLog="false" rayleighStiffness="0.1" rayleighMass="0.1" />
<CGLinearSolver iterations="25" name="linear solver" tolerance="1.0e-9" threshold="1.0e-9" />

<QuadSetTopologyContainer name="Quad_topo" src="@../grid" />
<QuadSetTopologyModifier name="Modifier" />
<QuadSetGeometryAlgorithms name="GeomAlgo" template="Vec3d" />

<MechanicalObject name="TriangleFEMForceField_dof" />
<DiagonalMass totalMass="1.0" />
<FixedConstraint indices="0 39" />

<Node name="T">
<TriangleSetTopologyContainer name="Triangle_topo" />
<TriangleSetTopologyModifier name="Modifier" />
<TriangleSetGeometryAlgorithms name="GeomAlgo" template="Vec3d" />

<Quad2TriangleTopologicalMapping name="TriangleFEMForceField_mapTopo" input="@../Quad_topo" output="@Triangle_topo" />
<TriangleFEMForceField name="FEM" youngModulus="1000" poissonRatio="0.3" method="large" />

<Node name="Visu">
<OglModel name="TriangleFEMForceField_visu" color="green" />
<IdentityMapping name="TriangleFEMForceField_mapping" input="@../.." output="@." />
</Node>
</Node>

</Node>

</Node>
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,7 @@ class TriangularFEMForceFieldOptim : public core::behavior::ForceField<DataTypes

protected:
typedef type::Mat<2, 3, Real > Transformation; ///< matrix for rigid transformations like rotations
typedef type::Mat<3, 3, Real> MaterialStiffness;
enum { DerivSize = DataTypes::deriv_total_size };
typedef type::Mat<DerivSize, DerivSize, Real> MatBloc;

Expand All @@ -102,6 +103,9 @@ class TriangularFEMForceFieldOptim : public core::behavior::ForceField<DataTypes

virtual ~TriangularFEMForceFieldOptim();
public:
Real getPoisson() { return d_poisson.getValue(); }
Real getYoung() { return d_young.getValue(); }

void init() override;
void reinit() override;
void addForce(const core::MechanicalParams* mparams, DataVecDeriv& f, const DataVecCoord& x, const DataVecDeriv& v) override;
Expand Down Expand Up @@ -264,6 +268,13 @@ class TriangularFEMForceFieldOptim : public core::behavior::ForceField<DataTypes
void getTriangleVonMisesStress(Index i, Real& stressValue);
void getTrianglePrincipalStress(Index i, Real& stressValue, Deriv& stressDirection, Real& stressValue2, Deriv& stressDirection2);

/// Public methods to access FEM information per element. Those method should not be used internally as they add check on element id.
type::fixed_array <Coord, 3> getRotatedInitialElement(Index elemId);
Transformation getRotationMatrix(Index elemId);
MaterialStiffness getMaterialStiffness(Index elemId);
type::Vec3 getStrainDisplacementFactors(Index elemId);
Real getTriangleFactor(Index elemId);

public:

/// Forcefield intern paramaters
Expand Down
Loading

0 comments on commit ba9b76b

Please sign in to comment.