refactor(quadrature): Rename glpyr.jl → gl_pyramids.jl with modern API

- Deleted: src/quadrature/glpyr.jl (45 lines, old zip-based API)
- Created: src/quadrature/gl_pyramids.jl (86 lines, modern type-based API)

Key improvements:
- Two pyramid rule variants: GaussLegendre{2}() (default) and GaussLegendre{2,:B}()
- QuadraturePoint{3} with Vec{3} coordinates (not tuples)
- SVector return type (zero allocation)
- Comprehensive documentation explaining pyramid singularity at apex
- Inline functions for performance
- Legacy Val{:GLPYR5} compatibility maintained

Technical details:
- Default variant: 4 base points + 1 elevated, non-uniform weights
- Variant B: 5 points with uniform weights (a = 2/15)
- Both exact for degree 1 polynomials
- Pyramid: quad base at z=0, apex at (0,0,1)

Net: +41 lines (added extensive documentation and modern types)
This commit is contained in:
Jukka Aho
2025-11-21 00:37:54 +02:00
parent 9f242df9db
commit 0b2195c7c2
2 changed files with 86 additions and 45 deletions
+86
View File
@@ -0,0 +1,86 @@
# 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 pyramid elements.
This file implements quadrature rules for 3D pyramid reference elements
with a quadrilateral base and triangular apex.
# Available Rules
- `GaussLegendre{2}()`: 5-point rule (default variant), exact for linear
- `GaussLegendre{2,:B}()`: 5-point alternative variant, exact for linear
# Notes
- Pyramid has quadrilateral base at z=0 and apex at (0,0,1)
- Quadrature for pyramids is non-trivial due to singularity at apex
- Both variants have order 1 accuracy
- Higher order rules for pyramids are complex and rarely implemented
See also: [`GaussLegendre`](@ref), [`QuadraturePoint`](@ref)
"""
# ============================================================================
# GaussLegendre{2}: 5-point rules (two variants)
# ============================================================================
"""
get_quadrature_points(::Type{Pyramid}, ::GaussLegendre{2})
5-point Gauss-Legendre rule for pyramid (default variant).
Exact for linear polynomials (degree 1).
4 points on base plane + 1 point elevated toward apex.
"""
@inline function get_quadrature_points(::Type{Pyramid}, ::GaussLegendre{2,:default})
g1 = 0.5842373946721771876874344
g2 = -2/3
g3 = 2/5
w1 = 81/100
w2 = 125/27
return SVector(
QuadraturePoint(Vec{3}(-g1, -g1, g2), w1),
QuadraturePoint(Vec{3}( g1, -g1, g2), w1),
QuadraturePoint(Vec{3}( g1, g1, g2), w1),
QuadraturePoint(Vec{3}(-g1, g1, g2), w1),
QuadraturePoint(Vec{3}(0.0, 0.0, g3), w2)
)
end
"""
get_quadrature_points(::Type{Pyramid}, ::GaussLegendre{2,:B})
5-point Gauss-Legendre rule for pyramid (variant B).
Exact for linear polynomials (degree 1).
Alternative point distribution with uniform weights.
"""
@inline function get_quadrature_points(::Type{Pyramid}, ::GaussLegendre{2,:B})
a = 2/15
h1 = 0.1531754163448146
h2 = 0.6372983346207416
return SVector(
QuadraturePoint(Vec{3}( 0.5, 0.0, h1), a),
QuadraturePoint(Vec{3}( 0.0, 0.5, h1), a),
QuadraturePoint(Vec{3}(-0.5, 0.0, h1), a),
QuadraturePoint(Vec{3}( 0.0, -0.5, h1), a),
QuadraturePoint(Vec{3}( 0.0, 0.0, h2), a)
)
end
# ============================================================================
# Legacy symbol-based API (deprecated, kept for backwards compatibility)
# ============================================================================
# These map old :GLPYR symbols to new API
get_quadrature_points(::Type{Val{:GLPYR5}}) =
get_quadrature_points(Pyramid, GaussLegendre{2}())
get_quadrature_points(::Type{Val{:GLPYR5B}}) =
get_quadrature_points(Pyramid, GaussLegendre{2,:B}())
# Order queries (deprecated)
get_order(::Type{Val{:GLPYR5}}) = 1
get_order(::Type{Val{:GLPYR5B}}) = 1
-45
View File
@@ -1,45 +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 pyramid elements
""" Gauss-Legendre quadrature, 5 point rule on pyramid. """
function get_quadrature_points(::Type{Val{:GLPYR5}})
g1 = 0.5842373946721771876874344
g2 = -2.0/3.0
g3 = 2.0/5.0
w1 = 81.0/100.0
w2 = 125.0/27.0
weights = (w1, w1, w1, w1, w2)
points = (
(-g1, -g1, g2),
( g1, -g1, g2),
( g1, g1, g2),
(-g1, g1, g2),
(0.0, 0.0, g3))
return zip(weights, points)
end
function get_order(::Type{Val{:GLPYR5}})
return 1
end
""" Gauss-Legendre quadrature, 5 point rule on pyramid. """
function get_quadrature_points(::Type{Val{:GLPYR5B}})
a = 2.0/15.0
h1 = 0.1531754163448146
h2 = 0.6372983346207416
weights = (a, a, a, a, a)
points = (
(0.5, 0.0, h1),
(0.0, 0.5, h1),
(-0.5, 0.0, h1),
(0.0, -0.5, h1),
(0.0, 0.0, h2)
)
return zip(weights, points)
end
function get_order(::Type{Val{:GLPYR5B}})
return 1
end