Skip to content

[FEM] Introduce higher accuracy quadrature rules #6225

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
th-skam wants to merge 15 commits into
base: master
Choose a base branch
from
Open

[FEM] Introduce higher accuracy quadrature rules #6225

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
f70bd97
[FEM] Introduce higher order quadrature rule switch for Triangle
th-skam Aug 10, 2026
542d147
Edge
th-skam Aug 10, 2026
2f86678
Tetrahedra
th-skam Aug 10, 2026
e0a5955
Quad
th-skam Aug 10, 2026
39d2942
Hexa
th-skam Aug 10, 2026
4abd0a2
Use uint for degree
th-skam Aug 11, 2026
85b9f79
Adopt quadratureRule in prism and pyramid
th-skam Sep 4, 2026
f6b83fc
Include quadratureRule in FE test
th-skam Sep 4, 2026
b6569cd
Add polynomial order to linear elements
th-skam Sep 11, 2026
fe026d9
[FEM] Add QuadratureRules.h
th-skam Sep 11, 2026
1681cc5
Add the ElementQuadrature class
th-skam Sep 11, 2026
2f3c194
[FEM] Bridge reference domain for quadratic/linear edges
th-skam Sep 11, 2026
8b83e4e
[FEM] Adopt Quadrature rule dispatch
th-skam Sep 11, 2026
d461e81
Remove is-a and go for has-a
th-skam Sep 11, 2026
e8d49f6
remove comments
th-skam Sep 14, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions Sofa/framework/FEM/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
12 changes: 11 additions & 1 deletion Sofa/framework/FEM/src/sofa/fem/FiniteElement.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
******************************************************************************/
#pragma once
#include <sofa/fem/config.h>
#include <sofa/fem/QuadratureRules.h>
#include <sofa/core/topology/BaseMeshTopology.h>

namespace sofa::fem
Expand All @@ -29,7 +30,7 @@ namespace sofa::fem
template <class ElementType, class DataTypes>
struct FiniteElement;

#define FINITEELEMENT_HEADER(ElType, DataTypes, dimension) \
#define FINITEELEMENT_HEADER(ElType, DataTypes, dimension, defaultQuadratureDegree) \
using Coord = sofa::Coord_t<DataTypes>;\
using Real = sofa::Real_t<DataTypes>;\
using ElementType = ElType;\
Expand All @@ -41,6 +42,15 @@ struct FiniteElement;
using ShapeFunctionType = std::function<Real(const ReferenceCoord&)>;\
using QuadraturePoint = ReferenceCoord; \
using QuadraturePointAndWeight = std::pair<QuadraturePoint, Real>;\
using Quadrature = FiniteElementQuadrature<ElType, Real, defaultQuadratureDegree>;\
static constexpr sofa::Size InterpolationOrder = Quadrature::InterpolationOrder;\
static constexpr sofa::Size MinimumQuadratureDegree = Quadrature::MinimumQuadratureDegree;\
static constexpr sofa::Size DefaultQuadratureDegree = Quadrature::DefaultQuadratureDegree;\
template <sofa::Size Degree = Quadrature::DefaultQuadratureDegree>\
static constexpr auto quadraturePoints()\
{ return Quadrature::template quadraturePoints<Degree>(); }\
static auto quadratureRule(sofa::Size degree)\
{ return Quadrature::quadratureRule(degree); }\
using Helper = FiniteElementHelper<ElementType, DataTypes>


Expand Down
11 changes: 1 addition & 10 deletions Sofa/framework/FEM/src/sofa/fem/FiniteElement[Edge].h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ namespace sofa::fem
template <class DataTypes>
struct FiniteElement<sofa::geometry::Edge, DataTypes>
{
FINITEELEMENT_HEADER(sofa::geometry::Edge, DataTypes, 1);
FINITEELEMENT_HEADER(sofa::geometry::Edge, DataTypes, 1, 1);

constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{ReferenceCoord{-1}, ReferenceCoord{1}}};

Expand All @@ -54,15 +54,6 @@ struct FiniteElement<sofa::geometry::Edge, DataTypes>
SOFA_UNUSED(q);
return {{-static_cast<Real>(0.5)}, {static_cast<Real>(0.5)}};
}

