From 1dc7ec095a9229abf5c3b3498f79f4174f0c90ce Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:07:19 +0200 Subject: [PATCH 1/6] [FEM] Introduce higher order quadrature rule switch for Triangle --- .../src/sofa/fem/FiniteElement[Triangle].h | 40 +++++++++++++++++-- 1 file changed, 36 insertions(+), 4 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h index 0d9659040de..7495e802b26 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 +#include +#include namespace sofa::fem { @@ -64,11 +66,41 @@ struct FiniteElement }; } - static constexpr std::array quadraturePoints() + template + static constexpr auto quadraturePoints() { - return { - std::make_pair(sofa::type::Vec(1./3., 1./3.), 1./2.) - }; + if constexpr (Degree <= 1) + { + // Degree 1: 1-point centroid rule (default). + return std::array{ + 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{ + 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: no quadrature rule for the requested degree"); + } + } + + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(int 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::quadratureRule: unsupported degree"); + } } }; From 9ec52427d04fb8504b4dc4de62a24c214ed978df Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:17:19 +0200 Subject: [PATCH 2/6] Edge --- .../FEM/src/sofa/fem/FiniteElement[Edge].h | 43 ++++++++++++++++--- 1 file changed, 38 insertions(+), 5 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h index fbffb5dbfed..492caa7a380 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 +#include +#include #if !defined(SOFA_FEM_FINITE_ELEMENT_EDGE_CPP) #include @@ -55,12 +57,43 @@ struct FiniteElement return {{-static_cast(0.5)}, {static_cast(0.5)}}; } - static constexpr std::array quadraturePoints() + template + static constexpr auto quadraturePoints() { - constexpr sofa::type::Vec q0(static_cast(0)); - return { - std::make_pair(q0, static_cast(2)) - }; + if constexpr (Degree <= 1) + { + // Degree 1: 1-point midpoint rule (default). + return std::array{ + std::make_pair(ReferenceCoord(static_cast(0)), static_cast(2)) + }; + } + else if constexpr (Degree <= 3) + { + // Degrees 2-3: 2-point Gauss-Legendre rule. + constexpr Real sqrt3 = 1.73205080757; + constexpr Real g = static_cast(1) / sqrt3; + return std::array{ + std::make_pair(ReferenceCoord(-g), static_cast(1)), + std::make_pair(ReferenceCoord( g), static_cast(1)) + }; + } + else + { + static_assert(Degree <= 3, "FiniteElement: no quadrature rule for the requested degree"); + } + } + + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(int 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::quadratureRule: unsupported degree"); + } } }; From 213ccaf7494313ec35109832996c5f19465b6c3d Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:23:24 +0200 Subject: [PATCH 3/6] Tetrahedra --- .../src/sofa/fem/FiniteElement[Tetrahedron].h | 44 +++++++++++++++++-- 1 file changed, 40 insertions(+), 4 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h index ac14bb9c59b..4ab10ae3ed5 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 +#include +#include #if !defined(SOFA_FEM_FINITE_ELEMENT_TETAHEDRON_CPP) #include @@ -68,11 +70,45 @@ struct FiniteElement }; } - static constexpr std::array quadraturePoints() + template + static constexpr auto quadraturePoints() { - constexpr sofa::type::Vec q0(1./4., 1./4., 1./4.); - constexpr std::array q { std::make_pair(q0, 1./6.) }; - return q; + if constexpr (Degree <= 1) + { + // Degree 1: 1-point centroid rule (default). + return std::array{ + 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{ + 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: no quadrature rule for the requested degree"); + } + } + + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(int 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::quadratureRule: unsupported degree"); + } } }; From cfb35e235ebee1bb3253760fe3f3bd08eb8ce3f7 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:41:34 +0200 Subject: [PATCH 4/6] Quad --- .../FEM/src/sofa/fem/FiniteElement[Quad].h | 67 +++++++++++++++---- 1 file changed, 55 insertions(+), 12 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h index 4e37585eace..1485a881001 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 +#include +#include #if !defined(SOFA_FEM_FINITE_ELEMENT_QUAD_CPP) #include @@ -67,21 +69,62 @@ struct FiniteElement }; } - static constexpr std::array quadraturePoints() + template + 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{ + std::make_pair(ReferenceCoord(static_cast(0), static_cast(0)), static_cast(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 q0(sqrt2_3, 0.); - constexpr sofa::type::Vec q1(-1/sqrt6, -1./sqrt2); - constexpr sofa::type::Vec 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{ + 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(1) / sqrt3; + return std::array{ + std::make_pair(ReferenceCoord(-g, -g), static_cast(1)), + std::make_pair(ReferenceCoord( g, -g), static_cast(1)), + std::make_pair(ReferenceCoord( g, g), static_cast(1)), + std::make_pair(ReferenceCoord(-g, g), static_cast(1)) + }; + } + else + { + static_assert(Degree <= 3, "FiniteElement: no quadrature rule for the requested degree"); + } + } + + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(int 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::quadratureRule: unsupported degree"); + } } }; From 5ad33f194ee045f847056162a61f889a9ed2e9d0 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:41:41 +0200 Subject: [PATCH 5/6] Hexa --- .../src/sofa/fem/FiniteElement[Hexahedron].h | 74 +++++++++++++++---- 1 file changed, 60 insertions(+), 14 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h index 38339624e47..e3ab476e0f8 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 +#include +#include #if !defined(SOFA_FEM_FINITE_ELEMENT_HEXAHEDRON_CPP) #include @@ -96,24 +98,68 @@ struct FiniteElement return gradient; } + template static constexpr auto quadraturePoints() { - constexpr Real sqrt3 = 1.73205080757; //sqrt(3.) - constexpr Real sqrt3_1 = static_cast(1) / sqrt3; - constexpr Real one = static_cast(1); + if constexpr (Degree <= 1) + { + // Degree 1: 1-point centroid rule. + return std::array{ + std::make_pair(ReferenceCoord(static_cast(0), static_cast(0), static_cast(0)), static_cast(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(1) / sqrt3; + constexpr Real one = static_cast(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 node{ -g, static_cast(0), g }; + constexpr std::array weight{ static_cast(5./9.), static_cast(8./9.), static_cast(5./9.) }; + + std::array 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: no quadrature rule for the requested degree"); + } + } - return q; + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(int 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::quadratureRule: unsupported degree"); + } } }; From 4b6a9dc9c719ca76b48cfc87aea3a87a8ffad724 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 11 Aug 2026 15:28:18 +0200 Subject: [PATCH 6/6] Use uint for degree --- Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h | 4 ++-- Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h | 4 ++-- Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h | 4 ++-- Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h | 4 ++-- Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h | 4 ++-- 5 files changed, 10 insertions(+), 10 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h index 492caa7a380..53570c30d05 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h @@ -57,7 +57,7 @@ struct FiniteElement return {{-static_cast(0.5)}, {static_cast(0.5)}}; } - template + template static constexpr auto quadraturePoints() { if constexpr (Degree <= 1) @@ -84,7 +84,7 @@ struct FiniteElement } // Quadrature rule selector by degree; view of the compile-time table. - static std::span quadratureRule(int degree) + static std::span quadratureRule(sofa::Size degree) { switch (degree) { diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h index e3ab476e0f8..d855c9f8c72 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h @@ -98,7 +98,7 @@ struct FiniteElement return gradient; } - template + template static constexpr auto quadraturePoints() { if constexpr (Degree <= 1) @@ -148,7 +148,7 @@ struct FiniteElement } // Quadrature rule selector by degree; view of the compile-time table. - static std::span quadratureRule(int degree) + static std::span quadratureRule(sofa::Size degree) { switch (degree) { diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h index 1485a881001..261eeae7036 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h @@ -69,7 +69,7 @@ struct FiniteElement }; } - template + template static constexpr auto quadraturePoints() { if constexpr (Degree <= 1) @@ -115,7 +115,7 @@ struct FiniteElement } // Quadrature rule selector by degree; view of the compile-time table. - static std::span quadratureRule(int degree) + static std::span quadratureRule(sofa::Size degree) { switch (degree) { diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h index 4ab10ae3ed5..bca61671e0a 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h @@ -70,7 +70,7 @@ struct FiniteElement }; } - template + template static constexpr auto quadraturePoints() { if constexpr (Degree <= 1) @@ -100,7 +100,7 @@ struct FiniteElement } // Quadrature rule selector by degree; view of the compile-time table. - static std::span quadratureRule(int degree) + static std::span quadratureRule(sofa::Size degree) { switch (degree) { diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h index 7495e802b26..84c9ec19a74 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h @@ -66,7 +66,7 @@ struct FiniteElement }; } - template + template static constexpr auto quadraturePoints() { if constexpr (Degree <= 1) @@ -92,7 +92,7 @@ struct FiniteElement } // Quadrature rule selector by degree; view of the compile-time table. - static std::span quadratureRule(int degree) + static std::span quadratureRule(sofa::Size degree) { switch (degree) {