From 4cf1cda2375cdf47255ab8045b562f11f7c48dfe Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Fri, 21 Nov 2025 00:39:27 +0200 Subject: [PATCH] =?UTF-8?q?refactor(quadrature):=20Rename=20gltet.jl=20?= =?UTF-8?q?=E2=86=92=20gl=5Ftetrahedra.jl=20with=20modern=20API?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Deleted: src/quadrature/gltet.jl (67 lines, old zip-based API) - Created: src/quadrature/gl_tetrahedra.jl (160 lines, modern type-based API) Key improvements: - QuadraturePoint{3} with Vec{3} coordinates (not tuples) - SVector return types (zero allocation) - @inline directives for performance - Comprehensive documentation with warnings - Four rules from 1-point to 15-point (degrees 1-4) Implemented rules: - GaussLegendre{1}(): 1 point (centroid, degree 1) - GaussLegendre{2}(): 4 points (degree 2, symmetric placement) - GaussLegendre{4}(): 15 points (degree 4, all weights positive) Technical details: - Reference tetrahedron: vertices at (0,0,0), (1,0,0), (0,1,0), (0,0,1) - All weights sum to 1/6 (volume of reference tetrahedron) - Parametric domain constraint: ξ + η + ζ ≤ 1 - Legacy Val{:GLTET*} compatibility maintained Warning documented: - GaussLegendre{3}() has negative weight at centroid (-2/15) - May cause numerical issues, GaussLegendre{2}() or {4}() recommended Net: +93 lines (added extensive documentation, modern types, and safety warnings) --- src/quadrature/gl_tetrahedra.jl | 160 ++++++++++++++++++++++++++++++++ src/quadrature/gltet.jl | 67 ------------- 2 files changed, 160 insertions(+), 67 deletions(-) create mode 100644 src/quadrature/gl_tetrahedra.jl delete mode 100644 src/quadrature/gltet.jl diff --git a/src/quadrature/gl_tetrahedra.jl b/src/quadrature/gl_tetrahedra.jl new file mode 100644 index 0000000..970dc54 --- /dev/null +++ b/src/quadrature/gl_tetrahedra.jl @@ -0,0 +1,160 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +""" +Gauss-Legendre quadrature rules for tetrahedral elements. + +This file implements quadrature rules for 3D tetrahedral reference elements +in the parametric domain with vertices at (0,0,0), (1,0,0), (0,1,0), (0,0,1). + +# Available Rules +- `GaussLegendre{1}()`: 1-point rule (centroid), exact for linear +- `GaussLegendre{2}()`: 4-point rule, exact for quadratic +- `GaussLegendre{3}()`: 5-point rule, exact for cubic (has negative weight!) +- `GaussLegendre{4}()`: 15-point rule, exact for quartic + +# Notes +- Reference tetrahedron has constraint ξ + η + ζ ≤ 1 +- All weights sum to 1/6 (volume of reference tetrahedron) +- GLTET5 has one negative weight, may cause numerical issues +- Higher order rules (>4) are available in literature but rarely used + +See also: [`GaussLegendre`](@ref), [`QuadraturePoint`](@ref) +""" + +# ============================================================================ +# GaussLegendre{1}: 1-point rule (centroid) +# ============================================================================ + +""" + get_quadrature_points(::Type{Tetrahedron}, ::GaussLegendre{1}) + +1-point Gauss-Legendre rule for tetrahedron (centroid). +Exact for linear polynomials (degree 1). +""" +@inline function get_quadrature_points(::Type{Tetrahedron}, ::GaussLegendre{1,V}) where V + return SVector( + QuadraturePoint(Vec{3}(1 / 4, 1 / 4, 1 / 4), 1 / 6) + ) +end + +# ============================================================================ +# GaussLegendre{2}: 4-point rule +# ============================================================================ + +""" + get_quadrature_points(::Type{Tetrahedron}, ::GaussLegendre{2}) + +4-point Gauss-Legendre rule for tetrahedron. +Exact for quadratic polynomials (degree 2). + +Points are symmetrically placed inside the tetrahedron. +""" +@inline function get_quadrature_points(::Type{Tetrahedron}, ::GaussLegendre{2,V}) where V + a = (5 + 3 * sqrt(5)) / 20 + b = (5 - sqrt(5)) / 20 + w = 1 / 24 + + return SVector( + QuadraturePoint(Vec{3}(a, b, b), w), + QuadraturePoint(Vec{3}(b, a, b), w), + QuadraturePoint(Vec{3}(b, b, a), w), + QuadraturePoint(Vec{3}(b, b, b), w) + ) +end + +# ============================================================================ +# GaussLegendre{3}: 5-point rule +# ============================================================================ + +""" + get_quadrature_points(::Type{Tetrahedron}, ::GaussLegendre{3}) + +5-point Gauss-Legendre rule for tetrahedron. +Exact for cubic polynomials (degree 3). + +**WARNING:** This rule has one negative weight (-2/15), which may cause +numerical issues. Consider using GaussLegendre{2}() or {4}() instead. +""" +@inline function get_quadrature_points(::Type{Tetrahedron}, ::GaussLegendre{3,V}) where V + a = 1 / 4 + b = 1 / 6 + c = 1 / 2 + + return SVector( + QuadraturePoint(Vec{3}(a, a, a), -2 / 15), # Negative weight! + QuadraturePoint(Vec{3}(b, b, b), 3 / 40), + QuadraturePoint(Vec{3}(b, b, c), 3 / 40), + QuadraturePoint(Vec{3}(b, c, b), 3 / 40), + QuadraturePoint(Vec{3}(c, b, b), 3 / 40) + ) +end + +# ============================================================================ +# GaussLegendre{4}: 15-point rule +# ============================================================================ + +""" + get_quadrature_points(::Type{Tetrahedron}, ::GaussLegendre{4}) + +15-point Gauss-Legendre rule for tetrahedron. +Exact for quartic polynomials (degree 4). + +This is a high-accuracy rule suitable for quadratic basis functions. +All weights are positive. +""" +@inline function get_quadrature_points(::Type{Tetrahedron}, ::GaussLegendre{4,V}) where V + a = 1 / 4 + b1 = (7 + sqrt(15)) / 34 + b2 = (7 - sqrt(15)) / 34 + c1 = (13 - 3 * sqrt(15)) / 34 + c2 = (13 + 3 * sqrt(15)) / 34 + d = (5 - sqrt(15)) / 20 + f = (5 + sqrt(15)) / 20 + + w1 = 8 / 405 + w2 = (2665 - 14 * sqrt(15)) / 226800 + w3 = (2665 + 14 * sqrt(15)) / 226800 + w4 = 5 / 567 + + return SVector( + QuadraturePoint(Vec{3}(a, a, a), w1), + QuadraturePoint(Vec{3}(b1, b1, b1), w2), + QuadraturePoint(Vec{3}(b1, b1, c1), w2), + QuadraturePoint(Vec{3}(b1, c1, b1), w2), + QuadraturePoint(Vec{3}(c1, b1, b1), w2), + QuadraturePoint(Vec{3}(b2, b2, b2), w3), + QuadraturePoint(Vec{3}(b2, b2, c2), w3), + QuadraturePoint(Vec{3}(b2, c2, b2), w3), + QuadraturePoint(Vec{3}(c2, b2, b2), w3), + QuadraturePoint(Vec{3}(d, d, f), w4), + QuadraturePoint(Vec{3}(d, f, d), w4), + QuadraturePoint(Vec{3}(f, d, d), w4), + QuadraturePoint(Vec{3}(d, f, f), w4), + QuadraturePoint(Vec{3}(f, d, f), w4), + QuadraturePoint(Vec{3}(f, f, d), w4) + ) +end + +# ============================================================================ +# Legacy symbol-based API (deprecated, kept for backwards compatibility) +# ============================================================================ + +# These map old :GLTET symbols to new API +get_quadrature_points(::Type{Val{:GLTET1}}) = + get_quadrature_points(Tetrahedron, GaussLegendre{1}()) + +get_quadrature_points(::Type{Val{:GLTET4}}) = + get_quadrature_points(Tetrahedron, GaussLegendre{2}()) + +get_quadrature_points(::Type{Val{:GLTET5}}) = + get_quadrature_points(Tetrahedron, GaussLegendre{3}()) + +get_quadrature_points(::Type{Val{:GLTET15}}) = + get_quadrature_points(Tetrahedron, GaussLegendre{4}()) + +# Order queries (deprecated) +get_order(::Type{Val{:GLTET1}}) = 1 +get_order(::Type{Val{:GLTET4}}) = 2 +get_order(::Type{Val{:GLTET5}}) = 3 +get_order(::Type{Val{:GLTET15}}) = 4 diff --git a/src/quadrature/gltet.jl b/src/quadrature/gltet.jl deleted file mode 100644 index 96b7c8e..0000000 --- a/src/quadrature/gltet.jl +++ /dev/null @@ -1,67 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE - -### Gauss quadrature rules for tetrahedrons - -""" Gauss-Legendre quadrature, 1 point rule on tetrahedron. """ -function get_quadrature_points(::Type{Val{:GLTET1}}) - weights = (1.0/6.0, ) - points = ((1.0/4.0, 1.0/4.0, 1.0/4.0), ) - return zip(weights, points) -end - -function get_order(::Type{Val{:GLTET1}}) - return 1 -end - -""" Gauss-Legendre quadrature, 4 point rule on tetrahedron. """ -function get_quadrature_points(::Type{Val{:GLTET4}}) - a = (5.0+3.0*sqrt(5.0))/20.0 - b = (5.0-sqrt(5.0))/20.0 - w = 1.0/24.0 - weights = (w, w, w, w) - points = ((a, b, b), (b, a, b), (b, b, a), (b, b, b)) - return zip(weights, points) -end - -function get_order(::Type{Val{:GLTET4}}) - return 2 -end - -""" Gauss-Legendre quadrature, 5 point rule on tetrahedron. """ -function get_quadrature_points(::Type{Val{:GLTET5}}) - a = 1.0/4.0 - b = 1.0/6.0 - c = 1.0/2.0 - weights = (-2.0/15.0, 3.0/40.0, 3.0/40.0, 3.0/40.0, 3.0/40.0) - points = ((a, a, a), (b, b, b), (b, b, c), (b, c, b), (c, b, b)) - return zip(weights, points) -end - -function get_order(::Type{Val{:GLTET5}}) - return 3 -end - -""" Gauss-Legendre quadrature, 15 point rule on tetrahedron. """ -function get_quadrature_points(::Type{Val{:GLTET15}}) - a = 1.0/4.0 - b1 = 1.0/34.0*(7.0 + sqrt(15.0)) - b2 = 1.0/34.0*(7.0 - sqrt(15.0)) - c1 = 1.0/34.0*(13.0 - 3.0*sqrt(15.0)) - c2 = 1.0/34.0*(13.0 + 3.0*sqrt(15.0)) - d = 1.0/20.0*(5.0 - sqrt(15.0)) - f = 1.0/20.0*(5.0 + sqrt(15.0)) - w1 = 8.0/405.0 - w2 = (2665.0 - 14.0*sqrt(15.0))/226800.0 - w3 = (2665.0 + 14.0*sqrt(15.0))/226800.0 - w4 = 5.0/567.0 - weights = (w1, w2, w2, w2, w2, w3, w3, w3, w3, w4, w4, w4, w4, w4, w4) - points = ((a, a, a), (b1, b1, b1), (b1, b1, c1), (b1, c1, b1), (c1, b1, b1), - (b2, b2, b2), (b2, b2, c2), (b2, c2, b2), (c2, b2, b2), (d, d, f), - (d, f, d), (f, d, d), (d, f, f), (f, d, f), (f, f, d)) - return zip(weights, points) -end - -function get_order(::Type{Val{:GLTET15}}) - return 4 -end