From f70bd97e4a871a388457b1e73100a80d4f6b8d3f Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:07:19 +0200 Subject: [PATCH 01/18] [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 542d1474aa79f4f3f6e9de90e7d7ebddbf12e5e3 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:17:19 +0200 Subject: [PATCH 02/18] 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 2f86678dc0348837e04b2b5689d6f838b2d723f0 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:23:24 +0200 Subject: [PATCH 03/18] 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 e0a5955911455b5f062c7d7e31c33dda36bc15d0 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:41:34 +0200 Subject: [PATCH 04/18] 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 39d29429cbe3bdbe54f4cbddc0b541b3ba8cd2a2 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 10 Aug 2026 17:41:41 +0200 Subject: [PATCH 05/18] 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 4abd0a26e271296c701b6c6b3255a8fb8470d464 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Tue, 11 Aug 2026 15:28:18 +0200 Subject: [PATCH 06/18] 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) { From 85b9f7907eb62ff2e5c2b73eeb223a3119dcf21f Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Fri, 4 Sep 2026 14:17:14 +0200 Subject: [PATCH 07/18] Adopt quadratureRule in prism and pyramid --- .../FEM/src/sofa/fem/FiniteElement[Prism].h | 44 ++++++++++++----- .../FEM/src/sofa/fem/FiniteElement[Pyramid].h | 48 ++++++++++++++----- 2 files changed, 68 insertions(+), 24 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h index 5a6731ff94c..8951dc385f3 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h @@ -21,6 +21,8 @@ ******************************************************************************/ #pragma once #include +#include +#include #if !defined(SOFA_FEM_FINITE_ELEMENT_PRISM_CPP) #include @@ -74,19 +76,39 @@ struct FiniteElement }; } - static constexpr std::array quadraturePoints() + template + static constexpr auto quadraturePoints() { - constexpr auto third = static_cast(1) / static_cast(3); - constexpr auto sqrt_3 = static_cast(0.57735026919); // 1/sqrt(3) - constexpr auto one = static_cast(1); - constexpr QuadraturePoint q0 {third, third, static_cast(0.5) * (one - sqrt_3)}; - constexpr QuadraturePoint q1 {third, third, static_cast(0.5) * (one + sqrt_3)}; + if constexpr (Degree <= 1) + { + // Degree 1: 2-point rule (default). + constexpr auto third = static_cast(1) / static_cast(3); + constexpr auto sqrt_3 = static_cast(0.57735026919); // 1/sqrt(3) + constexpr auto one = static_cast(1); + constexpr QuadraturePoint q0 {third, third, static_cast(0.5) * (one - sqrt_3)}; + constexpr QuadraturePoint q1 {third, third, static_cast(0.5) * (one + sqrt_3)}; - constexpr std::array q { - std::make_pair(q0, 1./4.), - std::make_pair(q1, 1./4.), - }; - return q; + constexpr std::array q { + std::make_pair(q0, 1./4.), + std::make_pair(q1, 1./4.), + }; + return q; + } + else + { + static_assert(Degree <= 1, "FiniteElement: no quadrature rule for the requested degree"); + } + } + + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(sofa::Size degree) + { + switch (degree) + { + case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; } + default: + throw std::invalid_argument("FiniteElement::quadratureRule: unsupported degree"); + } } }; diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h index 2efcb1613ab..422fcb920f3 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h @@ -22,6 +22,8 @@ #pragma once #include #include +#include +#include #if !defined(SOFA_FEM_FINITE_ELEMENT_PYRAMID_CPP) #include @@ -71,22 +73,42 @@ struct FiniteElement }; } + template static constexpr auto quadraturePoints() { - constexpr Real sqrt3_1 = static_cast(1) / static_cast(1.73205080757); - constexpr Real one = static_cast(1); + if constexpr (Degree <= 1) + { + // Degree 1: 8-point rule (default). + constexpr Real sqrt3_1 = static_cast(1) / static_cast(1.73205080757); + constexpr Real one = static_cast(1); - // We use the 8 Gauss points of the hexahedron, which exactly integrate the (1-z)^2 Jacobian of the pyramid. - return std::array { - std::pair{ReferenceCoord{-sqrt3_1, -sqrt3_1, -sqrt3_1}, one}, - std::pair{ReferenceCoord{ sqrt3_1, -sqrt3_1, -sqrt3_1}, one}, - std::pair{ReferenceCoord{ sqrt3_1, sqrt3_1, -sqrt3_1}, one}, - std::pair{ReferenceCoord{-sqrt3_1, sqrt3_1, -sqrt3_1}, one}, - std::pair{ReferenceCoord{-sqrt3_1, -sqrt3_1, sqrt3_1}, one}, - std::pair{ReferenceCoord{ sqrt3_1, -sqrt3_1, sqrt3_1}, one}, - std::pair{ReferenceCoord{ sqrt3_1, sqrt3_1, sqrt3_1}, one}, - std::pair{ReferenceCoord{-sqrt3_1, sqrt3_1, sqrt3_1}, one}, - }; + constexpr std::array q { + 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}, + }; + return q; + } + else + { + static_assert(Degree <= 1, "FiniteElement: no quadrature rule for the requested degree"); + } + } + + // Quadrature rule selector by degree; view of the compile-time table. + static std::span quadratureRule(sofa::Size degree) + { + switch (degree) + { + case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; } + default: + throw std::invalid_argument("FiniteElement::quadratureRule: unsupported degree"); + } } }; From f6b83fcc083a30797870f76d77b45b51337c5b42 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Fri, 4 Sep 2026 14:17:32 +0200 Subject: [PATCH 08/18] Include quadratureRule in FE test --- .../framework/FEM/test/FiniteElement_test.cpp | 39 ++++++++++--------- 1 file changed, 21 insertions(+), 18 deletions(-) diff --git a/Sofa/framework/FEM/test/FiniteElement_test.cpp b/Sofa/framework/FEM/test/FiniteElement_test.cpp index 0b2100347c3..c607d06fe15 100644 --- a/Sofa/framework/FEM/test/FiniteElement_test.cpp +++ b/Sofa/framework/FEM/test/FiniteElement_test.cpp @@ -26,71 +26,74 @@ namespace sofa { /** - * Computes the sum of the quadrature weights and compare it to an expected value + * For each quadrature degree in [1, MaxDegree], computes the sum of the quadrature weights and + * compare it to an expected value */ -template +template void testSumWeights(const sofa::Real_t expected) { using FE = sofa::fem::FiniteElement; - SReal weightSum = 0; - for (const auto& [q, w] : FE::quadraturePoints()) + for (sofa::Size degree = 1; 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(2); + testSumWeights(2); } TEST(FiniteElement, edge2dWeights) { - testSumWeights(2); + testSumWeights(2); } TEST(FiniteElement, edge3dWeights) { - testSumWeights(2); + testSumWeights(2); } TEST(FiniteElement, triangle2dWeights) { - testSumWeights(0.5); + testSumWeights(0.5); } TEST(FiniteElement, triangle3dWeights) { - testSumWeights(0.5); + testSumWeights(0.5); } TEST(FiniteElement, quad2dWeights) { - testSumWeights(4); + testSumWeights(4); } TEST(FiniteElement, quad3dWeights) { - testSumWeights(4); + testSumWeights(4); } TEST(FiniteElement, tetra3dWeights) { - testSumWeights(1 / 6.); + testSumWeights(1 / 6.); } TEST(FiniteElement, hexa3dWeights) { - testSumWeights(8); + testSumWeights(8); } TEST(FiniteElement, prism3dWeights) { - testSumWeights(1 / 2.); + testSumWeights(1 / 2.); } TEST(FiniteElement, pyramid3dWeights) { - testSumWeights(8); + testSumWeights(8); } TEST(FiniteElement, quadraticHexa3dWeights) From b6569cde651379d99c95e20744f0d32f67656477 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Fri, 11 Sep 2026 14:54:37 +0200 Subject: [PATCH 09/18] Add polynomial order to linear elements --- Sofa/framework/Geometry/src/sofa/geometry/Edge.h | 1 + Sofa/framework/Geometry/src/sofa/geometry/Hexahedron.h | 1 + Sofa/framework/Geometry/src/sofa/geometry/Prism.h | 1 + Sofa/framework/Geometry/src/sofa/geometry/Pyramid.h | 1 + Sofa/framework/Geometry/src/sofa/geometry/Quad.h | 1 + Sofa/framework/Geometry/src/sofa/geometry/Tetrahedron.h | 1 + Sofa/framework/Geometry/src/sofa/geometry/Triangle.h | 1 + 7 files changed, 7 insertions(+) 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; From fe026d9964dc186fbbf7a319d4d1611b5538fe0e Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Fri, 11 Sep 2026 15:33:58 +0200 Subject: [PATCH 10/18] [FEM] Add QuadratureRules.h --- Sofa/framework/FEM/CMakeLists.txt | 1 + .../FEM/src/sofa/fem/QuadratureRules.h | 489 ++++++++++++++++++ 2 files changed, 490 insertions(+) create mode 100644 Sofa/framework/FEM/src/sofa/fem/QuadratureRules.h 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/QuadratureRules.h b/Sofa/framework/FEM/src/sofa/fem/QuadratureRules.h new file mode 100644 index 00000000000..15f10e443f0 --- /dev/null +++ b/Sofa/framework/FEM/src/sofa/fem/QuadratureRules.h @@ -0,0 +1,489 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace sofa::fem +{ + +/// Mapping element types to their integration reference domain. +/// This domain is the same irrespective of the order of the element. +template +struct ReferenceDomain +{ + using type = ElementType; +}; + +template <> struct ReferenceDomain { using type = sofa::geometry::Edge; }; +template <> struct ReferenceDomain { using type = sofa::geometry::Triangle; }; +template <> struct ReferenceDomain { using type = sofa::geometry::Quad; }; +template <> struct ReferenceDomain { using type = sofa::geometry::Tetrahedron; }; +template <> struct ReferenceDomain { using type = sofa::geometry::Hexahedron; }; + +template +using ReferenceDomain_t = typename ReferenceDomain::type; + +/// Always false, but dependent on its argument: will trigger compilation failure in case no rule +/// has been tabulated for the requested degree. +template +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(), 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 +struct QuadratureRules; + +template +using QuadraturePointAndWeight_t = + std::pair::TopologicalDimension, Real>, Real>; + +/// Reference domain: the segment [-1, 1]. +template <> +struct QuadratureRules +{ + static constexpr sofa::Size TopologicalDimension = 1; + static constexpr std::array tabulatedDegrees { 1, 3 }; + + template + static constexpr auto points() + { + using Point = sofa::type::Vec; + using PointAndWeight = std::pair; + + if constexpr (Degree == 1) + { + // 1-point midpoint rule + return std::array{ + std::make_pair(Point(static_cast(0)), static_cast(2)) + }; + } + else if constexpr (Degree == 3) + { + // 2-point Gauss-Legendre rule + constexpr Real g = std::numbers::inv_sqrt3_v; // 1/sqrt(3) + return std::array{ + std::make_pair(Point(-g), static_cast(1)), + std::make_pair(Point( g), static_cast(1)) + }; + } + else + { + static_assert(degreeIsNotTabulated, "QuadratureRules: degree is not tabulated"); + } + } +}; + +/// Reference domain: the unit triangle (0,0), (1,0), (0,1). +template <> +struct QuadratureRules +{ + static constexpr sofa::Size TopologicalDimension = 2; + static constexpr std::array tabulatedDegrees { 1, 2 }; + + template + static constexpr auto points() + { + using Point = sofa::type::Vec; + using PointAndWeight = std::pair; + + constexpr Real third = static_cast(1) / static_cast(3); // 1/3 + constexpr Real sixth = static_cast(1) / static_cast(6); // 1/6 + + if constexpr (Degree == 1) + { + // 1-point centroid rule + return std::array{ + std::make_pair(Point(third, third), static_cast(0.5)) // area of the reference triangle + }; + } + else if constexpr (Degree == 2) + { + // 3-point interior rule + constexpr Real a = static_cast(2) / static_cast(3); // 2/3 + return std::array{ + 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, "QuadratureRules: degree is not tabulated"); + } + } +}; + +/// Reference domain: the square [-1, 1]^2. +template <> +struct QuadratureRules +{ + static constexpr sofa::Size TopologicalDimension = 2; + static constexpr std::array tabulatedDegrees { 1, 2, 3 }; + + template + static constexpr auto points() + { + using Point = sofa::type::Vec; + using PointAndWeight = std::pair; + + if constexpr (Degree == 1) + { + // 1-point centroid rule + return std::array{ + std::make_pair(Point(static_cast(0), static_cast(0)), static_cast(4)) + }; + } + else if constexpr (Degree == 2) + { + // 3-point rule + constexpr Real sqrt2_3 = std::numbers::sqrt2_v * std::numbers::inv_sqrt3_v; // sqrt(2/3) + constexpr Real inv_sqrt6 = std::numbers::inv_sqrt3_v / std::numbers::sqrt2_v; // 1/sqrt(6) + constexpr Real inv_sqrt2 = std::numbers::sqrt2_v / static_cast(2); // 1/sqrt(2) + constexpr Real w = static_cast(4) / static_cast(3); // 4/3 + + return std::array{ + std::make_pair(Point( sqrt2_3, static_cast(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; // 1/sqrt(3) + constexpr Real one = static_cast(1); + return std::array{ + 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, "QuadratureRules: degree is not tabulated"); + } + } +}; + +/// Reference domain: the unit tetrahedron (0,0,0), (1,0,0), (0,1,0), (0,0,1). +template <> +struct QuadratureRules +{ + static constexpr sofa::Size TopologicalDimension = 3; + static constexpr std::array tabulatedDegrees { 1, 2 }; + + template + static constexpr auto points() + { + using Point = sofa::type::Vec; + using PointAndWeight = std::pair; + + if constexpr (Degree == 1) + { + // 1-point centroid rule + constexpr Real quarter = static_cast(0.25); // 1/4 + return std::array{ + std::make_pair(Point(quarter, quarter, quarter), + static_cast(1) / static_cast(6)) // volume of the reference tetrahedron + }; + } + else if constexpr (Degree == 2) + { + // 4-point rule + constexpr Real sqrt5 = static_cast(2.23606797749978969640917366873); // sqrt(5) + constexpr Real a = (static_cast(5) - sqrt5) / static_cast(20); + constexpr Real b = (static_cast(5) + static_cast(3) * sqrt5) / static_cast(20); + constexpr Real w = static_cast(1) / static_cast(24); // 1/24 + + return std::array{ + 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, "QuadratureRules: degree is not tabulated"); + } + } +}; + +/// Reference domain: the cube [-1, 1]^3. +template <> +struct QuadratureRules +{ + static constexpr sofa::Size TopologicalDimension = 3; + static constexpr std::array tabulatedDegrees { 1, 3, 5 }; + + template + static constexpr auto points() + { + using Point = sofa::type::Vec; + using PointAndWeight = std::pair; + + if constexpr (Degree == 1) + { + // 1-point centroid rule + constexpr Real zero = static_cast(0); + return std::array{ + std::make_pair(Point(zero, zero, zero), static_cast(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; // 1/sqrt(3) + constexpr Real one = static_cast(1); + return std::array{ + 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(0.774596669241483377035853079956); // sqrt(3/5) + constexpr std::array node { -g, static_cast(0), g }; + constexpr std::array weight { + static_cast(5) / static_cast(9), // 5/9 + static_cast(8) / static_cast(9), // 8/9 + static_cast(5) / static_cast(9) // 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(Point(node[i], node[j], node[l]), + weight[i] * weight[j] * weight[l]); + return q; + } + else + { + static_assert(degreeIsNotTabulated, "QuadratureRules: degree is not tabulated"); + } + } +}; + +/// Reference domain: the unit triangle extruded along z over [0, 1]. +template <> +struct QuadratureRules +{ + static constexpr sofa::Size TopologicalDimension = 3; + static constexpr std::array tabulatedDegrees { 1 }; + + template + static constexpr auto points() + { + using Point = sofa::type::Vec; + using PointAndWeight = std::pair; + + if constexpr (Degree == 1) + { + // centroid of the triangle x 2-point Gauss-Legendre along the extrusion + constexpr Real third = static_cast(1) / static_cast(3); // 1/3 + constexpr Real g = std::numbers::inv_sqrt3_v; // 1/sqrt(3) + constexpr Real half = static_cast(0.5); + constexpr Real one = static_cast(1); + constexpr Real w = static_cast(0.25); // 1/4 + + return std::array{ + std::make_pair(Point(third, third, half * (one - g)), w), + std::make_pair(Point(third, third, half * (one + g)), w) + }; + } + else + { + static_assert(degreeIsNotTabulated, "QuadratureRules: degree is not tabulated"); + } + } +}; + +/// Reference domain: the cube [-1, 1]^3 collapsed onto the pyramid apex. +template <> +struct QuadratureRules +{ + static constexpr sofa::Size TopologicalDimension = 3; + static constexpr std::array tabulatedDegrees { 1 }; + + template + static constexpr auto points() + { + using Point = sofa::type::Vec; + using PointAndWeight = std::pair; + + if constexpr (Degree == 1) + { + // 2x2x2 Gauss-Legendre rule + constexpr Real g = std::numbers::inv_sqrt3_v; // 1/sqrt(3) + constexpr Real one = static_cast(1); + return std::array{ + 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, "QuadratureRules: degree is not tabulated"); + } + } +}; + +/** + * Selects a quadrature rule out of QuadratureRules, 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 +struct QuadratureDispatch +{ + using Rules = QuadratureRules; + static constexpr sofa::Size TopologicalDimension = Rules::TopologicalDimension; + + template + using PointAndWeight = QuadraturePointAndWeight_t; + + /// 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 + 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(); + } + + /// View of the table matching a degree known only at runtime. + template + static std::span> rule(sofa::Size requestedDegree) + { + return select(requestedDegree, + std::make_index_sequence{}); + } + +private: + /// The table of one tabulated degree. + template + static std::span> table() + { + static constexpr auto tabulated = Rules::template points(); + 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 + static std::span> select( + sofa::Size requestedDegree, std::index_sequence) + { + std::span> selected; + + const bool isTabulated = ((requestedDegree <= Rules::tabulatedDegrees[Index] + ? (selected = table(), 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; + } +}; + +} From 1681cc5247f0d47e6fefaf6df1c9340cec267943 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Fri, 11 Sep 2026 16:17:59 +0200 Subject: [PATCH 11/18] Add the ElementQuadrature class --- .../FEM/src/sofa/fem/QuadratureRules.h | 54 +++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/Sofa/framework/FEM/src/sofa/fem/QuadratureRules.h b/Sofa/framework/FEM/src/sofa/fem/QuadratureRules.h index 15f10e443f0..60417b88311 100644 --- a/Sofa/framework/FEM/src/sofa/fem/QuadratureRules.h +++ b/Sofa/framework/FEM/src/sofa/fem/QuadratureRules.h @@ -486,4 +486,58 @@ struct QuadratureDispatch } }; +/** + * 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 +struct FiniteElementQuadrature +{ + using Domain = ReferenceDomain_t; + using Dispatch = QuadratureDispatch; + using QuadraturePointAndWeight = QuadraturePointAndWeight_t; + + 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 + static constexpr auto quadraturePoints() + { + static_assert(Degree >= MinimumQuadratureDegree, + "quadrature degree too low for the interpolation order of this element"); + return Dispatch::template points(); + } + + /// Quadrature points and weights of a degree known only at runtime; view of the + /// compile-time table. + static std::span 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(degree); + } +}; + } From 2f3c194e493ddda032090062d2ab861d20b4549e Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Fri, 11 Sep 2026 16:20:32 +0200 Subject: [PATCH 12/18] [FEM] Bridge reference domain for quadratic/linear edges --- .../sofa/fem/FiniteElement[QuadraticEdge].h | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h index 4606613f486..ac31df9a355 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h @@ -35,9 +35,9 @@ struct FiniteElement FINITEELEMENT_HEADER(sofa::geometry::QuadraticEdge, DataTypes, 1); constexpr static std::array 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& getElementSequence(sofa::core::topology::BaseMeshTopology& topology) @@ -48,20 +48,22 @@ struct FiniteElement static constexpr sofa::type::Vec shapeFunctions(const sofa::type::Vec& q) { const Real xi = q[0]; + constexpr Real half = static_cast(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(1) - xi * xi // mid-edge node, at xi = 0 }; } static constexpr sofa::type::Mat gradientShapeFunctions(const sofa::type::Vec& q) { const Real xi = q[0]; + constexpr Real half = static_cast(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 }; } From 8b83e4e41d15baf0ca1a09d29546742dd771f7b5 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Fri, 11 Sep 2026 16:20:52 +0200 Subject: [PATCH 13/18] [FEM] Adopt Quadrature rule dispatch --- .../FEM/src/sofa/fem/FiniteElement.h | 1 + .../FEM/src/sofa/fem/FiniteElement[Edge].h | 43 +----------- .../src/sofa/fem/FiniteElement[Hexahedron].h | 67 +------------------ .../FEM/src/sofa/fem/FiniteElement[Prism].h | 38 +---------- .../FEM/src/sofa/fem/FiniteElement[Pyramid].h | 41 +----------- .../FEM/src/sofa/fem/FiniteElement[Quad].h | 61 +---------------- .../sofa/fem/FiniteElement[QuadraticEdge].h | 19 +----- .../fem/FiniteElement[QuadraticHexahedron].h | 28 +------- .../sofa/fem/FiniteElement[QuadraticQuad].h | 21 +----- .../fem/FiniteElement[QuadraticTetrahedron].h | 24 +------ .../fem/FiniteElement[QuadraticTriangle].h | 19 +----- .../src/sofa/fem/FiniteElement[Tetrahedron].h | 44 +----------- .../src/sofa/fem/FiniteElement[Triangle].h | 40 +---------- 13 files changed, 13 insertions(+), 433 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h index a490fddd6a3..c3f67148b0a 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 +#include #include namespace sofa::fem diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h index 53570c30d05..7d6fb368878 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h @@ -21,8 +21,6 @@ ******************************************************************************/ #pragma once #include -#include -#include #if !defined(SOFA_FEM_FINITE_ELEMENT_EDGE_CPP) #include @@ -33,6 +31,7 @@ namespace sofa::fem template struct FiniteElement + : FiniteElementQuadrature, 1> { FINITEELEMENT_HEADER(sofa::geometry::Edge, DataTypes, 1); @@ -56,46 +55,6 @@ struct FiniteElement SOFA_UNUSED(q); return {{-static_cast(0.5)}, {static_cast(0.5)}}; } - - template - static constexpr auto quadraturePoints() - { - 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(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::quadratureRule: unsupported degree"); - } - } - }; #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 d855c9f8c72..129c978c748 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h @@ -21,8 +21,6 @@ ******************************************************************************/ #pragma once #include -#include -#include #if !defined(SOFA_FEM_FINITE_ELEMENT_HEXAHEDRON_CPP) #include @@ -33,6 +31,7 @@ namespace sofa::fem template struct FiniteElement + : FiniteElementQuadrature, 3> { FINITEELEMENT_HEADER(sofa::geometry::Hexahedron, DataTypes, 3); static_assert(spatial_dimensions == 3, "Hexahedrons are only defined in 3D"); @@ -97,70 +96,6 @@ struct FiniteElement return gradient; } - - template - static constexpr auto quadraturePoints() - { - 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); - - 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"); - } - } - - // Quadrature rule selector by degree; view of the compile-time table. - static std::span 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::quadratureRule: unsupported degree"); - } - } }; #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 8951dc385f3..ab1f264bc76 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h @@ -21,8 +21,6 @@ ******************************************************************************/ #pragma once #include -#include -#include #if !defined(SOFA_FEM_FINITE_ELEMENT_PRISM_CPP) #include @@ -33,6 +31,7 @@ namespace sofa::fem template struct FiniteElement + : FiniteElementQuadrature, 1> { FINITEELEMENT_HEADER(sofa::geometry::Prism, DataTypes, 3); static_assert(spatial_dimensions == 3, "Prisms are only defined in 3D"); @@ -75,41 +74,6 @@ struct FiniteElement {0, q[2], q[1]}, }; } - - template - static constexpr auto quadraturePoints() - { - if constexpr (Degree <= 1) - { - // Degree 1: 2-point rule (default). - constexpr auto third = static_cast(1) / static_cast(3); - constexpr auto sqrt_3 = static_cast(0.57735026919); // 1/sqrt(3) - constexpr auto one = static_cast(1); - constexpr QuadraturePoint q0 {third, third, static_cast(0.5) * (one - sqrt_3)}; - constexpr QuadraturePoint q1 {third, third, static_cast(0.5) * (one + sqrt_3)}; - - constexpr std::array q { - std::make_pair(q0, 1./4.), - std::make_pair(q1, 1./4.), - }; - return q; - } - else - { - static_assert(Degree <= 1, "FiniteElement: no quadrature rule for the requested degree"); - } - } - - // Quadrature rule selector by degree; view of the compile-time table. - static std::span quadratureRule(sofa::Size degree) - { - switch (degree) - { - case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; } - default: - throw std::invalid_argument("FiniteElement::quadratureRule: unsupported degree"); - } - } }; #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 422fcb920f3..be118ca4171 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h @@ -22,8 +22,6 @@ #pragma once #include #include -#include -#include #if !defined(SOFA_FEM_FINITE_ELEMENT_PYRAMID_CPP) #include @@ -34,6 +32,7 @@ namespace sofa::fem template struct FiniteElement + : FiniteElementQuadrature, 1> { FINITEELEMENT_HEADER(sofa::geometry::Pyramid, DataTypes, 3); static_assert(spatial_dimensions == 3, "Pyramids are only defined in 3D"); @@ -72,44 +71,6 @@ struct FiniteElement { 0, 0, static_cast(0.5)} }; } - - template - static constexpr auto quadraturePoints() - { - if constexpr (Degree <= 1) - { - // Degree 1: 8-point rule (default). - constexpr Real sqrt3_1 = static_cast(1) / static_cast(1.73205080757); - constexpr Real one = static_cast(1); - - constexpr std::array q { - 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}, - }; - return q; - } - else - { - static_assert(Degree <= 1, "FiniteElement: no quadrature rule for the requested degree"); - } - } - - // Quadrature rule selector by degree; view of the compile-time table. - static std::span quadratureRule(sofa::Size degree) - { - switch (degree) - { - case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; } - default: - throw std::invalid_argument("FiniteElement::quadratureRule: unsupported degree"); - } - } }; #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 261eeae7036..0373a938559 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h @@ -21,8 +21,6 @@ ******************************************************************************/ #pragma once #include -#include -#include #if !defined(SOFA_FEM_FINITE_ELEMENT_QUAD_CPP) #include @@ -33,6 +31,7 @@ namespace sofa::fem template struct FiniteElement + : FiniteElementQuadrature, 2> { FINITEELEMENT_HEADER(sofa::geometry::Quad, DataTypes, 2); static_assert(spatial_dimensions > 1, "Quads cannot be defined in 1D"); @@ -68,64 +67,6 @@ struct FiniteElement {1 / static_cast(4) * (-static_cast(1) - q[1]), 1 / static_cast(4) * ( static_cast(1) - q[0])} }; } - - template - static constexpr auto quadraturePoints() - { - 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 ReferenceCoord q0(sqrt2_3, 0.); - constexpr ReferenceCoord q1(-1/sqrt6, -1./sqrt2); - constexpr ReferenceCoord q2(-1/sqrt6, 1./sqrt2); - - 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(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::quadratureRule: unsupported degree"); - } - } }; #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 ac31df9a355..23362a79857 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h @@ -31,6 +31,7 @@ namespace sofa::fem template struct FiniteElement + : FiniteElementQuadrature, 3> { FINITEELEMENT_HEADER(sofa::geometry::QuadraticEdge, DataTypes, 1); @@ -66,24 +67,6 @@ struct FiniteElement {-2 * xi} // mid-edge node }; } - - static constexpr std::array 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 q0(a); - constexpr sofa::type::Vec q1(b); - - constexpr std::array 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..8418c9e5eb4 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].h @@ -31,6 +31,7 @@ namespace sofa::fem template struct FiniteElement + : FiniteElementQuadrature, 5> { FINITEELEMENT_HEADER(sofa::geometry::QuadraticHexahedron, DataTypes, 3); static_assert(spatial_dimensions == 3, "Quadratic Hexahedrons are only defined in 3D"); @@ -219,33 +220,6 @@ struct FiniteElement 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 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 pts = {-a, 0.0, a}; - constexpr std::array wts = {w_a, w_0, w_a}; - - std::array 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(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..b70e32c97ef 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h @@ -31,6 +31,7 @@ namespace sofa::fem template struct FiniteElement + : FiniteElementQuadrature, 3> { FINITEELEMENT_HEADER(sofa::geometry::QuadraticQuad, DataTypes, 2); @@ -95,26 +96,6 @@ struct FiniteElement {-2 * xi * (1 - eta * eta), -2 * eta * (1 - xi * xi)} }; } - - static constexpr std::array quadraturePoints() - { - // constexpr Real a = 1. / std::sqrt(3.); - constexpr Real a { 0.57735026919 }; - constexpr Real w = 1.0; - - constexpr sofa::type::Vec q0(-a, -a); - constexpr sofa::type::Vec q1(a, -a); - constexpr sofa::type::Vec q2(a, a); - constexpr sofa::type::Vec q3(-a, a); - - constexpr std::array 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..c5cc0e2926a 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].h @@ -31,6 +31,7 @@ namespace sofa::fem template struct FiniteElement + : FiniteElementQuadrature, 2> { FINITEELEMENT_HEADER(sofa::geometry::QuadraticTetrahedron, DataTypes, 3); static_assert(spatial_dimensions == 3, "Quadratic Tetrahedrons are only defined in 3D"); @@ -104,29 +105,6 @@ struct FiniteElement {0, 4 * l3, 4 * l2} }; } - - static constexpr std::array 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 q0(a, b, b); - constexpr sofa::type::Vec q1(b, a, b); - constexpr sofa::type::Vec q2(b, b, a); - constexpr sofa::type::Vec q3(b, b, b); - - constexpr std::array 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..34e938ad97a 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].h @@ -31,6 +31,7 @@ namespace sofa::fem template struct FiniteElement + : FiniteElementQuadrature, 2> { FINITEELEMENT_HEADER(sofa::geometry::QuadraticTriangle, DataTypes, 2); @@ -85,24 +86,6 @@ struct FiniteElement {-4 * l2, 4 * (l0 - l2)} }; } - - static constexpr std::array quadraturePoints() - { - constexpr Real a = 2. / 3.; - constexpr Real b = 1. / 6.; - constexpr Real w = 1. / 6.; - - constexpr sofa::type::Vec q0(b, b); - constexpr sofa::type::Vec q1(a, b); - constexpr sofa::type::Vec q2(b, a); - - constexpr std::array 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 bca61671e0a..43db9b5c21c 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h @@ -21,8 +21,6 @@ ******************************************************************************/ #pragma once #include -#include -#include #if !defined(SOFA_FEM_FINITE_ELEMENT_TETAHEDRON_CPP) #include @@ -33,6 +31,7 @@ namespace sofa::fem template struct FiniteElement + : FiniteElementQuadrature, 1> { FINITEELEMENT_HEADER(sofa::geometry::Tetrahedron, DataTypes, 3); static_assert(spatial_dimensions == 3, "Tetrahedrons are only defined in 3D"); @@ -69,47 +68,6 @@ struct FiniteElement {0, 0, 1} }; } - - template - static constexpr auto quadraturePoints() - { - 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(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::quadratureRule: unsupported degree"); - } - } }; #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 84c9ec19a74..bfb4ba200d8 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h @@ -21,8 +21,6 @@ ******************************************************************************/ #pragma once #include -#include -#include namespace sofa::fem { @@ -33,6 +31,7 @@ namespace sofa::fem template struct FiniteElement + : FiniteElementQuadrature, 1> { FINITEELEMENT_HEADER(sofa::geometry::Triangle, DataTypes, 2); static_assert(spatial_dimensions > 1, "Triangles cannot be defined in 1D"); @@ -65,43 +64,6 @@ struct FiniteElement {0, 1} }; } - - template - static constexpr auto quadraturePoints() - { - 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(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::quadratureRule: unsupported degree"); - } - } }; #if !defined(SOFA_FEM_FINITE_ELEMENT_TRIANGLE_CPP) From d461e8178c1c8b1ad844b11295477766e59be46f Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Fri, 11 Sep 2026 16:35:37 +0200 Subject: [PATCH 14/18] Remove is-a and go for has-a --- Sofa/framework/FEM/src/sofa/fem/FiniteElement.h | 13 ++++++++++++- .../FEM/src/sofa/fem/FiniteElement[Edge].h | 3 +-- .../FEM/src/sofa/fem/FiniteElement[Hexahedron].h | 3 +-- .../FEM/src/sofa/fem/FiniteElement[Prism].h | 3 +-- .../FEM/src/sofa/fem/FiniteElement[Pyramid].h | 3 +-- .../FEM/src/sofa/fem/FiniteElement[Quad].h | 3 +-- .../FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h | 3 +-- .../sofa/fem/FiniteElement[QuadraticHexahedron].h | 3 +-- .../FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h | 3 +-- .../sofa/fem/FiniteElement[QuadraticTetrahedron].h | 3 +-- .../src/sofa/fem/FiniteElement[QuadraticTriangle].h | 3 +-- .../FEM/src/sofa/fem/FiniteElement[Tetrahedron].h | 3 +-- .../FEM/src/sofa/fem/FiniteElement[Triangle].h | 3 +-- Sofa/framework/FEM/test/FiniteElement_test.cpp | 10 ++++++---- 14 files changed, 30 insertions(+), 29 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h index c3f67148b0a..6e73274ff68 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h @@ -30,7 +30,7 @@ namespace sofa::fem template struct FiniteElement; -#define FINITEELEMENT_HEADER(ElType, DataTypes, dimension) \ +#define FINITEELEMENT_HEADER(ElType, DataTypes, dimension, defaultQuadratureDegree) \ using Coord = sofa::Coord_t;\ using Real = sofa::Real_t;\ using ElementType = ElType;\ @@ -42,6 +42,17 @@ struct FiniteElement; using ShapeFunctionType = std::function;\ using QuadraturePoint = ReferenceCoord; \ using QuadraturePointAndWeight = std::pair;\ + /* Quadrature rules of this element's reference domain. The two accessors are */\ + /* re-exported so callers keep writing FiniteElement::quadraturePoints(). */\ + using Quadrature = FiniteElementQuadrature;\ + static constexpr sofa::Size InterpolationOrder = Quadrature::InterpolationOrder;\ + static constexpr sofa::Size MinimumQuadratureDegree = Quadrature::MinimumQuadratureDegree;\ + static constexpr sofa::Size DefaultQuadratureDegree = Quadrature::DefaultQuadratureDegree;\ + template \ + static constexpr auto quadraturePoints()\ + { return Quadrature::template quadraturePoints(); }\ + static auto quadratureRule(sofa::Size degree)\ + { return Quadrature::quadratureRule(degree); }\ using Helper = FiniteElementHelper diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h index 7d6fb368878..2afaaf46be0 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h @@ -31,9 +31,8 @@ namespace sofa::fem template struct FiniteElement - : FiniteElementQuadrature, 1> { - FINITEELEMENT_HEADER(sofa::geometry::Edge, DataTypes, 1); + FINITEELEMENT_HEADER(sofa::geometry::Edge, DataTypes, 1, 1); constexpr static std::array referenceElementNodes {{ReferenceCoord{-1}, ReferenceCoord{1}}}; diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h index 129c978c748..bbb7a52111d 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h @@ -31,9 +31,8 @@ namespace sofa::fem template struct FiniteElement - : FiniteElementQuadrature, 3> { - 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: diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h index ab1f264bc76..f05ad349363 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h @@ -31,9 +31,8 @@ namespace sofa::fem template struct FiniteElement - : FiniteElementQuadrature, 1> { - 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 referenceElementNodes {{ diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h index be118ca4171..b83318b7a92 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h @@ -32,9 +32,8 @@ namespace sofa::fem template struct FiniteElement - : FiniteElementQuadrature, 1> { - 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 referenceElementNodes {{ diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h index 0373a938559..a7580dad7c7 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h @@ -31,9 +31,8 @@ namespace sofa::fem template struct FiniteElement - : FiniteElementQuadrature, 2> { - 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 referenceElementNodes {{ diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h index 23362a79857..f94d68691ac 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h @@ -31,9 +31,8 @@ namespace sofa::fem template struct FiniteElement - : FiniteElementQuadrature, 3> { - FINITEELEMENT_HEADER(sofa::geometry::QuadraticEdge, DataTypes, 1); + FINITEELEMENT_HEADER(sofa::geometry::QuadraticEdge, DataTypes, 1, 3); constexpr static std::array referenceElementNodes {{ ReferenceCoord{-1}, // vertex 0 diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].h index 8418c9e5eb4..875b02c4a26 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].h @@ -31,9 +31,8 @@ namespace sofa::fem template struct FiniteElement - : FiniteElementQuadrature, 5> { - 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 referenceElementNodes {{ diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h index b70e32c97ef..0cbd2a43f3a 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h @@ -31,9 +31,8 @@ namespace sofa::fem template struct FiniteElement - : FiniteElementQuadrature, 3> { - FINITEELEMENT_HEADER(sofa::geometry::QuadraticQuad, DataTypes, 2); + FINITEELEMENT_HEADER(sofa::geometry::QuadraticQuad, DataTypes, 2, 3); constexpr static std::array referenceElementNodes {{ {-1, -1}, // vertex 0 diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].h index c5cc0e2926a..f9f45b28a9a 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].h @@ -31,9 +31,8 @@ namespace sofa::fem template struct FiniteElement - : FiniteElementQuadrature, 2> { - 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 referenceElementNodes {{ diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].h index 34e938ad97a..3063a0679bb 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTriangle].h @@ -31,9 +31,8 @@ namespace sofa::fem template struct FiniteElement - : FiniteElementQuadrature, 2> { - FINITEELEMENT_HEADER(sofa::geometry::QuadraticTriangle, DataTypes, 2); + FINITEELEMENT_HEADER(sofa::geometry::QuadraticTriangle, DataTypes, 2, 2); constexpr static std::array referenceElementNodes {{ {0, 0}, // vertex 0 diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h index 43db9b5c21c..ec074d17e89 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h @@ -31,9 +31,8 @@ namespace sofa::fem template struct FiniteElement - : FiniteElementQuadrature, 1> { - 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 referenceElementNodes {{ diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h index bfb4ba200d8..83465f57ae0 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h @@ -31,9 +31,8 @@ namespace sofa::fem template struct FiniteElement - : FiniteElementQuadrature, 1> { - 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 referenceElementNodes {{ diff --git a/Sofa/framework/FEM/test/FiniteElement_test.cpp b/Sofa/framework/FEM/test/FiniteElement_test.cpp index c607d06fe15..a4ad8a5c44f 100644 --- a/Sofa/framework/FEM/test/FiniteElement_test.cpp +++ b/Sofa/framework/FEM/test/FiniteElement_test.cpp @@ -26,15 +26,17 @@ namespace sofa { /** - * For each quadrature degree in [1, MaxDegree], 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 void testSumWeights(const sofa::Real_t expected) { using FE = sofa::fem::FiniteElement; + static_assert(MaxDegree >= FE::MinimumQuadratureDegree); - for (sofa::Size degree = 1; degree <= MaxDegree; ++degree) + for (sofa::Size degree = FE::MinimumQuadratureDegree; degree <= MaxDegree; ++degree) { SReal weightSum = 0; for (const auto& [q, w] : FE::quadratureRule(degree)) @@ -98,7 +100,7 @@ TEST(FiniteElement, pyramid3dWeights) TEST(FiniteElement, quadraticHexa3dWeights) { - testSumWeights(8); + testSumWeights(8); } /** From e8d49f6f6ca10bbebba252568db6b0ec88b37339 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 14 Sep 2026 10:16:32 +0200 Subject: [PATCH 15/18] remove comments --- Sofa/framework/FEM/src/sofa/fem/FiniteElement.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h b/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h index 6e73274ff68..7723e86aec5 100644 --- a/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h +++ b/Sofa/framework/FEM/src/sofa/fem/FiniteElement.h @@ -42,8 +42,6 @@ struct FiniteElement; using ShapeFunctionType = std::function;\ using QuadraturePoint = ReferenceCoord; \ using QuadraturePointAndWeight = std::pair;\ - /* Quadrature rules of this element's reference domain. The two accessors are */\ - /* re-exported so callers keep writing FiniteElement::quadraturePoints(). */\ using Quadrature = FiniteElementQuadrature;\ static constexpr sofa::Size InterpolationOrder = Quadrature::InterpolationOrder;\ static constexpr sofa::Size MinimumQuadratureDegree = Quadrature::MinimumQuadratureDegree;\ From 53f50bf0c2eb5a3dbe3be6f6cc43a26033c886d1 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 1 Oct 2026 15:42:59 +0200 Subject: [PATCH 16/18] Use normalized difference for evaluating error in test --- .../component/mass/testing/MassTestCreation.h | 24 +++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) 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..1bfd9e36f96 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,9 @@ #include #include +#include +#include + namespace sofa::component::mass::testing { @@ -131,7 +134,10 @@ struct Mass_test : public sofa::testing::BaseSimulationTest, public sofa::testin sofa::type::vector expectedMv; sofa::testing::data_traits::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 +173,10 @@ struct Mass_test : public sofa::testing::BaseSimulationTest, public sofa::testin sofa::type::vector expectedMv; sofa::testing::data_traits::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 +190,11 @@ struct Mass_test : public sofa::testing::BaseSimulationTest, public sofa::testin vMv += static_cast(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 +248,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); From b159486fb1e5c4b5c2c9f562844816876c76e244 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 1 Oct 2026 15:59:39 +0200 Subject: [PATCH 17/18] Stricten errorMax as is affordable --- Sofa/Component/Mass/tests/MassTestCreation[FEMMass].cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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_testm_testAccFromF = false; this->m_testAddMToMatrix = false; - this->m_errorMax = 1e3_sreal; + this->m_errorMax = 1e2_sreal; auto topology = sofa::core::objectmodel::New(); this->m_node->addObject(topology); From f15663fe4645a4ee37b624c86a69c047c95e2b4b Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Thu, 1 Oct 2026 16:04:54 +0200 Subject: [PATCH 18/18] Remove useless include --- .../Testing/src/sofa/component/mass/testing/MassTestCreation.h | 1 - 1 file changed, 1 deletion(-) 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 1bfd9e36f96..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 @@ -33,7 +33,6 @@ #include #include -#include namespace sofa::component::mass::testing {