From 6ad8c43e7e2f295cb5d75aaca29e87fb99ce8324 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Fri, 21 Nov 2025 00:39:54 +0200 Subject: [PATCH] =?UTF-8?q?refactor(quadrature):=20Rename=20glwed.jl=20?= =?UTF-8?q?=E2=86=92=20gl=5Fwedges.jl=20with=20modern=20API?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Deleted: src/quadrature/glwed.jl (93 lines, old zip-based API) - Created: src/quadrature/gl_wedges.jl (142 lines, modern type-based API) Key improvements: - QuadraturePoint{3} with Vec{3} coordinates (not tuples) - SVector return types (zero allocation) - @inline directives for performance - Two 6-point variants (default + :B) - Comprehensive documentation explaining tensor product structure - Three rules: 6-point (2 variants) and 21-point Implemented rules: - GaussLegendre{2}(): 6 points (degree 1, default variant) → 3 triangle points × 2 segment points - GaussLegendre{2,:B}(): 6 points (degree 1, alternative triangle distribution) - GaussLegendre{5}(): 21 points (degree 5, high accuracy) → 7 triangle points × 3 segment points Technical details: - Wedge is tensor product of triangle (ξ,η) and segment (ζ) - Reference wedge: triangular cross-section extruded in z-direction - 6-point rules: 2D triangle rule × 2-point 1D Gauss - 21-point rule: 7-point triangle × 3-point 1D Gauss - Legacy Val{:GLWED*} compatibility maintained Net: +49 lines (added detailed documentation and modern type infrastructure) --- src/quadrature/gl_wedges.jl | 142 ++++++++++++++++++++++++++++++++++++ src/quadrature/glwed.jl | 93 ----------------------- 2 files changed, 142 insertions(+), 93 deletions(-) create mode 100644 src/quadrature/gl_wedges.jl delete mode 100644 src/quadrature/glwed.jl diff --git a/src/quadrature/gl_wedges.jl b/src/quadrature/gl_wedges.jl new file mode 100644 index 0000000..83cd0db --- /dev/null +++ b/src/quadrature/gl_wedges.jl @@ -0,0 +1,142 @@ +# 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 wedge (prism) elements. + +This file implements quadrature rules for 3D wedge/prism reference elements +formed by extruding a triangle in the z-direction. + +# Available Rules +- `GaussLegendre{2}()`: 6-point rule (default variant), exact for linear +- `GaussLegendre{2,:B}()`: 6-point alternative variant, exact for linear +- `GaussLegendre{5}()`: 21-point rule, exact for quintic + +# Notes +- Wedge is tensor product of triangle (ξ,η) and segment (ζ) +- Rules are typically products of triangle rules × 1D Gauss rules +- 6-point rules: 3 triangle points × 2 segment points +- 21-point rule: 7 triangle points × 3 segment points + +See also: [`GaussLegendre`](@ref), [`QuadraturePoint`](@ref) +""" + +# ============================================================================ +# GaussLegendre{2}: 6-point rules (two variants) +# ============================================================================ + +""" + get_quadrature_points(::Type{Wedge}, ::GaussLegendre{2}) + +6-point Gauss-Legendre rule for wedge (default variant). +Exact for linear polynomials (degree 1). + +Tensor product: 3 triangle points × 2 segment points. +""" +@inline function get_quadrature_points(::Type{Wedge}, ::GaussLegendre{2,:default}) + w = 1/6 + a = sqrt(1/3) + + return SVector( + QuadraturePoint(Vec{3}(0.5, 0.0, -a), w), + QuadraturePoint(Vec{3}(0.0, 0.5, -a), w), + QuadraturePoint(Vec{3}(0.5, 0.5, -a), w), + QuadraturePoint(Vec{3}(0.5, 0.0, a), w), + QuadraturePoint(Vec{3}(0.0, 0.5, a), w), + QuadraturePoint(Vec{3}(0.5, 0.5, a), w) + ) +end + +""" + get_quadrature_points(::Type{Wedge}, ::GaussLegendre{2,:B}) + +6-point Gauss-Legendre rule for wedge (variant B). +Exact for linear polynomials (degree 1). + +Alternative triangle point distribution. +""" +@inline function get_quadrature_points(::Type{Wedge}, ::GaussLegendre{2,:B}) + w = 1/6 + a = sqrt(1/3) + + return SVector( + QuadraturePoint(Vec{3}(2/3, 1/6, -a), w), + QuadraturePoint(Vec{3}(1/6, 2/3, -a), w), + QuadraturePoint(Vec{3}(1/6, 1/6, -a), w), + QuadraturePoint(Vec{3}(2/3, 1/6, a), w), + QuadraturePoint(Vec{3}(1/6, 2/3, a), w), + QuadraturePoint(Vec{3}(1/6, 1/6, a), w) + ) +end + +# ============================================================================ +# GaussLegendre{5}: 21-point rule +# ============================================================================ + +""" + get_quadrature_points(::Type{Wedge}, ::GaussLegendre{5}) + +21-point Gauss-Legendre rule for wedge. +Exact for quintic polynomials (degree 5). + +Tensor product: 7 triangle points × 3 segment points. +High-accuracy rule suitable for quadratic basis functions. +""" +@inline function get_quadrature_points(::Type{Wedge}, ::GaussLegendre{5,V}) where V + # 1D Gauss-Legendre 3-point rule + alpha = sqrt(3/5) + c1 = 5/9 + c2 = 8/9 + + # Triangle 7-point rule parameters + a = (6 + sqrt(15)) / 21 + b = (6 - sqrt(15)) / 21 + + return SVector( + # z = -alpha layer (7 points) + QuadraturePoint(Vec{3}(1/3, 1/3, -alpha), c1 * 9/80), + QuadraturePoint(Vec{3}(a, a, -alpha), c1 * (155 + sqrt(15))/2400), + QuadraturePoint(Vec{3}(1-2a, a, -alpha), c1 * (155 + sqrt(15))/2400), + QuadraturePoint(Vec{3}(a, 1-2a, -alpha), c1 * (155 + sqrt(15))/2400), + QuadraturePoint(Vec{3}(b, b, -alpha), c1 * (155 - sqrt(15))/2400), + QuadraturePoint(Vec{3}(1-2b, b, -alpha), c1 * (155 - sqrt(15))/2400), + QuadraturePoint(Vec{3}(b, 1-2b, -alpha), c1 * (155 - sqrt(15))/2400), + + # z = 0 layer (7 points) + QuadraturePoint(Vec{3}(1/3, 1/3, 0.0), c2 * 9/80), + QuadraturePoint(Vec{3}(a, a, 0.0), c2 * (155 + sqrt(15))/2400), + QuadraturePoint(Vec{3}(1-2a, a, 0.0), c2 * (155 + sqrt(15))/2400), + QuadraturePoint(Vec{3}(a, 1-2a, 0.0), c2 * (155 + sqrt(15))/2400), + QuadraturePoint(Vec{3}(b, b, 0.0), c2 * (155 - sqrt(15))/2400), + QuadraturePoint(Vec{3}(1-2b, b, 0.0), c2 * (155 - sqrt(15))/2400), + QuadraturePoint(Vec{3}(b, 1-2b, 0.0), c2 * (155 - sqrt(15))/2400), + + # z = alpha layer (7 points) + QuadraturePoint(Vec{3}(1/3, 1/3, alpha), c1 * 9/80), + QuadraturePoint(Vec{3}(a, a, alpha), c1 * (155 + sqrt(15))/2400), + QuadraturePoint(Vec{3}(1-2a, a, alpha), c1 * (155 + sqrt(15))/2400), + QuadraturePoint(Vec{3}(a, 1-2a, alpha), c1 * (155 + sqrt(15))/2400), + QuadraturePoint(Vec{3}(b, b, alpha), c1 * (155 - sqrt(15))/2400), + QuadraturePoint(Vec{3}(1-2b, b, alpha), c1 * (155 - sqrt(15))/2400), + QuadraturePoint(Vec{3}(b, 1-2b, alpha), c1 * (155 - sqrt(15))/2400) + ) +end + +# ============================================================================ +# Legacy symbol-based API (deprecated, kept for backwards compatibility) +# ============================================================================ + +# These map old :GLWED symbols to new API +get_quadrature_points(::Type{Val{:GLWED6}}) = + get_quadrature_points(Wedge, GaussLegendre{2}()) + +get_quadrature_points(::Type{Val{:GLWED6B}}) = + get_quadrature_points(Wedge, GaussLegendre{2,:B}()) + +get_quadrature_points(::Type{Val{:GLWED21}}) = + get_quadrature_points(Wedge, GaussLegendre{5}()) + +# Order queries (deprecated) +get_order(::Type{Val{:GLWED6}}) = 1 +get_order(::Type{Val{:GLWED6B}}) = 1 +get_order(::Type{Val{:GLWED21}}) = 5 diff --git a/src/quadrature/glwed.jl b/src/quadrature/glwed.jl deleted file mode 100644 index 2bb4ef0..0000000 --- a/src/quadrature/glwed.jl +++ /dev/null @@ -1,93 +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 prismatic elements (wedge) - -""" Gauss-Legendre quadrature, 6 point rule on wedge. """ -function get_quadrature_points(::Type{Val{:GLWED6}}) - w = 1.0/6.0 - a = sqrt(1.0/3.0) - weights = (w, w, w, w, w, w) - points = ((0.5, 0.0, -a), (0.0, 0.5, -a), (0.5, 0.5, -a), - (0.5, 0.0, a), (0.0, 0.5, a), (0.5, 0.5, a)) - return zip(weights, points) -end - -function get_order(::Type{Val{:GLWED6}}) - return 1 -end - -""" Gauss-Legendre quadrature, 6 point rule on wedge. """ -function get_quadrature_points(::Type{Val{:GLWED6B}}) - w = 1.0/6.0 - a = sqrt(1.0/3.0) - weights = (w, w, w, w, w, w) - points = ((2.0/3.0, 1.0/6.0, -a), (1.0/6.0, 2.0/3.0, -a), (1.0/6.0, 1.0/6.0, -a), - (2.0/3.0, 1.0/6.0, a), (1.0/6.0, 2.0/3.0, a), (1.0/6.0, 1.0/6.0, a)) - return zip(weights, points) -end - -function get_order(::Type{Val{:GLWED6B}}) - return 1 -end - -""" Gauss-Legendre quadrature, 21 point rule on wedge. """ -function get_quadrature_points(::Type{Val{:GLWED21}}) - alpha = sqrt(3/5) - c1 = 5/9 - c2 = 8/9 - a = (6+sqrt(15))/21 - b = (6-sqrt(15))/21 - - weights = ( - c1*9/80, - c1*((155+sqrt(15))/2400), - c1*((155+sqrt(15))/2400), - c1*((155+sqrt(15))/2400), - c1*((155-sqrt(15))/2400), - c1*((155-sqrt(15))/2400), - c1*((155-sqrt(15))/2400), - c2*9/80, - c2*((155+sqrt(15))/2400), - c2*((155+sqrt(15))/2400), - c2*((155+sqrt(15))/2400), - c2*((155-sqrt(15))/2400), - c2*((155-sqrt(15))/2400), - c2*((155-sqrt(15))/2400), - c1*9/80, - c1*((155+sqrt(15))/2400), - c1*((155+sqrt(15))/2400), - c1*((155+sqrt(15))/2400), - c1*((155-sqrt(15))/2400), - c1*((155-sqrt(15))/2400), - c1*((155-sqrt(15))/2400)) - - points = ( - (1/3, 1/3, -alpha), - (a, a, -alpha), - (1-2a, a, -alpha), - (a, 1-2a, -alpha), - (b, b, -alpha), - (1-2b, b, -alpha), - (b, 1-2b, -alpha), - (1/3, 1/3, 0), - (a, a, 0), - (1-2a, a, 0), - (a, 1-2a, 0), - (b, b, 0), - (1-2b, b, 0), - (b, 1-2b, 0), - (1/3, 1/3, alpha), - (a, a, alpha), - (1-2a, a, alpha), - (a, 1-2a, alpha), - (b, b, alpha), - (1-2b, b, alpha), - (b, 1-2b, alpha)) - - return zip(weights, points) -end - -function get_order(::Type{Val{:GLWED21}}) - return 5 -end