diff --git a/Sofa/Component/Mass/Testing/src/sofa/component/mass/testing/MassTestCreation.h b/Sofa/Component/Mass/Testing/src/sofa/component/mass/testing/MassTestCreation.h
index 0dbc331d251..e706e7999f7 100644
--- a/Sofa/Component/Mass/Testing/src/sofa/component/mass/testing/MassTestCreation.h
+++ b/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);
diff --git a/Sofa/Component/Mass/tests/MassTestCreation[FEMMass].cpp b/Sofa/Component/Mass/tests/MassTestCreation[FEMMass].cpp
index a4dadaaa65c..caa8d9bee0e 100644
--- a/Sofa/Component/Mass/tests/MassTestCreation[FEMMass].cpp
+++ b/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);
diff --git a/Sofa/framework/FEM/CMakeLists.txt b/Sofa/framework/FEM/CMakeLists.txt
index 183b38a9d87..a529117e814 100644
--- a/Sofa/framework/FEM/CMakeLists.txt
+++ b/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
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h
index a490fddd6a3..7723e86aec5 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h
+++ b/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>
 
 
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h
index fbffb5dbfed..2afaaf46be0 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h
+++ b/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)
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h
index 38339624e47..bbb7a52111d 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h
+++ b/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)
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h
index 5a6731ff94c..f05ad349363 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h
+++ b/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)
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h
index 2efcb1613ab..b83318b7a92 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h
+++ b/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)
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h
index 4e37585eace..a7580dad7c7 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h
+++ b/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)
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h
index 4606613f486..f94d68691ac 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h
+++ b/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)
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].h
index 733b39707ba..875b02c4a26 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].h
+++ b/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)
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h
index e1f312c21f5..0cbd2a43f3a 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h
+++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h
@@ -32,7 +32,7 @@ namespace sofa::fem
 template <class DataTypes>
 struct FiniteElement<sofa::geometry::QuadraticQuad, DataTypes>
 {
-    FINITEELEMENT_HEADER(sofa::geometry::QuadraticQuad, DataTypes, 2);
+    FINITEELEMENT_HEADER(sofa::geometry::QuadraticQuad, DataTypes, 2, 3);
 
     constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
         {-1, -1},  // vertex 0
@@ -95,26 +95,6 @@ struct FiniteElement<sofa::geometry::QuadraticQuad, DataTypes>
             {-2 * xi * (1 - eta * eta), -2 * eta * (1 - xi * xi)}
         };
     }
-
-    static constexpr std::array<QuadraturePointAndWeight, 4> quadraturePoints()
-    {
-        // constexpr Real a = 1. / std::sqrt(3.);
-        constexpr Real a { 0.57735026919 };
-        constexpr Real w = 1.0;
-
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q0(-a, -a);
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q1(a, -a);
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q2(a, a);
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q3(-a, a);
-
-        constexpr std::array<QuadraturePointAndWeight, 4> q {
-            std::make_pair(q0, w),
-            std::make_pair(q1, w),
-            std::make_pair(q2, w),
-            std::make_pair(q3, w)
-        };
-        return q;
-    }
 };
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_QUAD_CPP)
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].h
index 2cd90ee677e..f9f45b28a9a 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].h
+++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].h
@@ -32,7 +32,7 @@ namespace sofa::fem
 template <class DataTypes>
 struct FiniteElement<sofa::geometry::QuadraticTetrahedron, DataTypes>
 {
-    FINITEELEMENT_HEADER(sofa::geometry::QuadraticTetrahedron, DataTypes, 3);
+    FINITEELEMENT_HEADER(sofa::geometry::QuadraticTetrahedron, DataTypes, 3, 2);
     static_assert(spatial_dimensions == 3, "Quadratic Tetrahedrons are only defined in 3D");
 
     constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
@@ -104,29 +104,6 @@ struct FiniteElement<sofa::geometry::QuadraticTetrahedron, DataTypes>
             {0, 4 * l3, 4 * l2}
         };
     }
-
-    static constexpr std::array<QuadraturePointAndWeight, 4> quadraturePoints()
-    {
-        // 4-point quadrature rule for quadratic tetrahedron
-        // constexpr Real a = (5. + 3. * std::sqrt(5.)) / 20.;
-        constexpr Real a { 0.585410196625 };
-        // constexpr Real b = (5. - std::sqrt(5.)) / 20.;
-        constexpr Real b { 0.138196601125 };
-        constexpr Real w = 1. / 24.;
-
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q0(a, b, b);
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q1(b, a, b);
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q2(b, b, a);
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q3(b, b, b);
-
-        constexpr std::array<QuadraturePointAndWeight, 4> q {
-            std::make_pair(q0, w),
-            std::make_pair(q1, w),
-            std::make_pair(q2, w),
-            std::make_pair(q3, w)
-        };
-        return q;
-    }
 };
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_TETAHEDRON_CPP)
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].h
index c77ebd73cac..3063a0679bb 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].h
+++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].h
@@ -32,7 +32,7 @@ namespace sofa::fem
 template <class DataTypes>
 struct FiniteElement<sofa::geometry::QuadraticTriangle, DataTypes>
 {
-    FINITEELEMENT_HEADER(sofa::geometry::QuadraticTriangle, DataTypes, 2);
+    FINITEELEMENT_HEADER(sofa::geometry::QuadraticTriangle, DataTypes, 2, 2);
 
     constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
         {0, 0},      // vertex 0
@@ -85,24 +85,6 @@ struct FiniteElement<sofa::geometry::QuadraticTriangle, DataTypes>
             {-4 * l2, 4 * (l0 - l2)}
         };
     }