static constexpr std::array<QuadraturePointAndWeight, 1> quadraturePoints()
{
constexpr sofa::type::Vec<TopologicalDimension, Real> q0(static_cast<Real>(0));
return {
std::make_pair(q0, static_cast<Real>(2))
};
}

};

#if !defined(SOFA_FEM_FINITE_ELEMENT_EDGE_CPP)
Expand Down
22 changes: 1 addition & 21 deletions Sofa/framework/FEM/src/sofa/fem/FiniteElement[Hexahedron].h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ namespace sofa::fem
template <class DataTypes>
struct FiniteElement<sofa::geometry::Hexahedron, DataTypes>
{
FINITEELEMENT_HEADER(sofa::geometry::Hexahedron, DataTypes, 3);
FINITEELEMENT_HEADER(sofa::geometry::Hexahedron, DataTypes, 3, 3);
static_assert(spatial_dimensions == 3, "Hexahedrons are only defined in 3D");

// Following the convention in sofa::geometry::Hexahedron:
Expand Down Expand Up @@ -95,26 +95,6 @@ struct FiniteElement<sofa::geometry::Hexahedron, DataTypes>

return gradient;
}

static constexpr auto quadraturePoints()
{
constexpr Real sqrt3 = 1.73205080757; //sqrt(3.)
constexpr Real sqrt3_1 = static_cast<Real>(1) / sqrt3;
constexpr Real one = static_cast<Real>(1);

constexpr std::array q {
std::pair{referenceElementNodes[0] * sqrt3_1, one},
std::pair{referenceElementNodes[1] * sqrt3_1, one},
std::pair{referenceElementNodes[2] * sqrt3_1, one},
std::pair{referenceElementNodes[3] * sqrt3_1, one},
std::pair{referenceElementNodes[4] * sqrt3_1, one},
std::pair{referenceElementNodes[5] * sqrt3_1, one},
std::pair{referenceElementNodes[6] * sqrt3_1, one},
std::pair{referenceElementNodes[7] * sqrt3_1, one},
};

return q;
}
};

#if !defined(SOFA_FEM_FINITE_ELEMENT_HEXAHEDRON_CPP)
Expand Down
17 changes: 1 addition & 16 deletions Sofa/framework/FEM/src/sofa/fem/FiniteElement[Prism].h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ namespace sofa::fem
template <class DataTypes>
struct FiniteElement<sofa::geometry::Prism, DataTypes>
{
FINITEELEMENT_HEADER(sofa::geometry::Prism, DataTypes, 3);
FINITEELEMENT_HEADER(sofa::geometry::Prism, DataTypes, 3, 1);
static_assert(spatial_dimensions == 3, "Prisms are only defined in 3D");

constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
Expand Down Expand Up @@ -73,21 +73,6 @@ struct FiniteElement<sofa::geometry::Prism, DataTypes>
{0, q[2], q[1]},
};
}

static constexpr std::array<QuadraturePointAndWeight, 2> quadraturePoints()
{
constexpr auto third = static_cast<Real>(1) / static_cast<Real>(3);
constexpr auto sqrt_3 = static_cast<Real>(0.57735026919); // 1/sqrt(3)
constexpr auto one = static_cast<Real>(1);
constexpr QuadraturePoint q0 {third, third, static_cast<Real>(0.5) * (one - sqrt_3)};
constexpr QuadraturePoint q1 {third, third, static_cast<Real>(0.5) * (one + sqrt_3)};

constexpr std::array<QuadraturePointAndWeight, 2> q {
std::make_pair(q0, 1./4.),
std::make_pair(q1, 1./4.),
};
return q;
}
};

