[FEM] Introduce higher accuracy quadrature rules

Sofa/Component/Mass/Testing/src/sofa/component/mass/testing/MassTestCreation.h

@@ -32,6 +32,8 @@
 #include <sofa/testing/BaseSimulationTest.h>
 #include <sofa/testing/NumericTest.h>
 
+#include <algorithm>
+
 namespace sofa::component::mass::testing
 {
 
@@ -131,7 +133,10 @@ struct Mass_test : public sofa::testing::BaseSimulationTest, public sofa::testin
         sofa::type::vector<Real> expectedMv;
         sofa::testing::data_traits<DataTypes>::VecDeriv_to_Vector(expectedMv, Mv);
 
-        EXPECT_LT(this->vectorMaxDiff(eigenMv, expectedMv), m_errorMax * this->epsilon())
+        const SReal tolerance = m_errorMax * (SReal)this->epsilon()
+            * std::max((SReal)1, (SReal)this->vectorMaxAbs(expectedMv));
+
+        EXPECT_LT((SReal)this->vectorMaxDiff(eigenMv, expectedMv), tolerance)
             << "addMToMatrix inconsistent with addMDx";
     }
 
@@ -167,7 +172,10 @@ struct Mass_test : public sofa::testing::BaseSimulationTest, public sofa::testin
         sofa::type::vector<Real> expectedMv;
         sofa::testing::data_traits<DataTypes>::VecDeriv_to_Vector(expectedMv, Mv);
 
-        EXPECT_LT((SReal)this->vectorMaxDiff(eigenMv, expectedMv), (SReal)(m_errorMax * this->epsilon()))
+        const SReal tolerance = m_errorMax * (SReal)this->epsilon()
+            * std::max((SReal)1, (SReal)this->vectorMaxAbs(expectedMv));
+
+        EXPECT_LT((SReal)this->vectorMaxDiff(eigenMv, expectedMv), tolerance)
             << "buildMassMatrix inconsistent with addMDx";
     }
 
@@ -181,8 +189,11 @@ struct Mass_test : public sofa::testing::BaseSimulationTest, public sofa::testin
             vMv += static_cast<SReal>(sofa::type::dot(v[i], Mv[i]));
         }
 
+        const SReal tolerance = m_errorMax * (SReal)this->epsilon()
+            * std::max((SReal)1, std::abs(vMv));
+
         // v * Mv should be the 2 * kinetic energy
-        EXPECT_LT(std::abs(vMv - 2 * ke), (SReal)(m_errorMax * this->epsilon()))
+        EXPECT_LT(std::abs(vMv - 2 * ke), tolerance)
             << "Kinetic energy inconsistent with addMDx (vMv = " << vMv << ", kinetic energy = " << (Real)ke << ")";
     }
     /**
@@ -236,7 +247,11 @@ struct Mass_test : public sofa::testing::BaseSimulationTest, public sofa::testin
             VecDeriv a;
             sofa::testing::copyFromData(a, m_dof->readForces());
 
-            EXPECT_LT((SReal)this->vectorMaxDiff(a, v), (SReal)(m_errorMax * this->epsilon())) << "accFromF inconsistent with addMDx (M^-1 * M * v != v)";
+            const SReal tolerance = m_errorMax * (SReal)this->epsilon()
+                * std::max((SReal)1, (SReal)this->vectorMaxAbs(v));
+
+            EXPECT_LT((SReal)this->vectorMaxDiff(a, v), tolerance)
+                << "accFromF inconsistent with addMDx (M^-1 * M * v != v)";
         }
 
         sofa::SignedIndex matrixSize = (sofa::SignedIndex)(n * DataTypes::deriv_total_size);

Sofa/Component/Mass/tests/MassTestCreation[FEMMass].cpp

@@ -48,7 +48,7 @@ struct ElementMass_template_test : public Mass_test<FEMMass<typename MassParam::
     {
         this->m_testAccFromF = false;
         this->m_testAddMToMatrix = false;
-        this->m_errorMax = 1e3_sreal;
+        this->m_errorMax = 1e2_sreal;
 
         auto topology = sofa::core::objectmodel::New<MeshTopology>();
         this->m_node->addObject(topology);

Sofa/framework/FEM/CMakeLists.txt

@@ -9,6 +9,7 @@ set(HEADER_FILES
 
     ${SOFAFEMSRC_ROOT}/FiniteElement.h
     ${SOFAFEMSRC_ROOT}/FiniteElement[all].h
+    ${SOFAFEMSRC_ROOT}/QuadratureRules.h
 
     ${SOFAFEMSRC_ROOT}/FiniteElement[Edge].h
     ${SOFAFEMSRC_ROOT}/FiniteElement[Hexahedron].h

Sofa/framework/FEM/src/sofa/fem/FiniteElement.h

@@ -21,6 +21,7 @@
 ******************************************************************************/
 #pragma once
 #include <sofa/fem/config.h>
