diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h
index fbffb5dbfed..53570c30d05 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h
+++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h
@@ -21,6 +21,8 @@
 ******************************************************************************/
 #pragma once
 #include <sofa/fem/FiniteElement.h>
+#include <span>
+#include <stdexcept>
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_EDGE_CPP)
 #include <sofa/defaulttype/VecTypes.h>
@@ -55,12 +57,43 @@ struct FiniteElement<sofa::geometry::Edge, DataTypes>
         return {{-static_cast<Real>(0.5)}, {static_cast<Real>(0.5)}};
     }
 
-    static constexpr std::array<QuadraturePointAndWeight, 1> quadraturePoints()
+    template <sofa::Size Degree = 1>
+    static constexpr auto quadraturePoints()
     {
-        constexpr sofa::type::Vec<TopologicalDimension, Real> q0(static_cast<Real>(0));
-        return {
-            std::make_pair(q0, static_cast<Real>(2))
-        };
+        if constexpr (Degree <= 1)
+        {
+            // Degree 1: 1-point midpoint rule (default).
+            return std::array<QuadraturePointAndWeight, 1>{
+                std::make_pair(ReferenceCoord(static_cast<Real>(0)), static_cast<Real>(2))
+            };
+        }
+        else if constexpr (Degree <= 3)
+        {
+            // Degrees 2-3: 2-point Gauss-Legendre rule.
+            constexpr Real sqrt3 = 1.73205080757;
+            constexpr Real g = static_cast<Real>(1) / sqrt3;
+            return std::array<QuadraturePointAndWeight, 2>{
+                std::make_pair(ReferenceCoord(-g), static_cast<Real>(1)),
+                std::make_pair(ReferenceCoord( g), static_cast<Real>(1))
+            };
+        }
+        else
+        {
+            static_assert(Degree <= 3, "FiniteElement<Edge>: no quadrature rule for the requested degree");
+        }
+    }
+
+    // Quadrature rule selector by degree; view of the compile-time table.
+    static std::span<const QuadraturePointAndWeight> quadratureRule(sofa::Size degree)
+    {
+        switch (degree)
+        {
+            case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; }
+            case 2:
+            case 3: { static constexpr auto rule = quadraturePoints<3>(); return rule; }
+            default:
+                throw std::invalid_argument("FiniteElement<Edge>::quadratureRule: unsupported degree");
+        }
     }
 
 };
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h
index 38339624e47..d855c9f8c72 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h
+++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h
@@ -21,6 +21,8 @@
 ******************************************************************************/
 #pragma once
 #include <sofa/fem/FiniteElement.h>
+#include <span>
+#include <stdexcept>
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_HEXAHEDRON_CPP)
 #include <sofa/defaulttype/VecTypes.h>
@@ -96,24 +98,68 @@ struct FiniteElement<sofa::geometry::Hexahedron, DataTypes>
         return gradient;
     }
 
+    template <sofa::Size Degree = 3>
     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);
+        if constexpr (Degree <= 1)
+        {
+            // Degree 1: 1-point centroid rule.
+            return std::array<QuadraturePointAndWeight, 1>{
+                std::make_pair(ReferenceCoord(static_cast<Real>(0), static_cast<Real>(0), static_cast<Real>(0)), static_cast<Real>(8))
+            };
+        }
+        else if constexpr (Degree <= 3)
+        {
+            // Degrees 2-3: 2x2x2 Gauss-Legendre rule (default).
+            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 std::array {
+                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},
+            };
+        }
+        else if constexpr (Degree <= 5)
+        {
+            // Degrees 4-5: 3x3x3 Gauss-Legendre rule.
+            constexpr Real g = 0.77459666924; //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./9.), static_cast<Real>(8./9.), static_cast<Real>(5./9.) };
+
+            std::array<QuadraturePointAndWeight, 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(ReferenceCoord(node[i], node[j], node[l]), weight[i] * weight[j] * weight[l]);
+            return q;
+        }
+        else
+        {
+            static_assert(Degree <= 5, "FiniteElement<Hexahedron>: no quadrature rule for the requested degree");
+        }
+    }
 
-        return q;
+    // Quadrature rule selector by degree; view of the compile-time table.
+    static std::span<const QuadraturePointAndWeight> quadratureRule(sofa::Size degree)
+    {
+        switch (degree)
+        {
+            case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; }
+            case 2:
+            case 3: { static constexpr auto rule = quadraturePoints<3>(); return rule; }
+            case 4:
+            case 5: { static constexpr auto rule = quadraturePoints<5>(); return rule; }
+            default:
+                throw std::invalid_argument("FiniteElement<Hexahedron>::quadratureRule: unsupported degree");
+        }
     }
 };
 
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h
index 4e37585eace..261eeae7036 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h
+++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h
@@ -21,6 +21,8 @@
 ******************************************************************************/
 #pragma once
 #include <sofa/fem/FiniteElement.h>
+#include <span>
+#include <stdexcept>
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_QUAD_CPP)
 #include <sofa/defaulttype/VecTypes.h>
@@ -67,21 +69,62 @@ struct FiniteElement<sofa::geometry::Quad, DataTypes>
         };
     }
 
-    static constexpr std::array<QuadraturePointAndWeight, 3> quadraturePoints()
+    template <sofa::Size Degree = 2>
+    static constexpr auto quadraturePoints()
     {
-        constexpr Real sqrt2_3 = 0.816496580928; //sqrt(2./3.)
-        constexpr Real sqrt6 = 2.44948974278; //sqrt(6.)
-        constexpr Real sqrt2 = 1.41421356237; //sqrt(2.)
+        if constexpr (Degree <= 1)
+        {
+            // Degree 1: 1-point centroid rule.
+            return std::array<QuadraturePointAndWeight, 1>{
+                std::make_pair(ReferenceCoord(static_cast<Real>(0), static_cast<Real>(0)), static_cast<Real>(4))
+            };
+        }
+        else if constexpr (Degree <= 2)
+        {
+            // Degree 2: 3-point rule (default).
+            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);
+            constexpr ReferenceCoord q0(sqrt2_3, 0.);
+            constexpr ReferenceCoord q1(-1/sqrt6, -1./sqrt2);
+            constexpr ReferenceCoord q2(-1/sqrt6, 1./sqrt2);
 
-        return {
-            std::make_pair(q0, 4./3.),
-            std::make_pair(q1, 4./3.),
-            std::make_pair(q2, 4./3.),
-        };
+            return std::array<QuadraturePointAndWeight, 3>{
+                std::make_pair(q0, 4./3.),
+                std::make_pair(q1, 4./3.),
+                std::make_pair(q2, 4./3.)
+            };
+        }
+        else if constexpr (Degree <= 3)
+        {
+            // Degree 3: 2x2 Gauss-Legendre rule.
+            constexpr Real sqrt3 = 1.73205080757; //sqrt(3.)
+            constexpr Real g = static_cast<Real>(1) / sqrt3;
+            return std::array<QuadraturePointAndWeight, 4>{
+                std::make_pair(ReferenceCoord(-g, -g), static_cast<Real>(1)),
+                std::make_pair(ReferenceCoord( g, -g), static_cast<Real>(1)),
+                std::make_pair(ReferenceCoord( g,  g), static_cast<Real>(1)),
+                std::make_pair(ReferenceCoord(-g,  g), static_cast<Real>(1))
+            };
+        }
+        else
+        {
+            static_assert(Degree <= 3, "FiniteElement<Quad>: no quadrature rule for the requested degree");
+        }
+    }
+
+    // Quadrature rule selector by degree; view of the compile-time table.
+    static std::span<const QuadraturePointAndWeight> quadratureRule(sofa::Size degree)
+    {
+        switch (degree)
+        {
+            case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; }
+            case 2: { static constexpr auto rule = quadraturePoints<2>(); return rule; }
+            case 3: { static constexpr auto rule = quadraturePoints<3>(); return rule; }
+            default:
+                throw std::invalid_argument("FiniteElement<Quad>::quadratureRule: unsupported degree");
+        }
     }
 };
 
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h
index ac14bb9c59b..bca61671e0a 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h
+++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h
@@ -21,6 +21,8 @@
 ******************************************************************************/
 #pragma once
 #include <sofa/fem/FiniteElement.h>
