diff --git a/Sofa/Component/Mass/src/sofa/component/mass/ElementFEMMass.cpp b/Sofa/Component/Mass/src/sofa/component/mass/ElementFEMMass.cpp index e2c0f520e30..aae57cce0ec 100644 --- a/Sofa/Component/Mass/src/sofa/component/mass/ElementFEMMass.cpp +++ b/Sofa/Component/Mass/src/sofa/component/mass/ElementFEMMass.cpp @@ -38,6 +38,7 @@ template class SOFA_COMPONENT_MASS_API ElementFEMMass; template class SOFA_COMPONENT_MASS_API ElementFEMMass; template class SOFA_COMPONENT_MASS_API ElementFEMMass; +template class SOFA_COMPONENT_MASS_API ElementFEMMass; void registerFEMMass(sofa::core::ObjectFactory* factory) { @@ -62,6 +63,9 @@ void registerFEMMass(sofa::core::ObjectFactory* factory) factory->registerObjects(sofa::core::ObjectRegistrationData("Finite-element mass (inertia and body force) defined on prisms") .add< ElementFEMMass >(true)); + + factory->registerObjects(sofa::core::ObjectRegistrationData("Finite-element mass (inertia and body force) defined on pyramids") + .add< ElementFEMMass >(true)); } } diff --git a/Sofa/Component/Mass/src/sofa/component/mass/ElementFEMMass.h b/Sofa/Component/Mass/src/sofa/component/mass/ElementFEMMass.h index f390189c375..32eb58459fa 100644 --- a/Sofa/Component/Mass/src/sofa/component/mass/ElementFEMMass.h +++ b/Sofa/Component/Mass/src/sofa/component/mass/ElementFEMMass.h @@ -278,6 +278,7 @@ template class SOFA_COMPONENT_MASS_API ElementFEMMass; template class SOFA_COMPONENT_MASS_API ElementFEMMass; template class SOFA_COMPONENT_MASS_API ElementFEMMass; +template class SOFA_COMPONENT_MASS_API ElementFEMMass; #endif } // namespace sofa::component::mass 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 fd38c9e67bb..d9bd0a7df0c 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 @@ -52,6 +52,9 @@ void registerElementCorotationalFEMForceField(sofa::core::ObjectFactory* factory factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear prisms using the corotational approach") .add< ElementCorotationalFEMForceField >(true)); + + factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear pyramids using the corotational approach") + .add< ElementCorotationalFEMForceField >(true)); } // template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementCorotationalFEMForceField; @@ -64,5 +67,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 484a5c89697..f2dffaaf75e 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 @@ -200,6 +200,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 cb17452bdb2..ddc52d776fe 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 @@ -52,6 +52,9 @@ void registerElementLinearSmallStrainFEMForceField(sofa::core::ObjectFactory* fa factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear prisms assuming small strain") .add< ElementLinearSmallStrainFEMForceField >(true)); + + factory->registerObjects(sofa::core::ObjectRegistrationData("Hooke's law on linear pyramids assuming small strain") + .add< ElementLinearSmallStrainFEMForceField >(true)); } template class SOFA_COMPONENT_SOLIDMECHANICS_FEM_ELASTIC_API ElementLinearSmallStrainFEMForceField; @@ -64,5 +67,5 @@ 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 db8436d83fc..ee81ef8e6ac 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 @@ -107,6 +107,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/Core/src/sofa/core/visual/DrawMesh.h b/Sofa/framework/Core/src/sofa/core/visual/DrawMesh.h index d08445b1fb4..b0f2530fdbc 100644 --- a/Sofa/framework/Core/src/sofa/core/visual/DrawMesh.h +++ b/Sofa/framework/Core/src/sofa/core/visual/DrawMesh.h @@ -442,6 +442,85 @@ struct SOFA_CORE_API DrawElementMesh } }; +template<> +struct SOFA_CORE_API DrawElementMesh + : public BaseDrawMesh, 5> +{ + using ElementType = sofa::geometry::Pyramid; + friend BaseDrawMesh; + + static constexpr ColorContainer defaultColors { + sofa::type::RGBAColor::green(), + sofa::type::RGBAColor::teal(), + sofa::type::RGBAColor::navy(), + sofa::type::RGBAColor::gold(), + sofa::type::RGBAColor::purple() + }; + +private: + template + void doDraw( + sofa::helper::visual::DrawTool* drawTool, + const PositionContainer& position, + sofa::core::topology::BaseMeshTopology* topology, + const IndicesContainer& elementIndices, + const ColorContainer& colors) + { + if (!topology) + return; + + const auto& elements = topology->getPyramids(); + + // Allocate space for rendering points + // 1 Quad + 4 Triangles + renderedPoints[0].resize(elementIndices.size() * sofa::geometry::Quad::NumberOfNodes); + renderedPoints[1].resize(elementIndices.size() * sofa::geometry::Triangle::NumberOfNodes); + renderedPoints[2].resize(elementIndices.size() * sofa::geometry::Triangle::NumberOfNodes); + renderedPoints[3].resize(elementIndices.size() * sofa::geometry::Triangle::NumberOfNodes); + renderedPoints[4].resize(elementIndices.size() * sofa::geometry::Triangle::NumberOfNodes); + + std::array renderedPointId {}; + + for (auto i : elementIndices) + { + const auto& pyramid = elements[i]; + const auto center = this->elementCenter(position, pyramid); + + const auto drawQuad = [&](sofa::Index bufferId, sofa::Index v0, sofa::Index v1, sofa::Index v2, sofa::Index v3) + { + const std::array vertexIndices { pyramid[v0], pyramid[v1], pyramid[v2], pyramid[v3] }; + for (std::size_t k = 0; k < sofa::geometry::Quad::NumberOfNodes; ++k) + { + const auto p = this->applyElementSpace(position[vertexIndices[k]], center); + renderedPoints[bufferId][renderedPointId[bufferId]++] = sofa::type::toVec3(p); + } + }; + + const auto drawTriangle = [&](sofa::Index bufferId, sofa::Index v0, sofa::Index v1, sofa::Index v2) + { + const std::array vertexIndices { pyramid[v0], pyramid[v1], pyramid[v2] }; + for (std::size_t k = 0; k < sofa::geometry::Triangle::NumberOfNodes; ++k) + { + const auto p = this->applyElementSpace(position[vertexIndices[k]], center); + renderedPoints[bufferId][renderedPointId[bufferId]++] = sofa::type::toVec3(p); + } + }; + + drawQuad(0, 0, 3, 2, 1); + drawTriangle(1, 0, 1, 4); + drawTriangle(2, 1, 2, 4); + drawTriangle(3, 3, 4, 2); + drawTriangle(4, 0, 4, 3); + } + + drawTool->drawQuads(renderedPoints[0], colors[0]); + drawTool->drawTriangles(renderedPoints[1], colors[1]); + drawTool->drawTriangles(renderedPoints[2], colors[2]); + drawTool->drawTriangles(renderedPoints[3], colors[3]); + drawTool->drawTriangles(renderedPoints[4], colors[4]); + } +}; + template<> struct SOFA_CORE_API DrawElementMesh : public BaseDrawMesh, 6> @@ -542,6 +621,7 @@ class SOFA_CORE_API DrawMesh drawElements(drawTool, position, topology); drawElements(drawTool, position, topology); drawElements(drawTool, position, topology); + drawElements(drawTool, position, topology); } template @@ -560,8 +640,9 @@ class SOFA_CORE_API DrawMesh const auto hasTetra = !topology->getTetrahedra().empty(); const auto hasHexa = !topology->getHexahedra().empty(); const auto hasPrism = !topology->getPrisms().empty(); + const auto hasPyramid = !topology->getPyramids().empty(); - const bool hasVolumeElements = hasTetra || hasHexa || hasPrism; + const bool hasVolumeElements = hasTetra || hasHexa || hasPrism || hasPyramid; if (!hasSurfaceElements && !hasVolumeElements) { @@ -587,7 +668,8 @@ class SOFA_CORE_API DrawMesh DrawElementMesh, DrawElementMesh, DrawElementMesh, - DrawElementMesh + DrawElementMesh, + DrawElementMesh > m_meshes; }; diff --git a/Sofa/framework/FEM/CMakeLists.txt b/Sofa/framework/FEM/CMakeLists.txt index 222d64f587a..5e88eb85484 100644 --- a/Sofa/framework/FEM/CMakeLists.txt +++ b/Sofa/framework/FEM/CMakeLists.txt @@ -12,6 +12,7 @@ set(HEADER_FILES ${SOFAFEMSRC_ROOT}/FiniteElement[Edge].h ${SOFAFEMSRC_ROOT}/FiniteElement[Hexahedron].h ${SOFAFEMSRC_ROOT}/FiniteElement[Prism].h + ${SOFAFEMSRC_ROOT}/FiniteElement[Pyramid].h ${SOFAFEMSRC_ROOT}/FiniteElement[Quad].h ${SOFAFEMSRC_ROOT}/FiniteElement[Tetrahedron].h ${SOFAFEMSRC_ROOT}/FiniteElement[Triangle].h @@ -23,6 +24,7 @@ set(SOURCE_FILES ${SOFAFEMSRC_ROOT}/FiniteElement[Edge].cpp ${SOFAFEMSRC_ROOT}/FiniteElement[Hexahedron].cpp ${SOFAFEMSRC_ROOT}/FiniteElement[Prism].cpp + ${SOFAFEMSRC_ROOT}/FiniteElement[Pyramid].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[Pyramid].cpp b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].cpp new file mode 100644 index 00000000000..1ad35c584f5 --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].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_PYRAMID_CPP +#include +#include + +namespace sofa::fem +{ + +template struct SOFA_FEM_API FiniteElement; + +} diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h new file mode 100644 index 00000000000..2efcb1613ab --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].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 +#include + +#if !defined(SOFA_FEM_FINITE_ELEMENT_PYRAMID_CPP) +#include +#endif + +namespace sofa::fem +{ + +template +struct FiniteElement +{ + FINITEELEMENT_HEADER(sofa::geometry::Pyramid, DataTypes, 3); + static_assert(spatial_dimensions == 3, "Pyramids are only defined in 3D"); + + constexpr static std::array referenceElementNodes {{ + {-1, -1, -1}, + { 1, -1, -1}, + { 1, 1, -1}, + {-1, 1, -1}, + { 0, 0, 1}, + }}; + + static const sofa::type::vector& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) + { + return topology.getPyramids(); + } + + static constexpr sofa::type::Vec shapeFunctions(const sofa::type::Vec& q) + { + return { + static_cast(0.125) * (1 - q[0]) * (1 - q[1]) * (1 - q[2]), + static_cast(0.125) * (1 + q[0]) * (1 - q[1]) * (1 - q[2]), + static_cast(0.125) * (1 + q[0]) * (1 + q[1]) * (1 - q[2]), + static_cast(0.125) * (1 - q[0]) * (1 + q[1]) * (1 - q[2]), + static_cast(0.5) * (1 + q[2]) + }; + } + + static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) + { + return { + {-static_cast(0.125) * (1 - q[1]) * (1 - q[2]), -static_cast(0.125) * (1 - q[0]) * (1 - q[2]), -static_cast(0.125) * (1 - q[0]) * (1 - q[1])}, + { static_cast(0.125) * (1 - q[1]) * (1 - q[2]), -static_cast(0.125) * (1 + q[0]) * (1 - q[2]), -static_cast(0.125) * (1 + q[0]) * (1 - q[1])}, + { static_cast(0.125) * (1 + q[1]) * (1 - q[2]), static_cast(0.125) * (1 + q[0]) * (1 - q[2]), -static_cast(0.125) * (1 + q[0]) * (1 + q[1])}, + {-static_cast(0.125) * (1 + q[1]) * (1 - q[2]), static_cast(0.125) * (1 - q[0]) * (1 - q[2]), -static_cast(0.125) * (1 - q[0]) * (1 + q[1])}, + { 0, 0, static_cast(0.5)} + }; + } + + static constexpr auto quadraturePoints() + { + constexpr Real sqrt3_1 = static_cast(1) / static_cast(1.73205080757); + constexpr Real one = static_cast(1); + + // We use the 8 Gauss points of the hexahedron, which exactly integrate the (1-z)^2 Jacobian of the pyramid. + return std::array { + std::pair{ReferenceCoord{-sqrt3_1, -sqrt3_1, -sqrt3_1}, one}, + std::pair{ReferenceCoord{ sqrt3_1, -sqrt3_1, -sqrt3_1}, one}, + std::pair{ReferenceCoord{ sqrt3_1, sqrt3_1, -sqrt3_1}, one}, + std::pair{ReferenceCoord{-sqrt3_1, sqrt3_1, -sqrt3_1}, one}, + std::pair{ReferenceCoord{-sqrt3_1, -sqrt3_1, sqrt3_1}, one}, + std::pair{ReferenceCoord{ sqrt3_1, -sqrt3_1, sqrt3_1}, one}, + std::pair{ReferenceCoord{ sqrt3_1, sqrt3_1, sqrt3_1}, one}, + std::pair{ReferenceCoord{-sqrt3_1, sqrt3_1, sqrt3_1}, one}, + }; + } +}; + +#if !defined(SOFA_FEM_FINITE_ELEMENT_PYRAMID_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 c9dd5281428..163df52198c 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[all].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[all].h @@ -24,6 +24,7 @@ #include #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 4665831e6b7..68744ae5f85 100644 --- a/Sofa/framework/FEM/test/FiniteElement_test.cpp +++ b/Sofa/framework/FEM/test/FiniteElement_test.cpp @@ -88,6 +88,11 @@ TEST(FiniteElement, prism3dWeights) testSumWeights(1 / 2.); } +TEST(FiniteElement, pyramid3dWeights) +{ + testSumWeights(8); +} + /** * Checks that the sum of the gradients of shape functions is zero at the evaluation point. */ @@ -169,4 +174,10 @@ TEST(FiniteElement, prism3dGradientShapeFunctions) testGradientShapeFunctions(sofa::type::Vec3(1., 1., 1.)); } +TEST(FiniteElement, pyramid3dGradientShapeFunctions) +{ + testGradientShapeFunctions(sofa::type::Vec3(0., 0., 0.)); + testGradientShapeFunctions(sofa::type::Vec3(0.5, 0.5, 0.5)); +} + } diff --git a/examples/Component/SolidMechanics/FEM/PyramidCorotationalFEMForceField.scn b/examples/Component/SolidMechanics/FEM/PyramidCorotationalFEMForceField.scn new file mode 100644 index 00000000000..7413e1d90a4 --- /dev/null +++ b/examples/Component/SolidMechanics/FEM/PyramidCorotationalFEMForceField.scn @@ -0,0 +1,52 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + +