Skip to content

Commit 1e9e280

Browse files
alxbilgerfredroy
alxbilger
and
fredroy
authored
[FEM] Add support for elasticity of prisms (#6062)
Co-authored-by: Frederick Roy <fredroy@users.noreply.github.com>
1 parent 532d785 commit 1e9e280

10 files changed

Lines changed: 180 additions & 10 deletions

File tree

‎Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.cpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,9 @@ void registerElementCorotationalFEMForceField(sofa::core::ObjectFactory* factory
4949

5050
factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear hexahedra using the corotational approach")
5151
.add< ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron> >(true));
52+
53+
factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear prisms using the corotational approach")
54+
.add< ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Prism> >(true));
5255
}
5356

5457
// template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec1Types, sofa::geometry::Edge>;
@@ -60,5 +63,6 @@ template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotational
6063
template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Quad>;
6164
template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron>;
6265
template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron>;
66+
template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Prism>;
6367

6468
}

‎Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -198,6 +198,7 @@ extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorot
198198
extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Quad>;
199199
extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron>;
200200
extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron>;
201+
extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Prism>;
201202
#endif
202203

203204
} // namespace sofa::component::solidmechanics::fem::elastic

‎Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.cpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,9 @@ void registerElementLinearSmallStrainFEMForceField(sofa::core::ObjectFactory* fa
4949

5050
factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear hexahedra assuming small strain")
5151
.add< ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron> >(true));
52+
53+
factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear prisms assuming small strain")
54+
.add< ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Prism> >(true));
5255
}
5356

5457
template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec1Types, sofa::geometry::Edge>;
@@ -60,5 +63,6 @@ template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallS
6063
template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Quad>;
6164
template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron>;
6265
template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron>;
66+
template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Prism>;
6367

6468
}

‎Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -105,6 +105,7 @@ extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinea
105105
extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Quad>;
106106
extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Tetrahedron>;
107107
extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Hexahedron>;
108+
extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField<sofa::defaulttype::Vec3Types, sofa::geometry::Prism>;
108109
#endif
109110

110111
} // namespace sofa::component::solidmechanics::fem::elastic