#if !defined(SOFA_FEM_FINITE_ELEMENT_PRISM_CPP)
Expand Down
20 changes: 1 addition & 19 deletions Sofa/framework/FEM/src/sofa/fem/FiniteElement[Pyramid].h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ namespace sofa::fem
template <class DataTypes>
struct FiniteElement<sofa::geometry::Pyramid, DataTypes>
{
FINITEELEMENT_HEADER(sofa::geometry::Pyramid, DataTypes, 3);
FINITEELEMENT_HEADER(sofa::geometry::Pyramid, DataTypes, 3, 1);
static_assert(spatial_dimensions == 3, "Pyramids are only defined in 3D");

constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
Expand Down Expand Up @@ -70,24 +70,6 @@ struct FiniteElement<sofa::geometry::Pyramid, DataTypes>
{ 0, 0, static_cast<Real>(0.5)}
};
}

static constexpr auto quadraturePoints()
{
constexpr Real sqrt3_1 = static_cast<Real>(1) / static_cast<Real>(1.73205080757);
constexpr Real one = static_cast<Real>(1);

// We use the 8 Gauss points of the hexahedron, which exactly integrate the (1-z)^2 Jacobian of the pyramid.
return std::array {
std::pair{ReferenceCoord{-sqrt3_1, -sqrt3_1, -sqrt3_1}, one},
std::pair{ReferenceCoord{ sqrt3_1, -sqrt3_1, -sqrt3_1}, one},
std::pair{ReferenceCoord{ sqrt3_1, sqrt3_1, -sqrt3_1}, one},
std::pair{ReferenceCoord{-sqrt3_1, sqrt3_1, -sqrt3_1}, one},
std::pair{ReferenceCoord{-sqrt3_1, -sqrt3_1, sqrt3_1}, one},
std::pair{ReferenceCoord{ sqrt3_1, -sqrt3_1, sqrt3_1}, one},
std::pair{ReferenceCoord{ sqrt3_1, sqrt3_1, sqrt3_1}, one},
std::pair{ReferenceCoord{-sqrt3_1, sqrt3_1, sqrt3_1}, one},
};
}
};

#if !defined(SOFA_FEM_FINITE_ELEMENT_PYRAMID_CPP)
Expand Down
19 changes: 1 addition & 18 deletions Sofa/framework/FEM/src/sofa/fem/FiniteElement[Quad].h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ namespace sofa::fem
template <class DataTypes>
struct FiniteElement<sofa::geometry::Quad, DataTypes>
{
FINITEELEMENT_HEADER(sofa::geometry::Quad, DataTypes, 2);
FINITEELEMENT_HEADER(sofa::geometry::Quad, DataTypes, 2, 2);
static_assert(spatial_dimensions > 1, "Quads cannot be defined in 1D");

constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
Expand Down Expand Up @@ -66,23 +66,6 @@ struct FiniteElement<sofa::geometry::Quad, DataTypes>
{1 / static_cast<Real>(4) * (-static_cast<Real>(1) - q[1]), 1 / static_cast<Real>(4) * ( static_cast<Real>(1) - q[0])}
};
}

static constexpr std::array<QuadraturePointAndWeight, 3> quadraturePoints()
{
constexpr Real sqrt2_3 = 0.816496580928; //sqrt(2./3.)
constexpr Real sqrt6 = 2.44948974278; //sqrt(6.)
constexpr Real sqrt2 = 1.41421356237; //sqrt(2.)

constexpr sofa::type::Vec<TopologicalDimension, Real> q0(sqrt2_3, 0.);
constexpr sofa::type::Vec<TopologicalDimension, Real> q1(-1/sqrt6, -1./sqrt2);
constexpr sofa::type::Vec<TopologicalDimension, Real> q2(-1/sqrt6, 1./sqrt2);

return {
std::make_pair(q0, 4./3.),
std::make_pair(q1, 4./3.),
std::make_pair(q2, 4./3.),
};
}
};

#if !defined(SOFA_FEM_FINITE_ELEMENT_QUAD_CPP)
Expand Down
38 changes: 11 additions & 27 deletions Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticEdge].h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,12 +32,12 @@ namespace sofa::fem
template <class DataTypes>
struct FiniteElement<sofa::geometry::QuadraticEdge, DataTypes>
{
FINITEELEMENT_HEADER(sofa::geometry::QuadraticEdge, DataTypes, 1);
FINITEELEMENT_HEADER(sofa::geometry::QuadraticEdge, DataTypes, 1, 3);

constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
ReferenceCoord{0}, // vertex 0
ReferenceCoord{-1}, // vertex 0
ReferenceCoord{1}, // vertex 1
ReferenceCoord{0.5} // mid-edge node
ReferenceCoord{0} // mid-edge node
}};