+#include <sofa/fem/QuadratureRules.h>
 #include <sofa/core/topology/BaseMeshTopology.h>
 
 namespace sofa::fem
@@ -29,7 +30,7 @@ namespace sofa::fem
 template <class ElementType, class DataTypes>
 struct FiniteElement;
 
-#define FINITEELEMENT_HEADER(ElType, DataTypes, dimension) \
+#define FINITEELEMENT_HEADER(ElType, DataTypes, dimension, defaultQuadratureDegree) \
     using Coord = sofa::Coord_t<DataTypes>;\
     using Real = sofa::Real_t<DataTypes>;\
     using ElementType = ElType;\
@@ -41,6 +42,15 @@ struct FiniteElement;
     using ShapeFunctionType = std::function<Real(const ReferenceCoord&)>;\
     using QuadraturePoint = ReferenceCoord; \
     using QuadraturePointAndWeight = std::pair<QuadraturePoint, Real>;\
+    using Quadrature = FiniteElementQuadrature<ElType, Real, defaultQuadratureDegree>;\
+    static constexpr sofa::Size InterpolationOrder = Quadrature::InterpolationOrder;\
+    static constexpr sofa::Size MinimumQuadratureDegree = Quadrature::MinimumQuadratureDegree;\
+    static constexpr sofa::Size DefaultQuadratureDegree = Quadrature::DefaultQuadratureDegree;\
+    template <sofa::Size Degree = Quadrature::DefaultQuadratureDegree>\
+    static constexpr auto quadraturePoints()\
+        { return Quadrature::template quadraturePoints<Degree>(); }\
+    static auto quadratureRule(sofa::Size degree)\
+        { return Quadrature::quadratureRule(degree); }\
     using Helper = FiniteElementHelper<ElementType, DataTypes>
 
 

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

@@ -32,7 +32,7 @@ namespace sofa::fem
 template <class DataTypes>
 struct FiniteElement<sofa::geometry::Edge, DataTypes>
 {
-    FINITEELEMENT_HEADER(sofa::geometry::Edge, DataTypes, 1);
+    FINITEELEMENT_HEADER(sofa::geometry::Edge, DataTypes, 1, 1);
 
     constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{ReferenceCoord{-1}, ReferenceCoord{1}}};
 
@@ -54,15 +54,6 @@ struct FiniteElement<sofa::geometry::Edge, DataTypes>
         SOFA_UNUSED(q);
         return {{-static_cast<Real>(0.5)}, {static_cast<Real>(0.5)}};
     }
-
-    static constexpr std::array<QuadraturePointAndWeight, 1> quadraturePoints()
-    {
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q0(static_cast<Real>(0));
-        return {
-            std::make_pair(q0, static_cast<Real>(2))
-        };
-    }
-
 };
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_EDGE_CPP)

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

@@ -32,7 +32,7 @@ namespace sofa::fem
 template <class DataTypes>
 struct FiniteElement<sofa::geometry::Hexahedron, DataTypes>
 {
-    FINITEELEMENT_HEADER(sofa::geometry::Hexahedron, DataTypes, 3);
+    FINITEELEMENT_HEADER(sofa::geometry::Hexahedron, DataTypes, 3, 3);
     static_assert(spatial_dimensions == 3, "Hexahedrons are only defined in 3D");
 
     // Following the convention in sofa::geometry::Hexahedron:
@@ -95,26 +95,6 @@ struct FiniteElement<sofa::geometry::Hexahedron, DataTypes>
 
         return gradient;
     }
-
-    static constexpr auto quadraturePoints()
-    {
-        constexpr Real sqrt3 = 1.73205080757; //sqrt(3.)
-        constexpr Real sqrt3_1 = static_cast<Real>(1) / sqrt3;
-        constexpr Real one = static_cast<Real>(1);
-
-        constexpr std::array q {
-            std::pair{referenceElementNodes[0] * sqrt3_1, one},
-            std::pair{referenceElementNodes[1] * sqrt3_1, one},
-            std::pair{referenceElementNodes[2] * sqrt3_1, one},
-            std::pair{referenceElementNodes[3] * sqrt3_1, one},
-            std::pair{referenceElementNodes[4] * sqrt3_1, one},
-            std::pair{referenceElementNodes[5] * sqrt3_1, one},
-            std::pair{referenceElementNodes[6] * sqrt3_1, one},
-            std::pair{referenceElementNodes[7] * sqrt3_1, one},
-        };
-
-        return q;
-    }
 };
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_HEXAHEDRON_CPP)

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

