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 6 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
6 commits
Select commit Hold shift + click to select a range
1dc7ec0
[FEM] Introduce higher order quadrature rule switch for Triangle
th-skam Aug 10, 2026
9ec5242
Edge
th-skam Aug 10, 2026
213ccaf
Tetrahedra
th-skam Aug 10, 2026
cfb35e2
Quad
th-skam Aug 10, 2026
5ad33f1
Hexa
th-skam Aug 10, 2026
4b6a9dc
Use uint for degree
th-skam Aug 11, 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
43 changes: 38 additions & 5 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 @@ -21,6 +21,8 @@
******************************************************************************/
#pragma once
#include <sofa/fem/FiniteElement.h>
#include <span>
#include <stdexcept>

#if !defined(SOFA_FEM_FINITE_ELEMENT_EDGE_CPP)
#include <sofa/defaulttype/VecTypes.h>
Expand Down Expand Up @@ -55,12 +57,43 @@ struct FiniteElement<sofa::geometry::Edge, DataTypes>
return {{-static_cast<Real>(0.5)}, {static_cast<Real>(0.5)}};
}

static constexpr std::array<QuadraturePointAndWeight, 1> quadraturePoints()
template <sofa::Size Degree = 1>
static constexpr auto quadraturePoints()
{
constexpr sofa::type::Vec<TopologicalDimension, Real> q0(static_cast<Real>(0));
return {
std::make_pair(q0, static_cast<Real>(2))
};
if constexpr (Degree <= 1)
{
// Degree 1: 1-point midpoint rule (default).
return std::array<QuadraturePointAndWeight, 1>{
std::make_pair(ReferenceCoord(static_cast<Real>(0)), static_cast<Real>(2))
};
}
else if constexpr (Degree <= 3)
{
// Degrees 2-3: 2-point Gauss-Legendre rule.
constexpr Real sqrt3 = 1.73205080757;
constexpr Real g = static_cast<Real>(1) / sqrt3;
return std::array<QuadraturePointAndWeight, 2>{
std::make_pair(ReferenceCoord(-g), static_cast<Real>(1)),
std::make_pair(ReferenceCoord( g), static_cast<Real>(1))
};
}
else
{
static_assert(Degree <= 3, "FiniteElement<Edge>: no quadrature rule for the requested degree");
}
}

// Quadrature rule selector by degree; view of the compile-time table.
static std::span<const QuadraturePointAndWeight> quadratureRule(sofa::Size degree)
{
switch (degree)
{
case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; }
case 2:
case 3: { static constexpr auto rule = quadraturePoints<3>(); return rule; }
default:
throw std::invalid_argument("FiniteElement<Edge>::quadratureRule: unsupported degree");
}
}

};
Expand Down
74 changes: 60 additions & 14 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 @@ -21,6 +21,8 @@
******************************************************************************/
#pragma once
#include <sofa/fem/FiniteElement.h>
#include <span>
#include <stdexcept>

#if !defined(SOFA_FEM_FINITE_ELEMENT_HEXAHEDRON_CPP)
#include <sofa/defaulttype/VecTypes.h>
Expand Down Expand Up @@ -96,24 +98,68 @@ struct FiniteElement<sofa::geometry::Hexahedron, DataTypes>
return gradient;
}

template <sofa::Size Degree = 3>
static constexpr auto quadraturePoints()
{
constexpr Real sqrt3 = 1.73205080757; //sqrt(3.)
constexpr Real sqrt3_1 = static_cast<Real>(1) / sqrt3;
constexpr Real one = static_cast<Real>(1);
if constexpr (Degree <= 1)
{
// Degree 1: 1-point centroid rule.
return std::array<QuadraturePointAndWeight, 1>{
std::make_pair(ReferenceCoord(static_cast<Real>(0), static_cast<Real>(0), static_cast<Real>(0)), static_cast<Real>(8))
};
}
else if constexpr (Degree <= 3)
{
// Degrees 2-3: 2x2x2 Gauss-Legendre rule (default).
constexpr Real sqrt3 = 1.73205080757; //sqrt(3.)
constexpr Real sqrt3_1 = static_cast<Real>(1) / sqrt3;
constexpr Real one = static_cast<Real>(1);

constexpr std::array q {
std::pair{referenceElementNodes[0] * sqrt3_1, one},
std::pair{referenceElementNodes[1] * sqrt3_1, one},
std::pair{referenceElementNodes[2] * sqrt3_1, one},
std::pair{referenceElementNodes[3] * sqrt3_1, one},
std::pair{referenceElementNodes[4] * sqrt3_1, one},
std::pair{referenceElementNodes[5] * sqrt3_1, one},
std::pair{referenceElementNodes[6] * sqrt3_1, one},
std::pair{referenceElementNodes[7] * sqrt3_1, one},
};
return std::array {
std::pair{referenceElementNodes[0] * sqrt3_1, one},
std::pair{referenceElementNodes[1] * sqrt3_1, one},
std::pair{referenceElementNodes[2] * sqrt3_1, one},
std::pair{referenceElementNodes[3] * sqrt3_1, one},
std::pair{referenceElementNodes[4] * sqrt3_1, one},
std::pair{referenceElementNodes[5] * sqrt3_1, one},
std::pair{referenceElementNodes[6] * sqrt3_1, one},
std::pair{referenceElementNodes[7] * sqrt3_1, one},
};
}
else if constexpr (Degree <= 5)
{
// Degrees 4-5: 3x3x3 Gauss-Legendre rule.
constexpr Real g = 0.77459666924; //sqrt(3./5.)
constexpr std::array<Real, 3> node{ -g, static_cast<Real>(0), g };
constexpr std::array<Real, 3> weight{ static_cast<Real>(5./9.), static_cast<Real>(8./9.), static_cast<Real>(5./9.) };

std::array<QuadraturePointAndWeight, 27> q{};
sofa::Size k = 0;
for (sofa::Size i = 0; i < 3; ++i)
for (sofa::Size j = 0; j < 3; ++j)
for (sofa::Size l = 0; l < 3; ++l)
q[k++] = std::make_pair(ReferenceCoord(node[i], node[j], node[l]), weight[i] * weight[j] * weight[l]);
return q;
}
else
{
static_assert(Degree <= 5, "FiniteElement<Hexahedron>: no quadrature rule for the requested degree");
}
}