static const sofa::type::vector<TopologyElement>& getElementSequence(sofa::core::topology::BaseMeshTopology& topology)
Expand All @@ -48,40 +48,24 @@ struct FiniteElement<sofa::geometry::QuadraticEdge, DataTypes>
static constexpr sofa::type::Vec<NumberOfNodesInElement, Real> shapeFunctions(const sofa::type::Vec<TopologicalDimension, Real>& q)
{
const Real xi = q[0];
constexpr Real half = static_cast<Real>(0.5);
return {
2 * xi * (xi - 0.5) + 1 - 2 * xi, // vertex 0: (2*xi - 1) * (xi - 1) = 2*xi^2 - 3*xi + 1
2 * xi * (xi - 0.5), // vertex 1: (2*xi - 1) * xi
4 * xi * (1 - xi) // mid-edge: 4*xi*(1-xi)
half * xi * (xi - 1), // vertex 0, at xi = -1
half * xi * (xi + 1), // vertex 1, at xi = 1
static_cast<Real>(1) - xi * xi // mid-edge node, at xi = 0
};
}

static constexpr sofa::type::Mat<NumberOfNodesInElement, TopologicalDimension, Real> gradientShapeFunctions(const sofa::type::Vec<TopologicalDimension, Real>& q)
{
const Real xi = q[0];
constexpr Real half = static_cast<Real>(0.5);
return {
{4 * xi - 3}, // vertex 0
{4 * xi - 1}, // vertex 1
{4 - 8 * xi} // mid-edge
{xi - half}, // vertex 0
{xi + half}, // vertex 1
{-2 * xi} // mid-edge node
};
}

static constexpr std::array<QuadraturePointAndWeight, 2> quadraturePoints()
{
// constexpr Real a = (1. - 1. / std::sqrt(3.)) / 2.;
constexpr Real a { 0.211324865405 };
// constexpr Real b = (1. + 1. / std::sqrt(3.)) / 2.;
constexpr Real b { 0.788675134595 };
constexpr Real w = 0.5;

constexpr sofa::type::Vec<TopologicalDimension, Real> q0(a);
constexpr sofa::type::Vec<TopologicalDimension, Real> q1(b);

constexpr std::array<QuadraturePointAndWeight, 2> q {
std::make_pair(q0, w),
std::make_pair(q1, w)
};
return q;
}
};

#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_EDGE_CPP)
Expand Down
29 changes: 1 addition & 28 deletions Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticHexahedron].h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ namespace sofa::fem
template <class DataTypes>
struct FiniteElement<sofa::geometry::QuadraticHexahedron, DataTypes>
{
FINITEELEMENT_HEADER(sofa::geometry::QuadraticHexahedron, DataTypes, 3);
FINITEELEMENT_HEADER(sofa::geometry::QuadraticHexahedron, DataTypes, 3, 5);
static_assert(spatial_dimensions == 3, "Quadratic Hexahedrons are only defined in 3D");

constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
Expand Down Expand Up @@ -219,33 +219,6 @@ struct FiniteElement<sofa::geometry::QuadraticHexahedron, DataTypes>
N_0, N_1, N_2, N_3, N_4, N_5, N_6, N_7, N_8, N_9, N_10, N_11, N_12, N_13, N_14, N_15,
N_16, N_17, N_18, N_19, N_20, N_21, N_22, N_23, N_24, N_25, N_26);
}

static constexpr std::array<QuadraturePointAndWeight, 27> quadraturePoints()
{
constexpr Real a = 0.7745966692414834; // sqrt(3/5)
constexpr Real w_a = 5.0 / 9.0;
constexpr Real w_0 = 8.0 / 9.0;

constexpr std::array<Real, 3> pts = {-a, 0.0, a};
constexpr std::array<Real, 3> wts = {w_a, w_0, w_a};

std::array<QuadraturePointAndWeight, 27> q {};
int index = 0;
for (int i = 0; i < 3; ++i)
{
for (int j = 0; j < 3; ++j)
{
for (int k = 0; k < 3; ++k)
{
q[index++] = std::make_pair(
sofa::type::Vec<TopologicalDimension, Real>(pts[i], pts[j], pts[k]),
wts[i] * wts[j] * wts[k]
);
}
}
}
return q;
}
};