@@ -32,7 +32,7 @@ namespace sofa::fem
 template <class DataTypes>
 struct FiniteElement<sofa::geometry::Prism, DataTypes>
 {
-    FINITEELEMENT_HEADER(sofa::geometry::Prism, DataTypes, 3);
+    FINITEELEMENT_HEADER(sofa::geometry::Prism, DataTypes, 3, 1);
     static_assert(spatial_dimensions == 3, "Prisms are only defined in 3D");
 
     constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
@@ -73,21 +73,6 @@ struct FiniteElement<sofa::geometry::Prism, DataTypes>
             {0, q[2], q[1]},
         };
     }
-
-    static constexpr std::array<QuadraturePointAndWeight, 2> quadraturePoints()
-    {
-        constexpr auto third = static_cast<Real>(1) / static_cast<Real>(3);
-        constexpr auto sqrt_3 = static_cast<Real>(0.57735026919); // 1/sqrt(3)
-        constexpr auto one = static_cast<Real>(1);
-        constexpr QuadraturePoint q0 {third, third, static_cast<Real>(0.5) * (one - sqrt_3)};
-        constexpr QuadraturePoint q1 {third, third, static_cast<Real>(0.5) * (one + sqrt_3)};
-
-        constexpr std::array<QuadraturePointAndWeight, 2> q {
-            std::make_pair(q0, 1./4.),
-            std::make_pair(q1, 1./4.),
-        };
-        return q;
-    }
 };
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_PRISM_CPP)

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

@@ -33,7 +33,7 @@ namespace sofa::fem
 template <class DataTypes>
 struct FiniteElement<sofa::geometry::Pyramid, DataTypes>
 {
-    FINITEELEMENT_HEADER(sofa::geometry::Pyramid, DataTypes, 3);
+    FINITEELEMENT_HEADER(sofa::geometry::Pyramid, DataTypes, 3, 1);
     static_assert(spatial_dimensions == 3, "Pyramids are only defined in 3D");
 
     constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
@@ -70,24 +70,6 @@ struct FiniteElement<sofa::geometry::Pyramid, DataTypes>
             { 0,                                                 0,                                                 static_cast<Real>(0.5)}
         };
     }
-
-    static constexpr auto quadraturePoints()
-    {
-        constexpr Real sqrt3_1 = static_cast<Real>(1) / static_cast<Real>(1.73205080757);
-        constexpr Real one = static_cast<Real>(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)

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

@@ -32,7 +32,7 @@ namespace sofa::fem
 template <class DataTypes>
 struct FiniteElement<sofa::geometry::Quad, DataTypes>
 {
-    FINITEELEMENT_HEADER(sofa::geometry::Quad, DataTypes, 2);
+    FINITEELEMENT_HEADER(sofa::geometry::Quad, DataTypes, 2, 2);
     static_assert(spatial_dimensions > 1, "Quads cannot be defined in 1D");
 
     constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
@@ -66,23 +66,6 @@ struct FiniteElement<sofa::geometry::Quad, DataTypes>
             {1 / static_cast<Real>(4) * (-static_cast<Real>(1) - q[1]), 1 / static_cast<Real>(4) * ( static_cast<Real>(1) - q[0])}
         };
     }
-
-    static constexpr std::array<QuadraturePointAndWeight, 3> quadraturePoints()
-    {
-        constexpr Real sqrt2_3 = 0.816496580928; //sqrt(2./3.)
-        constexpr Real sqrt6 = 2.44948974278; //sqrt(6.)
-        constexpr Real sqrt2 = 1.41421356237; //sqrt(2.)
-
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q0(sqrt2_3, 0.);
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q1(-1/sqrt6, -1./sqrt2);
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q2(-1/sqrt6, 1./sqrt2);
-
-        return {
-            std::make_pair(q0, 4./3.),
-            std::make_pair(q1, 4./3.),
-            std::make_pair(q2, 4./3.),
-        };
-    }
 };
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_QUAD_CPP)

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

@@ -32,12 +32,12 @@ namespace sofa::fem
 template <class DataTypes>
 struct FiniteElement<sofa::geometry::QuadraticEdge, DataTypes>
 {
-    FINITEELEMENT_HEADER(sofa::geometry::QuadraticEdge, DataTypes, 1);
+    FINITEELEMENT_HEADER(sofa::geometry::QuadraticEdge, DataTypes, 1, 3);
 
     constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
-        ReferenceCoord{0},    // vertex 0
+        ReferenceCoord{-1},   // vertex 0
         ReferenceCoord{1},    // vertex 1
-        ReferenceCoord{0.5}   // mid-edge node
+        ReferenceCoord{0}     // mid-edge node
     }};
 
     static const sofa::type::vector<TopologyElement>& getElementSequence(sofa::core::topology::BaseMeshTopology& topology)