-
-    static constexpr std::array<QuadraturePointAndWeight, 3> quadraturePoints()
-    {
-        constexpr Real a = 2. / 3.;
-        constexpr Real b = 1. / 6.;
-        constexpr Real w = 1. / 6.;
-
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q0(b, b);
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q1(a, b);
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q2(b, a);
-
-        constexpr std::array<QuadraturePointAndWeight, 3> q {
-            std::make_pair(q0, w),
-            std::make_pair(q1, w),
-            std::make_pair(q2, w)
-        };
-        return q;
-    }
 };
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_TRIANGLE_CPP)
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h
index ac14bb9c59b..ec074d17e89 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h
+++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h
@@ -32,7 +32,7 @@ namespace sofa::fem
 template <class DataTypes>
 struct FiniteElement<sofa::geometry::Tetrahedron, DataTypes>
 {
-    FINITEELEMENT_HEADER(sofa::geometry::Tetrahedron, DataTypes, 3);
+    FINITEELEMENT_HEADER(sofa::geometry::Tetrahedron, DataTypes, 3, 1);
     static_assert(spatial_dimensions == 3, "Tetrahedrons are only defined in 3D");
 
     constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
@@ -67,13 +67,6 @@ struct FiniteElement<sofa::geometry::Tetrahedron, DataTypes>
             {0, 0, 1}
         };
     }
-
-    static constexpr std::array<QuadraturePointAndWeight, 1> quadraturePoints()
-    {
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q0(1./4., 1./4., 1./4.);
-        constexpr std::array<QuadraturePointAndWeight, 1> q { std::make_pair(q0, 1./6.) };
-        return q;
-    }
 };
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_TETAHEDRON_CPP)
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h
index 0d9659040de..83465f57ae0 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h
+++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h
@@ -32,7 +32,7 @@ namespace sofa::fem
 template <class DataTypes>
 struct FiniteElement<sofa::geometry::Triangle, DataTypes>
 {
-    FINITEELEMENT_HEADER(sofa::geometry::Triangle, DataTypes, 2);
+    FINITEELEMENT_HEADER(sofa::geometry::Triangle, DataTypes, 2, 1);
     static_assert(spatial_dimensions > 1, "Triangles cannot be defined in 1D");
 
     constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
@@ -63,13 +63,6 @@ struct FiniteElement<sofa::geometry::Triangle, DataTypes>
             {0, 1}
         };
     }