#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_HEXAHEDRON_CPP)
Expand Down
22 changes: 1 addition & 21 deletions Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticQuad].h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ namespace sofa::fem
template <class DataTypes>
struct FiniteElement<sofa::geometry::QuadraticQuad, DataTypes>
{
FINITEELEMENT_HEADER(sofa::geometry::QuadraticQuad, DataTypes, 2);
FINITEELEMENT_HEADER(sofa::geometry::QuadraticQuad, DataTypes, 2, 3);

constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
{-1, -1}, // vertex 0
Expand Down Expand Up @@ -95,26 +95,6 @@ struct FiniteElement<sofa::geometry::QuadraticQuad, DataTypes>
{-2 * xi * (1 - eta * eta), -2 * eta * (1 - xi * xi)}
};
}

static constexpr std::array<QuadraturePointAndWeight, 4> quadraturePoints()
{
// constexpr Real a = 1. / std::sqrt(3.);
constexpr Real a { 0.57735026919 };
constexpr Real w = 1.0;

constexpr sofa::type::Vec<TopologicalDimension, Real> q0(-a, -a);
constexpr sofa::type::Vec<TopologicalDimension, Real> q1(a, -a);
constexpr sofa::type::Vec<TopologicalDimension, Real> q2(a, a);
constexpr sofa::type::Vec<TopologicalDimension, Real> q3(-a, a);

constexpr std::array<QuadraturePointAndWeight, 4> q {
std::make_pair(q0, w),
std::make_pair(q1, w),
std::make_pair(q2, w),
std::make_pair(q3, w)
};
return q;
}
};

#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_QUAD_CPP)
Expand Down
25 changes: 1 addition & 24 deletions Sofa/framework/FEM/src/sofa/fem/FiniteElement[QuadraticTetrahedron].h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ namespace sofa::fem
template <class DataTypes>
struct FiniteElement<sofa::geometry::QuadraticTetrahedron, DataTypes>
{
FINITEELEMENT_HEADER(sofa::geometry::QuadraticTetrahedron, DataTypes, 3);
FINITEELEMENT_HEADER(sofa::geometry::QuadraticTetrahedron, DataTypes, 3, 2);
static_assert(spatial_dimensions == 3, "Quadratic Tetrahedrons are only defined in 3D");

constexpr static std::array<ReferenceCoord, NumberOfNodesInElement> referenceElementNodes {{
Expand Down Expand Up @@ -104,29 +104,6 @@ struct FiniteElement<sofa::geometry::QuadraticTetrahedron, DataTypes>
{0, 4 * l3, 4 * l2}
};
}

static constexpr std::array<QuadraturePointAndWeight, 4> quadraturePoints()
{
// 4-point quadrature rule for quadratic tetrahedron
// constexpr Real a = (5. + 3. * std::sqrt(5.)) / 20.;
constexpr Real a { 0.585410196625 };
// constexpr Real b = (5. - std::sqrt(5.)) / 20.;
constexpr Real b { 0.138196601125 };
constexpr Real w = 1. / 24.;

constexpr sofa::type::Vec<TopologicalDimension, Real> q0(a, b, b);
constexpr sofa::type::Vec<TopologicalDimension, Real> q1(b, a, b);
constexpr sofa::type::Vec<TopologicalDimension, Real> q2(b, b, a);
constexpr sofa::type::Vec<TopologicalDimension, Real> q3(b, b, b);

constexpr std::array<QuadraturePointAndWeight, 4> q {
std::make_pair(q0, w),
std::make_pair(q1, w),
std::make_pair(q2, w),
std::make_pair(q3, w)
};
return q;
}
};

#if !defined(SOFA_FEM_FINITE_ELEMENT_QUADRATIC_TETAHEDRON_CPP)
Expand Down
Loading
Loading