return q;
// Quadrature rule selector by degree; view of the compile-time table.
static std::span<const QuadraturePointAndWeight> quadratureRule(sofa::Size degree)
{
switch (degree)
{
case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; }
case 2:
case 3: { static constexpr auto rule = quadraturePoints<3>(); return rule; }
case 4:
case 5: { static constexpr auto rule = quadraturePoints<5>(); return rule; }
default:
throw std::invalid_argument("FiniteElement<Hexahedron>::quadratureRule: unsupported degree");
}
}
};

Expand Down
67 changes: 55 additions & 12 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 @@ -21,6 +21,8 @@
******************************************************************************/
#pragma once
#include <sofa/fem/FiniteElement.h>
#include <span>
#include <stdexcept>

#if !defined(SOFA_FEM_FINITE_ELEMENT_QUAD_CPP)
#include <sofa/defaulttype/VecTypes.h>
Expand Down Expand Up @@ -67,21 +69,62 @@ struct FiniteElement<sofa::geometry::Quad, DataTypes>
};
}

static constexpr std::array<QuadraturePointAndWeight, 3> quadraturePoints()
template <sofa::Size Degree = 2>
static constexpr auto quadraturePoints()
{
constexpr Real sqrt2_3 = 0.816496580928; //sqrt(2./3.)
constexpr Real sqrt6 = 2.44948974278; //sqrt(6.)
constexpr Real sqrt2 = 1.41421356237; //sqrt(2.)
if constexpr (Degree <= 1)
{
// Degree 1: 1-point centroid rule.
return std::array<QuadraturePointAndWeight, 1>{
std::make_pair(ReferenceCoord(static_cast<Real>(0), static_cast<Real>(0)), static_cast<Real>(4))
};
}
else if constexpr (Degree <= 2)
{
// Degree 2: 3-point rule (default).
constexpr Real sqrt2_3 = 0.816496580928; //sqrt(2./3.)
constexpr Real sqrt6 = 2.44948974278; //sqrt(6.)
constexpr Real sqrt2 = 1.41421356237; //sqrt(2.)

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

return {
std::make_pair(q0, 4./3.),
std::make_pair(q1, 4./3.),
std::make_pair(q2, 4./3.),
};
return std::array<QuadraturePointAndWeight, 3>{
std::make_pair(q0, 4./3.),
std::make_pair(q1, 4./3.),
std::make_pair(q2, 4./3.)
};
}
else if constexpr (Degree <= 3)
{
// Degree 3: 2x2 Gauss-Legendre rule.
constexpr Real sqrt3 = 1.73205080757; //sqrt(3.)
constexpr Real g = static_cast<Real>(1) / sqrt3;
return std::array<QuadraturePointAndWeight, 4>{
std::make_pair(ReferenceCoord(-g, -g), static_cast<Real>(1)),
std::make_pair(ReferenceCoord( g, -g), static_cast<Real>(1)),
std::make_pair(ReferenceCoord( g, g), static_cast<Real>(1)),
std::make_pair(ReferenceCoord(-g, g), static_cast<Real>(1))
};
}
else
{
static_assert(Degree <= 3, "FiniteElement<Quad>: no quadrature rule for the requested degree");
}
}