@@ -48,40 +48,24 @@ struct FiniteElement<sofa::geometry::QuadraticEdge, DataTypes>
     static constexpr sofa::type::Vec<NumberOfNodesInElement, Real> shapeFunctions(const sofa::type::Vec<TopologicalDimension, Real>& q)
     {
         const Real xi = q[0];
+        constexpr Real half = static_cast<Real>(0.5);
         return {
-            2 * xi * (xi - 0.5) + 1 - 2 * xi,  // vertex 0: (2*xi - 1) * (xi - 1) = 2*xi^2 - 3*xi + 1
-            2 * xi * (xi - 0.5),      // vertex 1: (2*xi - 1) * xi
-            4 * xi * (1 - xi)         // mid-edge: 4*xi*(1-xi)
+            half * xi * (xi - 1),                  // vertex 0, at xi = -1
+            half * xi * (xi + 1),                  // vertex 1, at xi = 1
+            static_cast<Real>(1) - xi * xi         // mid-edge node, at xi = 0
         };
     }
 
     static constexpr sofa::type::Mat<NumberOfNodesInElement, TopologicalDimension, Real> gradientShapeFunctions(const sofa::type::Vec<TopologicalDimension, Real>& q)
     {
         const Real xi = q[0];
+        constexpr Real half = static_cast<Real>(0.5);
         return {
-            {4 * xi - 3},    // vertex 0
-            {4 * xi - 1},    // vertex 1
-            {4 - 8 * xi}     // mid-edge
+            {xi - half},                   // vertex 0
+            {xi + half},                   // vertex 1
+            {-2 * xi}                      // mid-edge node
         };
     }
-
-    static constexpr std::array<QuadraturePointAndWeight, 2> quadraturePoints()
-    {
-        // constexpr Real a = (1. - 1. / std::sqrt(3.)) / 2.;
-        constexpr Real a { 0.211324865405 };
-        // constexpr Real b = (1. + 1. / std::sqrt(3.)) / 2.;
-        constexpr Real b { 0.788675134595 };
-        constexpr Real w = 0.5;
-
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q0(a);
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q1(b);
-
-        constexpr std::array<QuadraturePointAndWeight, 2> q {
-            std::make_pair(q0, w),
-            std::make_pair(q1, w)
-        };
-        return q;
-    }
 };
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_EDGE_CPP)

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

@@ -32,7 +32,7 @@ namespace sofa::fem
 template <class DataTypes>
 struct FiniteElement<sofa::geometry::QuadraticHexahedron, DataTypes>
 {
-    FINITEELEMENT_HEADER(sofa::geometry::QuadraticHexahedron, DataTypes, 3);
+    FINITEELEMENT_HEADER(sofa::geometry::QuadraticHexahedron, DataTypes, 3, 5);
     static_assert(spatial_dimensions == 3, "Quadratic Hexahedrons are only defined in 3D");
 
     constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
@@ -219,33 +219,6 @@ struct FiniteElement<sofa::geometry::QuadraticHexahedron, DataTypes>
             N_0, N_1, N_2, N_3, N_4, N_5, N_6, N_7, N_8, N_9, N_10, N_11, N_12, N_13, N_14, N_15,
             N_16, N_17, N_18, N_19, N_20, N_21, N_22, N_23, N_24, N_25, N_26);
     }
-
-    static constexpr std::array<QuadraturePointAndWeight, 27> quadraturePoints()
-    {
-        constexpr Real a = 0.7745966692414834; // sqrt(3/5)
-        constexpr Real w_a = 5.0 / 9.0;
-        constexpr Real w_0 = 8.0 / 9.0;
-
-        constexpr std::array<Real, 3> pts = {-a, 0.0, a};
-        constexpr std::array<Real, 3> wts = {w_a, w_0, w_a};
-
-        std::array<QuadraturePointAndWeight, 27> q {};
-        int index = 0;
-        for (int i = 0; i < 3; ++i)
-        {
-            for (int j = 0; j < 3; ++j)
-            {
-                for (int k = 0; k < 3; ++k)
-                {
-                    q[index++] = std::make_pair(
-                        sofa::type::Vec<TopologicalDimension, Real>(pts[i], pts[j], pts[k]),
-                        wts[i] * wts[j] * wts[k]
-                    );
-                }
-            }
-        }
-        return q;
-    }
 };
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_HEXAHEDRON_CPP)