‎Sofa/framework/FEM/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@ set(HEADER_FILES
1111
${SOFAFEMSRC_ROOT}/FiniteElement[all].h
1212
${SOFAFEMSRC_ROOT}/FiniteElement[Edge].h
1313
${SOFAFEMSRC_ROOT}/FiniteElement[Hexahedron].h
14+
${SOFAFEMSRC_ROOT}/FiniteElement[Prism].h
1415
${SOFAFEMSRC_ROOT}/FiniteElement[Quad].h
1516
${SOFAFEMSRC_ROOT}/FiniteElement[Tetrahedron].h
1617
${SOFAFEMSRC_ROOT}/FiniteElement[Triangle].h
@@ -21,6 +22,7 @@ set(SOURCE_FILES
2122

2223
${SOFAFEMSRC_ROOT}/FiniteElement[Edge].cpp
2324
${SOFAFEMSRC_ROOT}/FiniteElement[Hexahedron].cpp
25+
${SOFAFEMSRC_ROOT}/FiniteElement[Prism].cpp
2426
${SOFAFEMSRC_ROOT}/FiniteElement[Quad].cpp
2527
${SOFAFEMSRC_ROOT}/FiniteElement[Tetrahedron].cpp
2628
${SOFAFEMSRC_ROOT}/FiniteElement[Triangle].cpp

‎Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].cpp‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
/******************************************************************************
2+
* SOFA, Simulation Open-Framework Architecture *
3+
* (c) 2006 INRIA, USTL, UJF, CNRS, MGH *
4+
* *
5+
* This program is free software; you can redistribute it and/or modify it *
6+
* under the terms of the GNU Lesser General Public License as published by *
7+
* the Free Software Foundation; either version 2.1 of the License, or (at *
8+
* your option) any later version. *
9+
* *
10+
* This program is distributed in the hope that it will be useful, but WITHOUT *
11+
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
12+
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
13+
* for more details. *
14+
* *
15+
* You should have received a copy of the GNU Lesser General Public License *
16+
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
17+
*******************************************************************************
18+
* Authors: The SOFA Team and external contributors (see Authors.txt) *
19+
* *
20+
* Contact information: contact@sofa-framework.org *
21+
******************************************************************************/
22+
#define SOFA_FEM_FINITE_ELEMENT_PRISM_CPP
23+
#include <sofa/fem/FiniteElement[Prism].h>
24+
#include <sofa/defaulttype/VecTypes.h>
25+
26+
namespace sofa::fem
27+
{
28+
29+
template struct SOFA_FEM_API FiniteElement<sofa::geometry::Prism, sofa::defaulttype::Vec3Types>;
30+
31+
}

‎Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h‎

Lines changed: 97 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,97 @@
1+
/******************************************************************************
2+
* SOFA, Simulation Open-Framework Architecture *
3+
* (c) 2006 INRIA, USTL, UJF, CNRS, MGH *
4+
* *
5+
* This program is free software; you can redistribute it and/or modify it *
6+
* under the terms of the GNU Lesser General Public License as published by *
7+
* the Free Software Foundation; either version 2.1 of the License, or (at *
8+
* your option) any later version. *
9+
* *
10+
* This program is distributed in the hope that it will be useful, but WITHOUT *
11+
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
12+
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
13+
* for more details. *
14+
* *
15+
* You should have received a copy of the GNU Lesser General Public License *
16+
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
17+
*******************************************************************************
18+
* Authors: The SOFA Team and external contributors (see Authors.txt) *
19+
* *
20+
* Contact information: contact@sofa-framework.org *
21+
******************************************************************************/
22+
#pragma once
23+
#include <sofa/fem/FiniteElement.h>
24+
25+
#if !defined(SOFA_FEM_FINITE_ELEMENT_PRISM_CPP)
26+
#include <sofa/defaulttype/VecTypes.h>
27+
#endif
28+
29+
namespace sofa::fem
30+
{
31+
32+
template <class DataTypes>
33+
struct FiniteElement<sofa::geometry::Prism, DataTypes>
34+
{
35+
FINITEELEMENT_HEADER(sofa::geometry::Prism, DataTypes, 3);
36+
static_assert(spatial_dimensions == 3, "Prisms are only defined in 3D");
37+
38+
constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
39+
{0, 0, 0},
40+
{1, 0, 0},
41+
{0, 1, 0},
42+
{0, 0, 1},
43+
{1, 0, 1},
44+
{0, 1, 1},
45+
}};
46+
47+
static const sofa::type::vector<TopologyElement>& getElementSequence(sofa::core::topology::BaseMeshTopology& topology)
48+
{
49+
return topology.getPrisms();
50+
}
51+
52+
static constexpr sofa::type::Vec<NumberOfNodesInElement, Real> shapeFunctions(const sofa::type::Vec<TopologicalDimension, Real>& q)
53+
{
54+
return {
55+
q[0] * q[2] - q[0] + q[1] * q[2] - q[1] - q[2] + 1,
56+
(1 - q[2]) * q[0],
57+
(1 - q[2]) * q[1],
58+
(-q[0] - q[1] + 1) * q[2],
59+
q[0] * q[2],
60+
q[1] * q[2]
61+
};
62+
}
63+
64+
static constexpr sofa::type::Mat<NumberOfNodesInElement, TopologicalDimension, Real> gradientShapeFunctions(const sofa::type::Vec<TopologicalDimension, Real>& q)
65+
{
66+
SOFA_UNUSED(q);
67+
return {
68+
{q[2] - 1, q[2] - 1, q[0] + q[1] - 1},
69+
{1 - q[2], 0, -q[0]},
70+
{0, 1 - q[2], -q[1]},
71+
{-q[2], -q[2], -q[0] - q[1] + 1},
72+
{q[2], 0, q[0]},
73+
{0, q[2], q[1]},
74+
};
75+
}
76+
77+
static constexpr std::array<QuadraturePointAndWeight, 2> quadraturePoints()
78+
{
79+
constexpr auto third = static_cast<Real>(1) / static_cast<Real>(3);
80+
constexpr auto sqrt_3 = static_cast<Real>(0.57735026919); // 1/sqrt(3)
81+
constexpr auto one = static_cast<Real>(1);
82+
constexpr QuadraturePoint q0 {third, third, static_cast<Real>(0.5) * (one - sqrt_3)};
83+
constexpr QuadraturePoint q1 {third, third, static_cast<Real>(0.5) * (one + sqrt_3)};
84+
85+
constexpr std::array<QuadraturePointAndWeight, 2> q {
86+
std::make_pair(q0, 1./4.),
87+
std::make_pair(q1, 1./4.),
88+
};
89+
return q;
90+
}
91+
};
92+
93+
#if !defined(SOFA_FEM_FINITE_ELEMENT_PRISM_CPP)
94+
extern template struct SOFA_FEM_API FiniteElement<sofa::geometry::Prism, sofa::defaulttype::Vec3Types>;
95+
#endif
96+
97+
}

‎Sofa/framework/FEM/src/sofa/fem/FiniteElement[all].h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,7 @@
2323

2424
#include <sofa/fem/FiniteElement[Edge].h>
2525
#include <sofa/fem/FiniteElement[Hexahedron].h>
26+
#include <sofa/fem/FiniteElement[Prism].h>
2627
#include <sofa/fem/FiniteElement[Quad].h>
2728
#include <sofa/fem/FiniteElement[Tetrahedron].h>
2829
#include <sofa/fem/FiniteElement[Triangle].h>

