Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-17 01:02:13 +00:00
Code Issues Packages Projects Releases Wiki Activity

refactor(quadrature): Rename gltet.jl → gl_tetrahedra.jl with modern API

Browse Source 
- 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)
This commit is contained in:
Jukka Aho
2025-11-21 00:39:27 +02:00
parent 8c443a98fb
commit 4cf1cda237
2 changed files with 160 additions and 67 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/quadrature/gl_tetrahedra.jl
+160
View File
@@ -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
src/quadrature/gltet.jl
-67
View File
@@ -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