mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-19 01:48:47 +00:00
refactor(quadrature): Rename glquad.jl → gl_tensor_product.jl with modern API
- Deleted: src/quadrature/glquad.jl (40 lines, programmatic generation only)
- Created: src/quadrature/gl_tensor_product.jl (344 lines, complete implementations)
Key improvements:
- Explicit implementations for Segment, Quadrilateral, Hexahedron
- All rules GaussLegendre{1} through GaussLegendre{5}
- QuadraturePoint{D} with Vec{D} coordinates (not tuples)
- SVector return types for zero allocation
- @inline directives for performance
- Comprehensive documentation per function
Implemented rules:
Segments (1D):
- GaussLegendre{1}(): 1 point (exact degree 1)
- GaussLegendre{2}(): 2 points (exact degree 3)
- GaussLegendre{3}(): 3 points (exact degree 5)
- GaussLegendre{4}(): 4 points (exact degree 7)
- GaussLegendre{5}(): 5 points (exact degree 9)
Quadrilaterals (2D tensor products):
- GaussLegendre{1}(): 1×1 = 1 point
- GaussLegendre{2}(): 2×2 = 4 points (standard Quad4)
- GaussLegendre{3}(): 3×3 = 9 points (standard Quad9)
- GaussLegendre{4}(): 4×4 = 16 points
- GaussLegendre{5}(): 5×5 = 25 points
Hexahedra (3D tensor products):
- GaussLegendre{1}(): 1×1×1 = 1 point
- GaussLegendre{2}(): 2×2×2 = 8 points (standard Hex8)
- GaussLegendre{3}(): 3×3×3 = 27 points (standard Hex27)
- GaussLegendre{4}(): 4×4×4 = 64 points
- GaussLegendre{5}(): 5×5×5 = 125 points
Legacy compatibility:
- Old Val{:GLSEG*}, Val{:GLQUAD*}, Val{:GLHEX*} symbols still work
- Uses _legacy_tensor_product() helper (kept old zip-based API)
Technical details:
- Rules 1-3: Hardcoded for optimal performance
- Rules 4-5: Generated from QUAD_DATA using ntuple for zero allocation
- All reference domains: [-1,1]^D
Net: +304 lines (programmatic generation → explicit implementations + docs)
This commit is contained in:
Jukka Aho
@@ -0,0 +1,344 @@
|
||||
# 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 tensor-product elements.
|
||||
|
||||
This file implements quadrature rules for elements formed by tensor products
|
||||
of 1D Gauss-Legendre rules:
|
||||
- **Segments** (1D): Direct 1D Gauss rules
|
||||
- **Quadrilaterals** (2D): Product of two 1D rules
|
||||
- **Hexahedra** (3D): Product of three 1D rules
|
||||
|
||||
# Available Rules
|
||||
|
||||
## Segments (1D)
|
||||
- `GaussLegendre{1}()`: 1-point, exact for linear
|
||||
- `GaussLegendre{2}()`: 2-point, exact for cubic
|
||||
- `GaussLegendre{3}()`: 3-point, exact for quintic
|
||||
- `GaussLegendre{4}()`: 4-point, exact for degree 7
|
||||
- `GaussLegendre{5}()`: 5-point, exact for degree 9
|
||||
|
||||
## Quadrilaterals (2D)
|
||||
- `GaussLegendre{1}()`: 1×1 = 1 point
|
||||
- `GaussLegendre{2}()`: 2×2 = 4 points (standard for bilinear)
|
||||
- `GaussLegendre{3}()`: 3×3 = 9 points (standard for biquadratic)
|
||||
- `GaussLegendre{4}()`: 4×4 = 16 points
|
||||
- `GaussLegendre{5}()`: 5×5 = 25 points
|
||||
|
||||
## Hexahedra (3D)
|
||||
- `GaussLegendre{1}()`: 1×1×1 = 1 point
|
||||
- `GaussLegendre{2}()`: 2×2×2 = 8 points (standard for trilinear)
|
||||
- `GaussLegendre{3}()`: 3×3×3 = 27 points (standard for triquadratic)
|
||||
- `GaussLegendre{4}()`: 4×4×4 = 64 points
|
||||
- `GaussLegendre{5}()`: 5×5×5 = 125 points
|
||||
|
||||
# Notes
|
||||
- Tensor products are formed from 1D Gauss-Legendre points in quaddata.jl
|
||||
- N points per dimension → N^D total points
|
||||
- Reference domains: Segment [-1,1], Quad [-1,1]², Hex [-1,1]³
|
||||
- All implementations use @inline for zero-allocation
|
||||
|
||||
See also: [`GaussLegendre`](@ref), [`QuadraturePoint`](@ref)
|
||||
"""
|
||||
|
||||
# ============================================================================
|
||||
# 1D Segments
|
||||
# ============================================================================
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Segment}, ::GaussLegendre{1})
|
||||
|
||||
1-point Gauss-Legendre rule for segment.
|
||||
Exact for linear polynomials (degree 1).
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Segment}, ::GaussLegendre{1,V}) where V
|
||||
return SVector(
|
||||
QuadraturePoint(Vec{1}(0.0), 2.0)
|
||||
)
|
||||
end
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Segment}, ::GaussLegendre{2})
|
||||
|
||||
2-point Gauss-Legendre rule for segment.
|
||||
Exact for cubic polynomials (degree 3).
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Segment}, ::GaussLegendre{2,V}) where V
|
||||
a = 1.0 / sqrt(3.0)
|
||||
return SVector(
|
||||
QuadraturePoint(Vec{1}(-a), 1.0),
|
||||
QuadraturePoint(Vec{1}( a), 1.0)
|
||||
)
|
||||
end
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Segment}, ::GaussLegendre{3})
|
||||
|
||||
3-point Gauss-Legendre rule for segment.
|
||||
Exact for quintic polynomials (degree 5).
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Segment}, ::GaussLegendre{3,V}) where V
|
||||
a = sqrt(3.0 / 5.0)
|
||||
return SVector(
|
||||
QuadraturePoint(Vec{1}(-a), 5.0/9.0),
|
||||
QuadraturePoint(Vec{1}(0.0), 8.0/9.0),
|
||||
QuadraturePoint(Vec{1}( a), 5.0/9.0)
|
||||
)
|
||||
end
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Segment}, ::GaussLegendre{4})
|
||||
|
||||
4-point Gauss-Legendre rule for segment.
|
||||
Exact for degree 7 polynomials.
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Segment}, ::GaussLegendre{4,V}) where V
|
||||
pts, wts = QUAD_DATA[4]
|
||||
return SVector(
|
||||
QuadraturePoint(Vec{1}(pts[1]), wts[1]),
|
||||
QuadraturePoint(Vec{1}(pts[2]), wts[2]),
|
||||
QuadraturePoint(Vec{1}(pts[3]), wts[3]),
|
||||
QuadraturePoint(Vec{1}(pts[4]), wts[4])
|
||||
)
|
||||
end
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Segment}, ::GaussLegendre{5})
|
||||
|
||||
5-point Gauss-Legendre rule for segment.
|
||||
Exact for degree 9 polynomials.
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Segment}, ::GaussLegendre{5,V}) where V
|
||||
pts, wts = QUAD_DATA[5]
|
||||
return SVector(
|
||||
QuadraturePoint(Vec{1}(pts[1]), wts[1]),
|
||||
QuadraturePoint(Vec{1}(pts[2]), wts[2]),
|
||||
QuadraturePoint(Vec{1}(pts[3]), wts[3]),
|
||||
QuadraturePoint(Vec{1}(pts[4]), wts[4]),
|
||||
QuadraturePoint(Vec{1}(pts[5]), wts[5])
|
||||
)
|
||||
end
|
||||
|
||||
# ============================================================================
|
||||
# 2D Quadrilaterals (tensor products)
|
||||
# ============================================================================
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Quadrilateral}, ::GaussLegendre{1})
|
||||
|
||||
1×1 = 1-point Gauss-Legendre rule for quadrilateral.
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Quadrilateral}, ::GaussLegendre{1,V}) where V
|
||||
return SVector(
|
||||
QuadraturePoint(Vec{2}(0.0, 0.0), 4.0)
|
||||
)
|
||||
end
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Quadrilateral}, ::GaussLegendre{2})
|
||||
|
||||
2×2 = 4-point Gauss-Legendre rule for quadrilateral.
|
||||
Standard rule for bilinear elements (Quad4).
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Quadrilateral}, ::GaussLegendre{2,V}) where V
|
||||
a = 1.0 / sqrt(3.0)
|
||||
return SVector(
|
||||
QuadraturePoint(Vec{2}(-a, -a), 1.0),
|
||||
QuadraturePoint(Vec{2}( a, -a), 1.0),
|
||||
QuadraturePoint(Vec{2}(-a, a), 1.0),
|
||||
QuadraturePoint(Vec{2}( a, a), 1.0)
|
||||
)
|
||||
end
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Quadrilateral}, ::GaussLegendre{3})
|
||||
|
||||
3×3 = 9-point Gauss-Legendre rule for quadrilateral.
|
||||
Standard rule for biquadratic elements (Quad9).
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Quadrilateral}, ::GaussLegendre{3,V}) where V
|
||||
a = sqrt(3.0 / 5.0)
|
||||
w1 = 5.0 / 9.0
|
||||
w2 = 8.0 / 9.0
|
||||
|
||||
return SVector(
|
||||
QuadraturePoint(Vec{2}(-a, -a), w1*w1),
|
||||
QuadraturePoint(Vec{2}(0.0, -a), w2*w1),
|
||||
QuadraturePoint(Vec{2}( a, -a), w1*w1),
|
||||
QuadraturePoint(Vec{2}(-a, 0.0), w1*w2),
|
||||
QuadraturePoint(Vec{2}(0.0, 0.0), w2*w2),
|
||||
QuadraturePoint(Vec{2}( a, 0.0), w1*w2),
|
||||
QuadraturePoint(Vec{2}(-a, a), w1*w1),
|
||||
QuadraturePoint(Vec{2}(0.0, a), w2*w1),
|
||||
QuadraturePoint(Vec{2}( a, a), w1*w1)
|
||||
)
|
||||
end
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Quadrilateral}, ::GaussLegendre{4})
|
||||
|
||||
4×4 = 16-point Gauss-Legendre rule for quadrilateral.
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Quadrilateral}, ::GaussLegendre{4,V}) where V
|
||||
pts, wts = QUAD_DATA[4]
|
||||
|
||||
result = ntuple(16) do i
|
||||
ix = (i - 1) % 4 + 1
|
||||
iy = div(i - 1, 4) + 1
|
||||
QuadraturePoint(Vec{2}(pts[ix], pts[iy]), wts[ix] * wts[iy])
|
||||
end
|
||||
|
||||
return SVector(result)
|
||||
end
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Quadrilateral}, ::GaussLegendre{5})
|
||||
|
||||
5×5 = 25-point Gauss-Legendre rule for quadrilateral.
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Quadrilateral}, ::GaussLegendre{5,V}) where V
|
||||
pts, wts = QUAD_DATA[5]
|
||||
|
||||
result = ntuple(25) do i
|
||||
ix = (i - 1) % 5 + 1
|
||||
iy = div(i - 1, 5) + 1
|
||||
QuadraturePoint(Vec{2}(pts[ix], pts[iy]), wts[ix] * wts[iy])
|
||||
end
|
||||
|
||||
return SVector(result)
|
||||
end
|
||||
|
||||
# ============================================================================
|
||||
# 3D Hexahedra (tensor products)
|
||||
# ============================================================================
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Hexahedron}, ::GaussLegendre{1})
|
||||
|
||||
1×1×1 = 1-point Gauss-Legendre rule for hexahedron.
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Hexahedron}, ::GaussLegendre{1,V}) where V
|
||||
return SVector(
|
||||
QuadraturePoint(Vec{3}(0.0, 0.0, 0.0), 8.0)
|
||||
)
|
||||
end
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Hexahedron}, ::GaussLegendre{2})
|
||||
|
||||
2×2×2 = 8-point Gauss-Legendre rule for hexahedron.
|
||||
Standard rule for trilinear elements (Hex8).
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Hexahedron}, ::GaussLegendre{2,V}) where V
|
||||
a = 1.0 / sqrt(3.0)
|
||||
|
||||
return SVector(
|
||||
QuadraturePoint(Vec{3}(-a, -a, -a), 1.0),
|
||||
QuadraturePoint(Vec{3}( a, -a, -a), 1.0),
|
||||
QuadraturePoint(Vec{3}(-a, a, -a), 1.0),
|
||||
QuadraturePoint(Vec{3}( a, a, -a), 1.0),
|
||||
QuadraturePoint(Vec{3}(-a, -a, a), 1.0),
|
||||
QuadraturePoint(Vec{3}( a, -a, a), 1.0),
|
||||
QuadraturePoint(Vec{3}(-a, a, a), 1.0),
|
||||
QuadraturePoint(Vec{3}( a, a, a), 1.0)
|
||||
)
|
||||
end
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Hexahedron}, ::GaussLegendre{3})
|
||||
|
||||
3×3×3 = 27-point Gauss-Legendre rule for hexahedron.
|
||||
Standard rule for triquadratic elements (Hex27).
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Hexahedron}, ::GaussLegendre{3,V}) where V
|
||||
a = sqrt(3.0 / 5.0)
|
||||
coords_1d = ((-a, 5.0/9.0), (0.0, 8.0/9.0), (a, 5.0/9.0))
|
||||
|
||||
result = ntuple(27) do i
|
||||
ix = (i - 1) % 3 + 1
|
||||
iy = div(i - 1, 3) % 3 + 1
|
||||
iz = div(i - 1, 9) + 1
|
||||
|
||||
x, wx = coords_1d[ix]
|
||||
y, wy = coords_1d[iy]
|
||||
z, wz = coords_1d[iz]
|
||||
|
||||
QuadraturePoint(Vec{3}(x, y, z), wx * wy * wz)
|
||||
end
|
||||
|
||||
return SVector(result)
|
||||
end
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Hexahedron}, ::GaussLegendre{4})
|
||||
|
||||
4×4×4 = 64-point Gauss-Legendre rule for hexahedron.
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Hexahedron}, ::GaussLegendre{4,V}) where V
|
||||
pts, wts = QUAD_DATA[4]
|
||||
|
||||
result = ntuple(64) do i
|
||||
ix = (i - 1) % 4 + 1
|
||||
iy = div(i - 1, 4) % 4 + 1
|
||||
iz = div(i - 1, 16) + 1
|
||||
QuadraturePoint(Vec{3}(pts[ix], pts[iy], pts[iz]), wts[ix] * wts[iy] * wts[iz])
|
||||
end
|
||||
|
||||
return SVector(result)
|
||||
end
|
||||
|
||||
"""
|
||||
get_quadrature_points(::Type{Hexahedron}, ::GaussLegendre{5})
|
||||
|
||||
5×5×5 = 125-point Gauss-Legendre rule for hexahedron.
|
||||
"""
|
||||
@inline function get_quadrature_points(::Type{Hexahedron}, ::GaussLegendre{5,V}) where V
|
||||
pts, wts = QUAD_DATA[5]
|
||||
|
||||
result = ntuple(125) do i
|
||||
ix = (i - 1) % 5 + 1
|
||||
iy = div(i - 1, 5) % 5 + 1
|
||||
iz = div(i - 1, 25) + 1
|
||||
QuadraturePoint(Vec{3}(pts[ix], pts[iy], pts[iz]), wts[ix] * wts[iy] * wts[iz])
|
||||
end
|
||||
|
||||
return SVector(result)
|
||||
end
|
||||
|
||||
# ============================================================================
|
||||
# Legacy symbol-based API (deprecated, kept for backwards compatibility)
|
||||
# ============================================================================
|
||||
|
||||
# Helper function for old API (generates zip of tuples)
|
||||
function _legacy_tensor_product(w::Tuple, p::Tuple, dim::Int)
|
||||
@assert length(w) == length(p)
|
||||
N = length(w)
|
||||
weights = Float64[]
|
||||
points = NTuple{dim, Float64}[]
|
||||
for i in CartesianIndices(ntuple(i -> 1:N, dim))
|
||||
push!(weights, prod(w[k] for k in Tuple(i)))
|
||||
push!(points, ntuple(k -> p[i[k]], dim))
|
||||
end
|
||||
return zip(Tuple(weights), Tuple(points))
|
||||
end
|
||||
|
||||
# Generate legacy symbol-based rules programmatically
|
||||
const _LEGACY_NAMES = [:GLSEG, :GLQUAD, :GLHEX]
|
||||
const _LEGACY_DIMS = (1, 2, 3)
|
||||
|
||||
for n in 1:length(QUAD_DATA)
|
||||
points, weights = QUAD_DATA[n]
|
||||
order = 2(n-1) + 1
|
||||
|
||||
for (dim, name) in zip(_LEGACY_DIMS, _LEGACY_NAMES)
|
||||
n_points = length(points)^dim
|
||||
quadname = QuoteNode(Symbol(string(name, n_points)))
|
||||
z = _legacy_tensor_product(weights, points, dim)
|
||||
|
||||
@eval begin
|
||||
get_quadrature_points(::Type{Val{$(quadname)}}) = $z
|
||||
get_order(::Type{Val{$(quadname)}}) = $order
|
||||
end
|
||||
end
|
||||
end
|
||||
@@ -1,40 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE
|
||||
|
||||
# Tensorial quadrature rules in 1, 2, 3 dimensions
|
||||
|
||||
function tensor_product(w::Tuple, p::Tuple, dim::Int)
|
||||
@assert length(w) == length(p)
|
||||
|
||||
N = length(w)
|
||||
weights = Float64[]
|
||||
points = NTuple{dim, Float64}[]
|
||||
for i in CartesianIndices(ntuple(i -> 1:N, dim))
|
||||
push!(weights, prod(w[k] for k in Tuple(i)))
|
||||
push!(points, ntuple(k -> p[i[k]], dim))
|
||||
end
|
||||
return zip(Tuple(weights), Tuple(points))
|
||||
end
|
||||
|
||||
|
||||
names = [:GLSEG, :GLQUAD, :GLHEX]
|
||||
names2 = ["segment", "quadrilateral", "hexahedron"]
|
||||
|
||||
for n in 1:length(QUAD_DATA)
|
||||
points, weights = QUAD_DATA[n]
|
||||
order = 2(n-1)+1
|
||||
for dim in (1, 2, 3)
|
||||
n_points = length(points)^dim
|
||||
quadname = QuoteNode(Symbol(string(names[dim], n_points)))
|
||||
z = tensor_product(weights, points, dim)
|
||||
@eval begin
|
||||
@doc """
|
||||
get_quadrature_points(::Type{Val{:$($(quadname))})
|
||||
|
||||
Gauss-Legendre quadrature, $($(n_points)) point rule on $($(names2[dim]))."""
|
||||
get_quadrature_points(::Type{Val{$(quadname)}}) = $z
|
||||
end
|
||||
@eval get_order(::Type{Val{$(quadname)}}) = $order
|
||||
end
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user