From 0b2195c7c2fd337f6a9fc7d0e16cd70d6bc49e9e Mon Sep 17 00:00:00 2001
From: Jukka Aho <ahojukka5@gmail.com>
Date: Fri, 21 Nov 2025 00:37:54 +0200
Subject: [PATCH] =?UTF-8?q?refactor(quadrature):=20Rename=20glpyr.jl=20?=
 =?UTF-8?q?=E2=86=92=20gl=5Fpyramids.jl=20with=20modern=20API?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

- 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)
---
 src/quadrature/gl_pyramids.jl | 86 +++++++++++++++++++++++++++++++++++
 src/quadrature/glpyr.jl       | 45 ------------------
 2 files changed, 86 insertions(+), 45 deletions(-)
 create mode 100644 src/quadrature/gl_pyramids.jl
 delete mode 100644 src/quadrature/glpyr.jl

diff --git a/src/quadrature/gl_pyramids.jl b/src/quadrature/gl_pyramids.jl
new file mode 100644
index 0000000..300dd87
--- /dev/null
+++ b/src/quadrature/gl_pyramids.jl
@@ -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
diff --git a/src/quadrature/glpyr.jl b/src/quadrature/glpyr.jl
deleted file mode 100644
index f49da64..0000000
--- a/src/quadrature/glpyr.jl
+++ /dev/null
@@ -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