+#include <span>
+#include <stdexcept>
 
 #if !defined(SOFA_FEM_FINITE_ELEMENT_TETAHEDRON_CPP)
 #include <sofa/defaulttype/VecTypes.h>
@@ -68,11 +70,45 @@ struct FiniteElement<sofa::geometry::Tetrahedron, DataTypes>
         };
     }
 
-    static constexpr std::array<QuadraturePointAndWeight, 1> quadraturePoints()
+    template <sofa::Size Degree = 1>
+    static constexpr auto 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 constexpr (Degree <= 1)
+        {
+            // Degree 1: 1-point centroid rule (default).
+            return std::array<QuadraturePointAndWeight, 1>{
+                std::make_pair(ReferenceCoord(1./4., 1./4., 1./4.), Real(1./6.))
+            };
+        }
+        else if constexpr (Degree <= 2)
+        {
+            // Degree 2: 4-point rule.
+            constexpr Real sqrt5 = 2.2360679774997896;
+            constexpr Real a = (5. - sqrt5) / 20.;
+            constexpr Real b = (5. + 3. * sqrt5) / 20.;
+            return std::array<QuadraturePointAndWeight, 4>{
+                std::make_pair(ReferenceCoord(a, a, a), Real(1./24.)),
+                std::make_pair(ReferenceCoord(b, a, a), Real(1./24.)),
+                std::make_pair(ReferenceCoord(a, b, a), Real(1./24.)),
+                std::make_pair(ReferenceCoord(a, a, b), Real(1./24.))
+            };
+        }
+        else
+        {
+            static_assert(Degree <= 2, "FiniteElement<Tetrahedron>: no quadrature rule for the requested degree");
+        }
+    }
+
+    // Quadrature rule selector by degree; view of the compile-time table.
+    static std::span<const QuadraturePointAndWeight> quadratureRule(sofa::Size degree)
+    {
+        switch (degree)
+        {
+            case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; }
+            case 2: { static constexpr auto rule = quadraturePoints<2>(); return rule; }
+            default:
+                throw std::invalid_argument("FiniteElement<Tetrahedron>::quadratureRule: unsupported degree");
+        }
     }
 };
 
diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h
index 0d9659040de..84c9ec19a74 100644
--- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h
+++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h
@@ -21,6 +21,8 @@
 ******************************************************************************/
 #pragma once
 #include <sofa/fem/FiniteElement.h>
+#include <span>
+#include <stdexcept>
 
 namespace sofa::fem
 {
@@ -64,11 +66,41 @@ struct FiniteElement<sofa::geometry::Triangle, DataTypes>
         };
     }
 
-    static constexpr std::array<QuadraturePointAndWeight, 1> quadraturePoints()
+    template <sofa::Size Degree = 1>
+    static constexpr auto quadraturePoints()
     {
-        return {
-            std::make_pair(sofa::type::Vec<TopologicalDimension, Real>(1./3., 1./3.), 1./2.)
-        };
+        if constexpr (Degree <= 1)
+        {
+            // Degree 1: 1-point centroid rule (default).
+            return std::array<QuadraturePointAndWeight, 1>{
+                std::make_pair(ReferenceCoord(1./3., 1./3.), Real(1./2.))
+            };
+        }
+        else if constexpr (Degree <= 2)
+        {
+            // Degree 2: 3-point interior rule.
+            return std::array<QuadraturePointAndWeight, 3>{
+                std::make_pair(ReferenceCoord(1./6., 1./6.), Real(1./6.)),
+                std::make_pair(ReferenceCoord(2./3., 1./6.), Real(1./6.)),
+                std::make_pair(ReferenceCoord(1./6., 2./3.), Real(1./6.))
+            };
+        }
+        else
+        {
+            static_assert(Degree <= 2, "FiniteElement<Triangle>: no quadrature rule for the requested degree");
+        }
+    }
+
+    // Quadrature rule selector by degree; view of the compile-time table.
+    static std::span<const QuadraturePointAndWeight> quadratureRule(sofa::Size degree)
+    {
+        switch (degree)
+        {
+            case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; }
+            case 2: { static constexpr auto rule = quadraturePoints<2>(); return rule; }
+            default:
+                throw std::invalid_argument("FiniteElement<Triangle>::quadratureRule: unsupported degree");
+        }
     }
 };