mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-19 01:48:47 +00:00
refactor(quadrature): Rename glwed.jl → gl_wedges.jl with modern API
- 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)
This commit is contained in:
Jukka Aho
@@ -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
|
||||
@@ -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
|
||||
Reference in New Issue
Block a user