‎Sofa/framework/FEM/test/FiniteElement_test.cpp‎

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -83,6 +83,11 @@ TEST(FiniteElement, hexa3dWeights)
8383
testSumWeights<sofa::geometry::Hexahedron, sofa::defaulttype::Vec3Types>(8);
8484
}
8585

86+
TEST(FiniteElement, prism3dWeights)
87+
{
88+
testSumWeights<sofa::geometry::Prism, sofa::defaulttype::Vec3Types>(1 / 2.);
89+
}
90+
8691
/**
8792
* Checks that the sum of the gradients of shape functions is zero at the evaluation point.
8893
*/
@@ -158,6 +163,10 @@ TEST(FiniteElement, hexa3dGradientShapeFunctions)
158163
testGradientShapeFunctions<sofa::geometry::Hexahedron, sofa::defaulttype::Vec3Types>(sofa::type::Vec3(1., 1., 1.));
159164
}
160165

161-
166+
TEST(FiniteElement, prism3dGradientShapeFunctions)
167+
{
168+
testGradientShapeFunctions<sofa::geometry::Prism, sofa::defaulttype::Vec3Types>(sofa::type::Vec3(0., 0., 0.));
169+
testGradientShapeFunctions<sofa::geometry::Prism, sofa::defaulttype::Vec3Types>(sofa::type::Vec3(1., 1., 1.));
170+
}
162171

163172
}

‎examples/Component/Topology/Mapping/Hexa2PrismTopologicalMapping.scn‎

Lines changed: 29 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,12 @@
11
<?xml version="1.0"?>
2-
<Node name="root" dt="0.01" gravity="0 -9 0">
2+
<Node name="root" dt="0.01" gravity="0 -9.81 0">
3+
34
<Node name="plugins">
5+
<RequiredPlugin pluginName="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedProjectiveConstraint] -->
6+
<RequiredPlugin pluginName="Sofa.Component.Engine.Select"/> <!-- Needed to use components [BoxROI] -->
7+
<RequiredPlugin pluginName="Sofa.Component.LinearSolver.Direct"/> <!-- Needed to use components [SparseLDLSolver] -->
8+
<RequiredPlugin pluginName="Sofa.Component.Mass"/> <!-- Needed to use components [MeshMatrixMass] -->
9+
<RequiredPlugin pluginName="Sofa.Component.SolidMechanics.FEM.Elastic"/>
410
<RequiredPlugin pluginName="Sofa.Component.StateContainer"/> <!-- Needed to use components [MechanicalObject] -->
511
<RequiredPlugin pluginName="Sofa.Component.Topology.Container.Constant"/> <!-- Needed to use components [MeshTopology] -->
612
<RequiredPlugin pluginName="Sofa.Component.Topology.Container.Grid"/> <!-- Needed to use components [RegularGridTopology] -->
@@ -9,16 +15,30 @@
915
</Node>
1016

1117
<DefaultAnimationLoop/>
18+
<VisualStyle displayFlags="showBehaviorModels showForceFields" />
19+
20+
<VisualGrid size="0.1"/>
21+
<LineAxis size="0.1"/>
22+
<OglSceneFrame/>
23+
24+
<EulerImplicitSolver name="backward Euler" rayleighStiffness="0.1" rayleighMass="0.1" />
25+
<SparseLDLSolver/>
1226

13-
<Node name="grid">
14-
<RegularGridTopology name="grid" min="-5 -5 0" max="5 5 40" n="5 5 20"/>
15-
<MechanicalObject template="Vec3" name="state" position="@grid.position"/>
27+
<RegularGridTopology name="grid" min="-0.01 -0.01 0" max="0.01 0.01 0.2" n="5 5 30"/>
28+
<MechanicalObject template="Vec3" name="state" showObject="true"/>
1629

17-
<Node name="prisms">
18-
<MeshTopology name="prism_topology"/>
19-
<Hexa2PrismTopologicalMapping input="@grid" output="@prism_topology" />
20-
<VisualMesh position="@../state.position" topology="@prism_topology" enable="true"/>
21-
</Node>
30+
<MeshMatrixMass massDensity="1100"/>
2231

32+
<Node name="prisms">
33+
<MeshTopology name="prism_topology"/>
34+
<Hexa2PrismTopologicalMapping input="@grid" output="@prism_topology" />
35+
36+
<PrismCorotationalFEMForceField name="FEM" youngModulus="2e6" poissonRatio="0.45" topology="@prism_topology"
37+
rotationMethod="polar" computeForceStrategy="sequenced" computeForceDerivStrategy="sequenced"/>
38+
39+
<VisualMesh position="@../state.position" topology="@prism_topology" enable="true"/>
2340
</Node>
41+
42+
<BoxROI template="Vec3" name="box_roi" box="-0.011 -0.011 -0.0001 0.011 0.011 0.0001" drawBoxes="1" />
43+
<FixedProjectiveConstraint template="Vec3" indices="@box_roi.indices" />
2444
</Node>

0 commit comments

Comments
 (0)