// Quadrature rule selector by degree; view of the compile-time table.
static std::span<const QuadraturePointAndWeight> quadratureRule(sofa::Size degree)
{
switch (degree)
{
case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; }
case 2: { static constexpr auto rule = quadraturePoints<2>(); return rule; }
case 3: { static constexpr auto rule = quadraturePoints<3>(); return rule; }
default:
throw std::invalid_argument("FiniteElement<Quad>::quadratureRule: unsupported degree");
}
}
};

Expand Down
44 changes: 40 additions & 4 deletions Sofa/framework/FEM/src/sofa/fem/FiniteElement[Tetrahedron].h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,8 @@
******************************************************************************/
#pragma once
#include <sofa/fem/FiniteElement.h>
#include <span>
#include <stdexcept>

#if !defined(SOFA_FEM_FINITE_ELEMENT_TETAHEDRON_CPP)
#include <sofa/defaulttype/VecTypes.h>
Expand Down Expand Up @@ -68,11 +70,45 @@ struct FiniteElement<sofa::geometry::Tetrahedron, DataTypes>
};
}

static constexpr std::array<QuadraturePointAndWeight, 1> quadraturePoints()
template <sofa::Size Degree = 1>
static constexpr auto quadraturePoints()
{
constexpr sofa::type::Vec<TopologicalDimension, Real> q0(1./4., 1./4., 1./4.);
constexpr std::array<QuadraturePointAndWeight, 1> q { std::make_pair(q0, 1./6.) };
return q;
if constexpr (Degree <= 1)
{
// Degree 1: 1-point centroid rule (default).
return std::array<QuadraturePointAndWeight, 1>{
std::make_pair(ReferenceCoord(1./4., 1./4., 1./4.), Real(1./6.))
};
}
else if constexpr (Degree <= 2)
{
// Degree 2: 4-point rule.
constexpr Real sqrt5 = 2.2360679774997896;
constexpr Real a = (5. - sqrt5) / 20.;
constexpr Real b = (5. + 3. * sqrt5) / 20.;
return std::array<QuadraturePointAndWeight, 4>{
std::make_pair(ReferenceCoord(a, a, a), Real(1./24.)),
std::make_pair(ReferenceCoord(b, a, a), Real(1./24.)),
std::make_pair(ReferenceCoord(a, b, a), Real(1./24.)),
std::make_pair(ReferenceCoord(a, a, b), Real(1./24.))
};
}
else
{
static_assert(Degree <= 2, "FiniteElement<Tetrahedron>: no quadrature rule for the requested degree");
}
}

// Quadrature rule selector by degree; view of the compile-time table.
static std::span<const QuadraturePointAndWeight> quadratureRule(sofa::Size degree)
{
switch (degree)
{
case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; }
case 2: { static constexpr auto rule = quadraturePoints<2>(); return rule; }
default:
throw std::invalid_argument("FiniteElement<Tetrahedron>::quadratureRule: unsupported degree");
}
}
};

Expand Down
40 changes: 36 additions & 4 deletions Sofa/framework/FEM/src/sofa/fem/FiniteElement[Triangle].h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,8 @@
******************************************************************************/
#pragma once
#include <sofa/fem/FiniteElement.h>
#include <span>
#include <stdexcept>

namespace sofa::fem
{
Expand Down Expand Up @@ -64,11 +66,41 @@ struct FiniteElement<sofa::geometry::Triangle, DataTypes>
};
}

static constexpr std::array<QuadraturePointAndWeight, 1> quadraturePoints()
template <sofa::Size Degree = 1>
static constexpr auto quadraturePoints()
{
return {
std::make_pair(sofa::type::Vec<TopologicalDimension, Real>(1./3., 1./3.), 1./2.)
};
if constexpr (Degree <= 1)
{
// Degree 1: 1-point centroid rule (default).
return std::array<QuadraturePointAndWeight, 1>{
std::make_pair(ReferenceCoord(1./3., 1./3.), Real(1./2.))
};
}
else if constexpr (Degree <= 2)
{
// Degree 2: 3-point interior rule.
return std::array<QuadraturePointAndWeight, 3>{
std::make_pair(ReferenceCoord(1./6., 1./6.), Real(1./6.)),
std::make_pair(ReferenceCoord(2./3., 1./6.), Real(1./6.)),
std::make_pair(ReferenceCoord(1./6., 2./3.), Real(1./6.))
};
}
else
{
static_assert(Degree <= 2, "FiniteElement<Triangle>: no quadrature rule for the requested degree");
}
}

// Quadrature rule selector by degree; view of the compile-time table.
static std::span<const QuadraturePointAndWeight> quadratureRule(sofa::Size degree)
{
switch (degree)
{
case 1: { static constexpr auto rule = quadraturePoints<1>(); return rule; }
case 2: { static constexpr auto rule = quadraturePoints<2>(); return rule; }
default:
throw std::invalid_argument("FiniteElement<Triangle>::quadratureRule: unsupported degree");
}
}
};

Expand Down
Loading