diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.cpp index 8ac1f227d65..fd38c9e67bb 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.cpp +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.cpp @@ -49,6 +49,9 @@ void registerElementCorotationalFEMForceField(sofa::core::ObjectFactory* factory factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear hexahedra using the corotational approach") .add< ElementCorotationalFEMForceField >(true)); + + factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear prisms using the corotational approach") + .add< ElementCorotationalFEMForceField >(true)); } // template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; @@ -60,5 +63,6 @@ template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotational template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; } diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h index 7bfc3c303f0..8184a56dffe 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementCorotationalFEMForceField.h @@ -198,6 +198,7 @@ extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorot extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; #endif } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.cpp b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.cpp index 97dcff5f064..cb17452bdb2 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.cpp +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.cpp @@ -49,6 +49,9 @@ void registerElementLinearSmallStrainFEMForceField(sofa::core::ObjectFactory* fa factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear hexahedra assuming small strain") .add< ElementLinearSmallStrainFEMForceField >(true)); + + factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear prisms assuming small strain") + .add< ElementLinearSmallStrainFEMForceField >(true)); } template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; @@ -60,5 +63,6 @@ template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallS template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; } diff --git a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h index 74d17801b29..d61a536ec05 100644 --- a/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h +++ b/Sofa/Component/SolidMechanics/FEM/Elastic/src/sofa/component/solidmechanics/fem/elastic/ElementLinearSmallStrainFEMForceField.h @@ -105,6 +105,7 @@ extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinea extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; +extern template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; #endif } // namespace sofa::component::solidmechanics::fem::elastic diff --git a/Sofa/framework/FEM/CMakeLists.txt b/Sofa/framework/FEM/CMakeLists.txt index 1ebf55e9e25..222d64f587a 100644 --- a/Sofa/framework/FEM/CMakeLists.txt +++ b/Sofa/framework/FEM/CMakeLists.txt @@ -11,6 +11,7 @@ set(HEADER_FILES ${SOFAFEMSRC_ROOT}/FiniteElement[all].h ${SOFAFEMSRC_ROOT}/FiniteElement[Edge].h ${SOFAFEMSRC_ROOT}/FiniteElement[Hexahedron].h + ${SOFAFEMSRC_ROOT}/FiniteElement[Prism].h ${SOFAFEMSRC_ROOT}/FiniteElement[Quad].h ${SOFAFEMSRC_ROOT}/FiniteElement[Tetrahedron].h ${SOFAFEMSRC_ROOT}/FiniteElement[Triangle].h @@ -21,6 +22,7 @@ set(SOURCE_FILES ${SOFAFEMSRC_ROOT}/FiniteElement[Edge].cpp ${SOFAFEMSRC_ROOT}/FiniteElement[Hexahedron].cpp + ${SOFAFEMSRC_ROOT}/FiniteElement[Prism].cpp ${SOFAFEMSRC_ROOT}/FiniteElement[Quad].cpp ${SOFAFEMSRC_ROOT}/FiniteElement[Tetrahedron].cpp ${SOFAFEMSRC_ROOT}/FiniteElement[Triangle].cpp diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].cpp b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].cpp new file mode 100644 index 00000000000..27c68a3fb65 --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].cpp @@ -0,0 +1,31 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define SOFA_FEM_FINITE_ELEMENT_PRISM_CPP +#include +#include + +namespace sofa::fem +{ + +template struct SOFA_FEM_API FiniteElement; + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h new file mode 100644 index 00000000000..5a6731ff94c --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h @@ -0,0 +1,97 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include + +#if !defined(SOFA_FEM_FINITE_ELEMENT_PRISM_CPP) +#include +#endif + +namespace sofa::fem +{ + +template +struct FiniteElement +{ + FINITEELEMENT_HEADER(sofa::geometry::Prism, DataTypes, 3); + static_assert(spatial_dimensions == 3, "Prisms are only defined in 3D"); + + constexpr static std::array referenceElementNodes {{ + {0, 0, 0}, + {1, 0, 0}, + {0, 1, 0}, + {0, 0, 1}, + {1, 0, 1}, + {0, 1, 1}, + }}; + + static const sofa::type::vector& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) + { + return topology.getPrisms(); + } + + static constexpr sofa::type::Vec shapeFunctions(const sofa::type::Vec& q) + { + return { + q[0] * q[2] - q[0] + q[1] * q[2] - q[1] - q[2] + 1, + (1 - q[2]) * q[0], + (1 - q[2]) * q[1], + (-q[0] - q[1] + 1) * q[2], + q[0] * q[2], + q[1] * q[2] + }; + } + + static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) + { + SOFA_UNUSED(q); + return { + {q[2] - 1, q[2] - 1, q[0] + q[1] - 1}, + {1 - q[2], 0, -q[0]}, + {0, 1 - q[2], -q[1]}, + {-q[2], -q[2], -q[0] - q[1] + 1}, + {q[2], 0, q[0]}, + {0, q[2], q[1]}, + }; + } + + static constexpr std::array quadraturePoints() + { + constexpr auto third = static_cast(1) / static_cast(3); + constexpr auto sqrt_3 = static_cast(0.57735026919); // 1/sqrt(3) + constexpr auto one = static_cast(1); + constexpr QuadraturePoint q0 {third, third, static_cast(0.5) * (one - sqrt_3)}; + constexpr QuadraturePoint q1 {third, third, static_cast(0.5) * (one + sqrt_3)}; + + constexpr std::array q { + std::make_pair(q0, 1./4.), + std::make_pair(q1, 1./4.), + }; + return q; + } +}; + +#if !defined(SOFA_FEM_FINITE_ELEMENT_PRISM_CPP) +extern template struct SOFA_FEM_API FiniteElement; +#endif + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[all].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[all].h index eafac60913c..c9dd5281428 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[all].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[all].h @@ -23,6 +23,7 @@ #include #include +#include #include #include #include diff --git a/Sofa/framework/FEM/test/FiniteElement_test.cpp b/Sofa/framework/FEM/test/FiniteElement_test.cpp index dfe53b92218..4665831e6b7 100644 --- a/Sofa/framework/FEM/test/FiniteElement_test.cpp +++ b/Sofa/framework/FEM/test/FiniteElement_test.cpp @@ -83,6 +83,11 @@ TEST(FiniteElement, hexa3dWeights) testSumWeights(8); } +TEST(FiniteElement, prism3dWeights) +{ + testSumWeights(1 / 2.); +} + /** * Checks that the sum of the gradients of shape functions is zero at the evaluation point. */ @@ -158,6 +163,10 @@ TEST(FiniteElement, hexa3dGradientShapeFunctions) testGradientShapeFunctions(sofa::type::Vec3(1., 1., 1.)); } - +TEST(FiniteElement, prism3dGradientShapeFunctions) +{ + testGradientShapeFunctions(sofa::type::Vec3(0., 0., 0.)); + testGradientShapeFunctions(sofa::type::Vec3(1., 1., 1.)); +} } diff --git a/examples/Component/Topology/Mapping/Hexa2PrismTopologicalMapping.scn b/examples/Component/Topology/Mapping/Hexa2PrismTopologicalMapping.scn index 00196cbadc0..76442156aa3 100644 --- a/examples/Component/Topology/Mapping/Hexa2PrismTopologicalMapping.scn +++ b/examples/Component/Topology/Mapping/Hexa2PrismTopologicalMapping.scn @@ -1,24 +1,49 @@ - + + + + + + + + + + + - + - - + + + + + + + + + + + - - - + + - - - - - + + + + + + + + + + +