-
-    static constexpr std::array<QuadraturePointAndWeight, 1> quadraturePoints()
-    {
-        return {
-            std::make_pair(sofa::type::Vec<TopologicalDimension, Real>(1./3., 1./3.), 1./2.)
-        };
-    }
 };
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_TRIANGLE_CPP)
diff --git a/Sofa/framework/FEM/src/sofa/fem/QuadratureRules.h b/Sofa/framework/FEM/src/sofa/fem/QuadratureRules.h
new file mode 100644
index 00000000000..60417b88311
--- /dev/null
+++ b/Sofa/framework/FEM/src/sofa/fem/QuadratureRules.h
@@ -0,0 +1,543 @@
+/******************************************************************************
+*                 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 <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Authors: The SOFA Team and external contributors (see Authors.txt)          *
+*                                                                             *
+* Contact information: contact@sofa-framework.org                             *
+******************************************************************************/
+#pragma once
+#include <sofa/fem/config.h>
+
+#include <sofa/geometry/Edge.h>
+#include <sofa/geometry/Triangle.h>
+#include <sofa/geometry/Quad.h>
+#include <sofa/geometry/Tetrahedron.h>
+#include <sofa/geometry/Hexahedron.h>
+#include <sofa/geometry/Prism.h>
+#include <sofa/geometry/Pyramid.h>
+#include <sofa/geometry/QuadraticElements.h>
+#include <sofa/type/Vec.h>
+
+#include <array>
+#include <numbers>
+#include <span>
+#include <stdexcept>
+#include <string>
+#include <utility>
+
+namespace sofa::fem
+{
+
+/// Mapping element types to their integration reference domain. 
+/// This domain is the same irrespective of the order of the element.
+template <class ElementType>
+struct ReferenceDomain
+{
+    using type = ElementType;
+};
+
+template <> struct ReferenceDomain<sofa::geometry::QuadraticEdge>        { using type = sofa::geometry::Edge; };
+template <> struct ReferenceDomain<sofa::geometry::QuadraticTriangle>    { using type = sofa::geometry::Triangle; };
+template <> struct ReferenceDomain<sofa::geometry::QuadraticQuad>        { using type = sofa::geometry::Quad; };
+template <> struct ReferenceDomain<sofa::geometry::QuadraticTetrahedron> { using type = sofa::geometry::Tetrahedron; };
+template <> struct ReferenceDomain<sofa::geometry::QuadraticHexahedron>  { using type = sofa::geometry::Hexahedron; };
+
+template <class ElementType>
+using ReferenceDomain_t = typename ReferenceDomain<ElementType>::type;
+
+/// Always false, but dependent on its argument: will trigger compilation failure in case no rule
+/// has been tabulated for the requested degree.
+template <sofa::Size>
+inline constexpr bool degreeIsNotTabulated = false;
+
+/**
+ * Quadrature rules over the reference domain of a topological element.
+ *
+ * An incomplete type, followed by specializations. Each specialization exposes:
+ *  - TopologicalDimension, the dimension of the reference domain;
+ *  - tabulatedDegrees, the polynomial degrees for which a table is provided, ascending;
+ *  - points<Degree, Real>(), the table for one of those degrees.
+ *
+ * A request for a degree that is not tabulated is served by the smallest tabulated
+ * degree above it; see QuadratureDispatch.
+ */
+template <class Domain>
+struct QuadratureRules;
+
+template <class Domain, class Real>
+using QuadraturePointAndWeight_t =
+    std::pair<sofa::type::Vec<QuadratureRules<Domain>::TopologicalDimension, Real>, Real>;
+
+/// Reference domain: the segment [-1, 1].
+template <>
+struct QuadratureRules<sofa::geometry::Edge>
+{
+    static constexpr sofa::Size TopologicalDimension = 1;
+    static constexpr std::array<sofa::Size, 2> tabulatedDegrees { 1, 3 };
+
+    template <sofa::Size Degree, class Real>
+    static constexpr auto points()
+    {
+        using Point = sofa::type::Vec<TopologicalDimension, Real>;
+        using PointAndWeight = std::pair<Point, Real>;
+
+        if constexpr (Degree == 1)
+        {
+            // 1-point midpoint rule
+            return std::array<PointAndWeight, 1>{
+                std::make_pair(Point(static_cast<Real>(0)), static_cast<Real>(2))
+            };
+        }
+        else if constexpr (Degree == 3)
+        {
+            // 2-point Gauss-Legendre rule
+            constexpr Real g = std::numbers::inv_sqrt3_v<Real>; // 1/sqrt(3)
+            return std::array<PointAndWeight, 2>{
+                std::make_pair(Point(-g), static_cast<Real>(1)),
+                std::make_pair(Point( g), static_cast<Real>(1))
+            };
+        }
+        else
+        {
+            static_assert(degreeIsNotTabulated<Degree>, "QuadratureRules<Edge>: degree is not tabulated");
+        }
+    }
+};
+
+/// Reference domain: the unit triangle (0,0), (1,0), (0,1).
+template <>
+struct QuadratureRules<sofa::geometry::Triangle>
+{
+    static constexpr sofa::Size TopologicalDimension = 2;
+    static constexpr std::array<sofa::Size, 2> tabulatedDegrees { 1, 2 };
+
+    template <sofa::Size Degree, class Real>
+    static constexpr auto points()
+    {
+        using Point = sofa::type::Vec<TopologicalDimension, Real>;
+        using PointAndWeight = std::pair<Point, Real>;
+
+        constexpr Real third = static_cast<Real>(1) / static_cast<Real>(3); // 1/3
+        constexpr Real sixth = static_cast<Real>(1) / static_cast<Real>(6); // 1/6
+
+        if constexpr (Degree == 1)
+        {
+            // 1-point centroid rule
+            return std::array<PointAndWeight, 1>{
+                std::make_pair(Point(third, third), static_cast<Real>(0.5)) // area of the reference triangle
+            };
+        }
+        else if constexpr (Degree == 2)
+        {
+            // 3-point interior rule
+            constexpr Real a = static_cast<Real>(2) / static_cast<Real>(3); // 2/3
+            return std::array<PointAndWeight, 3>{
+                std::make_pair(Point(sixth, sixth), sixth),
+                std::make_pair(Point(a,     sixth), sixth),
+                std::make_pair(Point(sixth, a),     sixth)
+            };
+        }
+        else
+        {
+            static_assert(degreeIsNotTabulated<Degree>, "QuadratureRules<Triangle>: degree is not tabulated");
+        }
+    }
+};
+
+/// Reference domain: the square [-1, 1]^2.
+template <>
+struct QuadratureRules<sofa::geometry::Quad>
+{
+    static constexpr sofa::Size TopologicalDimension = 2;
+    static constexpr std::array<sofa::Size, 3> tabulatedDegrees { 1, 2, 3 };
+
+    template <sofa::Size Degree, class Real>
+    static constexpr auto points()
+    {
+        using Point = sofa::type::Vec<TopologicalDimension, Real>;
+        using PointAndWeight = std::pair<Point, Real>;
+
+        if constexpr (Degree == 1)
+        {
+            // 1-point centroid rule
+            return std::array<PointAndWeight, 1>{
+                std::make_pair(Point(static_cast<Real>(0), static_cast<Real>(0)), static_cast<Real>(4))
+            };
+        }
+        else if constexpr (Degree == 2)
+        {
+            // 3-point rule
+            constexpr Real sqrt2_3 = std::numbers::sqrt2_v<Real> * std::numbers::inv_sqrt3_v<Real>; // sqrt(2/3)
+            constexpr Real inv_sqrt6 = std::numbers::inv_sqrt3_v<Real> / std::numbers::sqrt2_v<Real>; // 1/sqrt(6)
+            constexpr Real inv_sqrt2 = std::numbers::sqrt2_v<Real> / static_cast<Real>(2); // 1/sqrt(2)
+            constexpr Real w = static_cast<Real>(4) / static_cast<Real>(3); // 4/3
+
+            return std::array<PointAndWeight, 3>{
+                std::make_pair(Point( sqrt2_3,    static_cast<Real>(0)), w),
+                std::make_pair(Point(-inv_sqrt6, -inv_sqrt2),            w),
+                std::make_pair(Point(-inv_sqrt6,  inv_sqrt2),            w)
+            };
+        }
+        else if constexpr (Degree == 3)
+        {
+            // 2x2 Gauss-Legendre rule
+            constexpr Real g = std::numbers::inv_sqrt3_v<Real>; // 1/sqrt(3)
+            constexpr Real one = static_cast<Real>(1);
+            return std::array<PointAndWeight, 4>{
+                std::make_pair(Point(-g, -g), one),
+                std::make_pair(Point( g, -g), one),
+                std::make_pair(Point( g,  g), one),
+                std::make_pair(Point(-g,  g), one)
+            };
+        }
+        else
+        {
+            static_assert(degreeIsNotTabulated<Degree>, "QuadratureRules<Quad>: degree is not tabulated");
+        }
+    }
+};
+
+/// Reference domain: the unit tetrahedron (0,0,0), (1,0,0), (0,1,0), (0,0,1).
+template <>
+struct QuadratureRules<sofa::geometry::Tetrahedron>
+{
+    static constexpr sofa::Size TopologicalDimension = 3;
+    static constexpr std::array<sofa::Size, 2> tabulatedDegrees { 1, 2 };
+
+    template <sofa::Size Degree, class Real>
+    static constexpr auto points()
+    {
+        using Point = sofa::type::Vec<TopologicalDimension, Real>;
+        using PointAndWeight = std::pair<Point, Real>;
+
+        if constexpr (Degree == 1)
+        {
+            // 1-point centroid rule
+            constexpr Real quarter = static_cast<Real>(0.25); // 1/4
+            return std::array<PointAndWeight, 1>{
+                std::make_pair(Point(quarter, quarter, quarter),
+                               static_cast<Real>(1) / static_cast<Real>(6)) // volume of the reference tetrahedron
+            };
+        }
+        else if constexpr (Degree == 2)
+        {
+            // 4-point rule
+            constexpr Real sqrt5 = static_cast<Real>(2.23606797749978969640917366873); // sqrt(5)
+            constexpr Real a = (static_cast<Real>(5) - sqrt5) / static_cast<Real>(20);
+            constexpr Real b = (static_cast<Real>(5) + static_cast<Real>(3) * sqrt5) / static_cast<Real>(20);
+            constexpr Real w = static_cast<Real>(1) / static_cast<Real>(24); // 1/24
+
+            return std::array<PointAndWeight, 4>{
+                std::make_pair(Point(a, a, a), w),
+                std::make_pair(Point(b, a, a), w),
+                std::make_pair(Point(a, b, a), w),
+                std::make_pair(Point(a, a, b), w)
+            };
+        }
+        else
+        {
+            static_assert(degreeIsNotTabulated<Degree>, "QuadratureRules<Tetrahedron>: degree is not tabulated");
+        }
+    }
+};
+
+/// Reference domain: the cube [-1, 1]^3.
+template <>
+struct QuadratureRules<sofa::geometry::Hexahedron>
+{
+    static constexpr sofa::Size TopologicalDimension = 3;
+    static constexpr std::array<sofa::Size, 3> tabulatedDegrees { 1, 3, 5 };
+
+    template <sofa::Size Degree, class Real>
+    static constexpr auto points()
+    {
+        using Point = sofa::type::Vec<TopologicalDimension, Real>;
+        using PointAndWeight = std::pair<Point, Real>;
+
+        if constexpr (Degree == 1)
+        {
+            // 1-point centroid rule
+            constexpr Real zero = static_cast<Real>(0);
+            return std::array<PointAndWeight, 1>{
+                std::make_pair(Point(zero, zero, zero), static_cast<Real>(8))
+            };
+        }
+        else if constexpr (Degree == 3)
+        {
+            // 2x2x2 Gauss-Legendre rule, in the node order of sofa::geometry::Hexahedron
+            constexpr Real g = std::numbers::inv_sqrt3_v<Real>; // 1/sqrt(3)
+            constexpr Real one = static_cast<Real>(1);
+            return std::array<PointAndWeight, 8>{
+                std::make_pair(Point(-g, -g, -g), one),
+                std::make_pair(Point( g, -g, -g), one),
+                std::make_pair(Point( g,  g, -g), one),
+                std::make_pair(Point(-g,  g, -g), one),
+                std::make_pair(Point(-g, -g,  g), one),
+                std::make_pair(Point( g, -g,  g), one),
+                std::make_pair(Point( g,  g,  g), one),
+                std::make_pair(Point(-g,  g,  g), one)
+            };
+        }
+        else if constexpr (Degree == 5)
+        {
+            // 3x3x3 Gauss-Legendre rule
+            constexpr Real g = static_cast<Real>(0.774596669241483377035853079956); // sqrt(3/5)
+            constexpr std::array<Real, 3> node { -g, static_cast<Real>(0), g };
+            constexpr std::array<Real, 3> weight {
+                static_cast<Real>(5) / static_cast<Real>(9), // 5/9
+                static_cast<Real>(8) / static_cast<Real>(9), // 8/9
+                static_cast<Real>(5) / static_cast<Real>(9)  // 5/9
+            };
+
+            std::array<PointAndWeight, 27> q {};
+            sofa::Size k = 0;
+            for (sofa::Size i = 0; i < 3; ++i)
+                for (sofa::Size j = 0; j < 3; ++j)
+                    for (sofa::Size l = 0; l < 3; ++l)
+                        q[k++] = std::make_pair(Point(node[i], node[j], node[l]),
+                                                weight[i] * weight[j] * weight[l]);
+            return q;
+        }
+        else
+        {
+            static_assert(degreeIsNotTabulated<Degree>, "QuadratureRules<Hexahedron>: degree is not tabulated");
+        }
+    }
+};
+
+/// Reference domain: the unit triangle extruded along z over [0, 1].
+template <>
+struct QuadratureRules<sofa::geometry::Prism>
+{
+    static constexpr sofa::Size TopologicalDimension = 3;
+    static constexpr std::array<sofa::Size, 1> tabulatedDegrees { 1 };
+
+    template <sofa::Size Degree, class Real>
+    static constexpr auto points()
+    {
+        using Point = sofa::type::Vec<TopologicalDimension, Real>;
+        using PointAndWeight = std::pair<Point, Real>;
+
+        if constexpr (Degree == 1)
+        {
+            // centroid of the triangle x 2-point Gauss-Legendre along the extrusion
+            constexpr Real third = static_cast<Real>(1) / static_cast<Real>(3); // 1/3
+            constexpr Real g = std::numbers::inv_sqrt3_v<Real>; // 1/sqrt(3)
+            constexpr Real half = static_cast<Real>(0.5);
+            constexpr Real one = static_cast<Real>(1);
+            constexpr Real w = static_cast<Real>(0.25); // 1/4
+
+            return std::array<PointAndWeight, 2>{
+                std::make_pair(Point(third, third, half * (one - g)), w),
+                std::make_pair(Point(third, third, half * (one + g)), w)
+            };
+        }
+        else
+        {
+            static_assert(degreeIsNotTabulated<Degree>, "QuadratureRules<Prism>: degree is not tabulated");
+        }
+    }
+};
+
+/// Reference domain: the cube [-1, 1]^3 collapsed onto the pyramid apex.
+template <>
+struct QuadratureRules<sofa::geometry::Pyramid>
+{
+    static constexpr sofa::Size TopologicalDimension = 3;
+    static constexpr std::array<sofa::Size, 1> tabulatedDegrees { 1 };
+
+    template <sofa::Size Degree, class Real>
+    static constexpr auto points()
+    {
+        using Point = sofa::type::Vec<TopologicalDimension, Real>;
+        using PointAndWeight = std::pair<Point, Real>;
+
+        if constexpr (Degree == 1)
+        {
+            // 2x2x2 Gauss-Legendre rule
+            constexpr Real g = std::numbers::inv_sqrt3_v<Real>; // 1/sqrt(3)
+            constexpr Real one = static_cast<Real>(1);
+            return std::array<PointAndWeight, 8>{
+                std::make_pair(Point(-g, -g, -g), one),
+                std::make_pair(Point( g, -g, -g), one),
+                std::make_pair(Point( g,  g, -g), one),
+                std::make_pair(Point(-g,  g, -g), one),
+                std::make_pair(Point(-g, -g,  g), one),
+                std::make_pair(Point( g, -g,  g), one),
+                std::make_pair(Point( g,  g,  g), one),
+                std::make_pair(Point(-g,  g,  g), one)
+            };
+        }
+        else
+        {
+            static_assert(degreeIsNotTabulated<Degree>, "QuadratureRules<Pyramid>: degree is not tabulated");
+        }
+    }
+};
+
+/**
+ * Selects a quadrature rule out of QuadratureRules<Domain>, at compile time or at runtime.
+ *
+ * A requested degree that is not tabulated is served by the smallest tabulated degree
+ * above it, so a caller states the accuracy it needs and never a table size. The
+ * selection is derived from Rules::tabulatedDegrees, which stays the single place where
+ * the available rules of a domain are listed.
+ */
+template <class Domain>
+struct QuadratureDispatch
+{
+    using Rules = QuadratureRules<Domain>;
+    static constexpr sofa::Size TopologicalDimension = Rules::TopologicalDimension;
+
+    template <class Real>
+    using PointAndWeight = QuadraturePointAndWeight_t<Domain, Real>;
+
+    /// Smallest tabulated degree integrating a polynomial of the requested degree
+    /// exactly, or 0 when the request exceeds every table of this domain.
+    static constexpr sofa::Size smallestSufficientDegree(sofa::Size requestedDegree)
+    {
+        for (const sofa::Size degree : Rules::tabulatedDegrees)
+        {
+            if (degree >= requestedDegree)
+            {
+                return degree;
+            }
+        }
+        return 0;
+    }
+
+    /// The table itself, selected at compile time.
+    template <sofa::Size RequestedDegree, class Real>
+    static constexpr auto points()
+    {
+        constexpr sofa::Size degree = smallestSufficientDegree(RequestedDegree);
+        static_assert(degree != 0,
+            "QuadratureDispatch: no quadrature rule is tabulated for the requested degree");
+        return Rules::template points<degree, Real>();
+    }
+
+    /// View of the table matching a degree known only at runtime.
+    template <class Real>
+    static std::span<const PointAndWeight<Real>> rule(sofa::Size requestedDegree)
+    {
+        return select<Real>(requestedDegree,
+            std::make_index_sequence<Rules::tabulatedDegrees.size()>{});
+    }
+
+private:
+    /// The table of one tabulated degree.
+    template <sofa::Size Degree, class Real>
+    static std::span<const PointAndWeight<Real>> table()
+    {
+        static constexpr auto tabulated = Rules::template points<Degree, Real>();
+        return tabulated;
+    }
+
+    /**
+     * Runtime counterpart of smallestSufficientDegree: keeps the table of the first
+     * tabulated degree that covers the request.
+     *
+     * Index is the pack 0, 1, ... over Rules::tabulatedDegrees, so the fold below expands
+     * to one test per tabulated degree. For a domain tabulating {1, 3, 5} it reads:
+     *
+     *     (requestedDegree <= 1 ? (selected = table<1, Real>(), true) : false) ||
+     *     (requestedDegree <= 3 ? (selected = table<3, Real>(), true) : false) ||
+     *     (requestedDegree <= 5 ? (selected = table<5, Real>(), true) : false)
+     *
+     * tabulatedDegrees is ascending and || short-circuits, so the first sufficient table
+     * is the one kept and the smaller ones are never evaluated. Each term yields a bool
+     * because a fold over || demands one; the comma operator performs the assignment and
+     * then reports true. The fold is false only when no degree matched, hence the throw.
+     */
+    template <class Real, std::size_t... Index>
+    static std::span<const PointAndWeight<Real>> select(
+        sofa::Size requestedDegree, std::index_sequence<Index...>)
+    {
+        std::span<const PointAndWeight<Real>> selected;
+
+        const bool isTabulated = ((requestedDegree <= Rules::tabulatedDegrees[Index]
+            ? (selected = table<Rules::tabulatedDegrees[Index], Real>(), true)
+            : false) || ...);
+
+        if (!isTabulated)
+        {
+            throw std::invalid_argument(
+                std::string("QuadratureDispatch<")
+                + sofa::geometry::elementTypeToString(Domain::Element_type)
+                + ">: no quadrature rule is tabulated for degree "
+                + std::to_string(requestedDegree));
+        }
+
+        return selected;
+    }
+};
+
+/**
+ * Quadrature interface of a finite element: the rules of its reference domain, restricted
+ * to the degrees its interpolation order admits.
+ *
+ * ElementType supplies the domain (through ReferenceDomain) and the interpolation order
+ * (through sofa::geometry), so an element only states the degree it wants by default.
+ * DefaultDegree is given explicitly rather than derived, because the accuracy an element
+ * ships with is a choice, not a consequence of its order.
+ */
+template <class ElementType, class Real, sofa::Size DefaultDegree>
+struct FiniteElementQuadrature
+{
+    using Domain = ReferenceDomain_t<ElementType>;
+    using Dispatch = QuadratureDispatch<Domain>;
+    using QuadraturePointAndWeight = QuadraturePointAndWeight_t<Domain, Real>;
+
+    static constexpr sofa::Size InterpolationOrder = ElementType::PolynomialOrder;
+
+    /// Lowest degree integrating the stiffness integrand of an order-p element exactly on
+    /// an affine mapping: the product of two shape function gradients is of degree 2(p-1).
+    /// Clamped to 1, since no degree-0 rule is tabulated.
+    static constexpr sofa::Size MinimumQuadratureDegree =
+        InterpolationOrder > 1 ? 2 * (InterpolationOrder - 1) : 1;
+
+    static constexpr sofa::Size DefaultQuadratureDegree = DefaultDegree;
+    static_assert(DefaultQuadratureDegree >= MinimumQuadratureDegree,
+        "the default quadrature degree is too low for the order of this element");
+
+    /// Quadrature points and weights of the given degree, resolved at compile time.
+    template <sofa::Size Degree = DefaultQuadratureDegree>
+    static constexpr auto quadraturePoints()
+    {
+        static_assert(Degree >= MinimumQuadratureDegree,
+            "quadrature degree too low for the interpolation order of this element");
+        return Dispatch::template points<Degree, Real>();
+    }
+
+    /// Quadrature points and weights of a degree known only at runtime; view of the
+    /// compile-time table.
+    static std::span<const QuadraturePointAndWeight> quadratureRule(sofa::Size degree)
+    {
+        if (degree < MinimumQuadratureDegree)
+        {
+            throw std::invalid_argument(
+                std::string("FiniteElement<")
+                + sofa::geometry::elementTypeToString(ElementType::Element_type)
+                + ">::quadratureRule: degree " + std::to_string(degree)
+                + " is too low for an element of interpolation order "
+                + std::to_string(InterpolationOrder));
+        }
+        return Dispatch::template rule<Real>(degree);
+    }
+};
+
+}
diff --git a/Sofa/framework/FEM/test/FiniteElement_test.cpp b/Sofa/framework/FEM/test/FiniteElement_test.cpp
index 0b2100347c3..a4ad8a5c44f 100644
--- a/Sofa/framework/FEM/test/FiniteElement_test.cpp
+++ b/Sofa/framework/FEM/test/FiniteElement_test.cpp
@@ -26,76 +26,81 @@ namespace sofa
 {
 
 /**
- * Computes the sum of the quadrature weights and compare it to an expected value
+ * For each quadrature degree the element admits, computes the sum of the quadrature weights
+ * and compares it to an expected value. Degrees below MinimumQuadratureDegree are rejected
+ * by the element, so the sweep starts there rather than at 1.
  */
-template <class ElementType, class DataTypes>
+template <class ElementType, class DataTypes, sofa::Size MaxDegree>
 void testSumWeights(const sofa::Real_t<DataTypes> expected)
 {
     using FE = sofa::fem::FiniteElement<ElementType, DataTypes>;
+    static_assert(MaxDegree >= FE::MinimumQuadratureDegree);
 
-    SReal weightSum = 0;
-    for (const auto& [q, w] : FE::quadraturePoints())
+    for (sofa::Size degree = FE::MinimumQuadratureDegree; degree <= MaxDegree; ++degree)
     {
-        weightSum += w;
+        SReal weightSum = 0;
+        for (const auto& [q, w] : FE::quadratureRule(degree))
+        {
+            weightSum += w;
+        }
+        EXPECT_DOUBLE_EQ(weightSum, expected) << "degree " << degree;
     }
-
-    EXPECT_DOUBLE_EQ(weightSum, expected);
 }
 
 TEST(FiniteElement, edge1dWeights)
 {
-    testSumWeights<sofa::geometry::Edge, sofa::defaulttype::Vec1Types>(2);
+    testSumWeights<sofa::geometry::Edge, sofa::defaulttype::Vec1Types, 3>(2);
 }
 TEST(FiniteElement, edge2dWeights)
 {
-    testSumWeights<sofa::geometry::Edge, sofa::defaulttype::Vec2Types>(2);
+    testSumWeights<sofa::geometry::Edge, sofa::defaulttype::Vec2Types, 3>(2);
 }
 TEST(FiniteElement, edge3dWeights)
 {
-    testSumWeights<sofa::geometry::Edge, sofa::defaulttype::Vec3Types>(2);
+    testSumWeights<sofa::geometry::Edge, sofa::defaulttype::Vec3Types, 3>(2);
 }
 
 TEST(FiniteElement, triangle2dWeights)
 {
-    testSumWeights<sofa::geometry::Triangle, sofa::defaulttype::Vec2Types>(0.5);
+    testSumWeights<sofa::geometry::Triangle, sofa::defaulttype::Vec2Types, 2>(0.5);
 }
 TEST(FiniteElement, triangle3dWeights)
 {
-    testSumWeights<sofa::geometry::Triangle, sofa::defaulttype::Vec3Types>(0.5);
+    testSumWeights<sofa::geometry::Triangle, sofa::defaulttype::Vec3Types, 2>(0.5);
 }
 
 TEST(FiniteElement, quad2dWeights)
 {
-    testSumWeights<sofa::geometry::Quad, sofa::defaulttype::Vec2Types>(4);
+    testSumWeights<sofa::geometry::Quad, sofa::defaulttype::Vec2Types, 3>(4);
 }
 TEST(FiniteElement, quad3dWeights)
 {
-    testSumWeights<sofa::geometry::Quad, sofa::defaulttype::Vec3Types>(4);
+    testSumWeights<sofa::geometry::Quad, sofa::defaulttype::Vec3Types, 3>(4);
 }
 
 TEST(FiniteElement, tetra3dWeights)
 {
-    testSumWeights<sofa::geometry::Tetrahedron, sofa::defaulttype::Vec3Types>(1 / 6.);
+    testSumWeights<sofa::geometry::Tetrahedron, sofa::defaulttype::Vec3Types, 2>(1 / 6.);
 }
 
 TEST(FiniteElement, hexa3dWeights)
 {
-    testSumWeights<sofa::geometry::Hexahedron, sofa::defaulttype::Vec3Types>(8);
+    testSumWeights<sofa::geometry::Hexahedron, sofa::defaulttype::Vec3Types, 5>(8);
 }
 
 TEST(FiniteElement, prism3dWeights)
 {
-    testSumWeights<sofa::geometry::Prism, sofa::defaulttype::Vec3Types>(1 / 2.);
+    testSumWeights<sofa::geometry::Prism, sofa::defaulttype::Vec3Types, 1>(1 / 2.);
 }
 
 TEST(FiniteElement, pyramid3dWeights)
 {
-    testSumWeights<sofa::geometry::Pyramid, sofa::defaulttype::Vec3Types>(8);
+    testSumWeights<sofa::geometry::Pyramid, sofa::defaulttype::Vec3Types, 1>(8);
 }
 
 TEST(FiniteElement, quadraticHexa3dWeights)
 {
-    testSumWeights<sofa::geometry::QuadraticHexahedron, sofa::defaulttype::Vec3Types>(8);
+    testSumWeights<sofa::geometry::QuadraticHexahedron, sofa::defaulttype::Vec3Types, 5>(8);
 }
 
 /**
diff --git a/Sofa/framework/Geometry/src/sofa/geometry/Edge.h b/Sofa/framework/Geometry/src/sofa/geometry/Edge.h
index 645a408642a..4cc01b08017 100644
--- a/Sofa/framework/Geometry/src/sofa/geometry/Edge.h
+++ b/Sofa/framework/Geometry/src/sofa/geometry/Edge.h
@@ -36,6 +36,7 @@ namespace sofa::geometry
 struct Edge
 {
     static constexpr sofa::Size NumberOfNodes = 2;
+    static constexpr sofa::Size PolynomialOrder = 1;
     static constexpr ElementType Element_type = ElementType::EDGE;
 
     Edge() = delete;
diff --git a/Sofa/framework/Geometry/src/sofa/geometry/Hexahedron.h b/Sofa/framework/Geometry/src/sofa/geometry/Hexahedron.h
index 161f19ea304..b1996c6ad60 100644
--- a/Sofa/framework/Geometry/src/sofa/geometry/Hexahedron.h
+++ b/Sofa/framework/Geometry/src/sofa/geometry/Hexahedron.h
@@ -38,6 +38,7 @@ namespace sofa::geometry
 struct Hexahedron
 {
     static constexpr sofa::Size NumberOfNodes = 8;
+    static constexpr sofa::Size PolynomialOrder = 1;
     static constexpr ElementType Element_type = ElementType::HEXAHEDRON;
 
     Hexahedron() = delete;
diff --git a/Sofa/framework/Geometry/src/sofa/geometry/Prism.h b/Sofa/framework/Geometry/src/sofa/geometry/Prism.h
index 9ae91e534ab..a019ac44d6a 100644
--- a/Sofa/framework/Geometry/src/sofa/geometry/Prism.h
+++ b/Sofa/framework/Geometry/src/sofa/geometry/Prism.h
@@ -30,6 +30,7 @@ namespace sofa::geometry
 struct Prism
 {
     static constexpr sofa::Size NumberOfNodes = 6;
+    static constexpr sofa::Size PolynomialOrder = 1;
     static constexpr ElementType Element_type = ElementType::PRISM;
 
     Prism() = delete;
diff --git a/Sofa/framework/Geometry/src/sofa/geometry/Pyramid.h b/Sofa/framework/Geometry/src/sofa/geometry/Pyramid.h
index 7ac44849080..92b10fecb68 100644
--- a/Sofa/framework/Geometry/src/sofa/geometry/Pyramid.h
+++ b/Sofa/framework/Geometry/src/sofa/geometry/Pyramid.h
@@ -30,6 +30,7 @@ namespace sofa::geometry
 struct Pyramid
 {
     static constexpr sofa::Size NumberOfNodes = 5;
+    static constexpr sofa::Size PolynomialOrder = 1;
     static constexpr ElementType Element_type = ElementType::PYRAMID;
 
     Pyramid() = delete;
diff --git a/Sofa/framework/Geometry/src/sofa/geometry/Quad.h b/Sofa/framework/Geometry/src/sofa/geometry/Quad.h
index e091d63d6eb..4de1635c325 100644
--- a/Sofa/framework/Geometry/src/sofa/geometry/Quad.h
+++ b/Sofa/framework/Geometry/src/sofa/geometry/Quad.h
@@ -31,6 +31,7 @@ namespace sofa::geometry
 struct Quad
 {
     static constexpr sofa::Size NumberOfNodes = 4;
+    static constexpr sofa::Size PolynomialOrder = 1;
     static constexpr ElementType Element_type = ElementType::QUAD;
 
     Quad() = delete;
diff --git a/Sofa/framework/Geometry/src/sofa/geometry/Tetrahedron.h b/Sofa/framework/Geometry/src/sofa/geometry/Tetrahedron.h
index 156a89a4a92..7fed8a2d398 100644
--- a/Sofa/framework/Geometry/src/sofa/geometry/Tetrahedron.h
+++ b/Sofa/framework/Geometry/src/sofa/geometry/Tetrahedron.h
@@ -35,6 +35,7 @@ namespace sofa::geometry
 struct Tetrahedron
 {
     static constexpr sofa::Size NumberOfNodes = 4;
+    static constexpr sofa::Size PolynomialOrder = 1;
     static constexpr ElementType Element_type = ElementType::TETRAHEDRON;
 
     Tetrahedron() = delete;
diff --git a/Sofa/framework/Geometry/src/sofa/geometry/Triangle.h b/Sofa/framework/Geometry/src/sofa/geometry/Triangle.h
index 63e6c26fbbc..3707c2f861b 100644
--- a/Sofa/framework/Geometry/src/sofa/geometry/Triangle.h
+++ b/Sofa/framework/Geometry/src/sofa/geometry/Triangle.h
@@ -31,6 +31,7 @@ namespace sofa::geometry
 struct Triangle
 {
     static constexpr sofa::Size NumberOfNodes = 3;
+    static constexpr sofa::Size PolynomialOrder = 1;
     static constexpr ElementType Element_type = ElementType::TRIANGLE;
 
     Triangle() = delete;