diff --git a/src/quadrature/gl_triangles.jl b/src/quadrature/gl_triangles.jl new file mode 100644 index 0000000..2475750 --- /dev/null +++ b/src/quadrature/gl_triangles.jl @@ -0,0 +1,245 @@ +# 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 triangular elements. + +This file implements quadrature rules for 2D triangular reference elements +in the parametric domain [0,1]² with constraint ξ + η ≤ 1. + +# Available Rules +- `GaussLegendre{1}()`: 1-point rule (centroid), exact for linear +- `GaussLegendre{2}()`: 3-point rule, exact for quadratic (default variant) +- `GaussLegendre{2,:B}()`: 3-point alternative (midpoint variant) +- `GaussLegendre{3}()`: 4-point rule, exact for cubic (default variant) +- `GaussLegendre{3,:B}()`: 4-point alternative (has negative weight!) +- `GaussLegendre{4}()`: 6-point rule, exact for quartic +- `GaussLegendre{5}()`: 7-point rule, exact for quintic +- `GaussLegendre{6}()`: 12-point rule, exact for 6th degree + +# Notes +- Reference triangle has vertices at (0,0), (1,0), (0,1) +- All weights sum to 0.5 (area of reference triangle) +- Variant :B rules exist for historical reasons, default is recommended +- Negative weights in GLTRI4B may cause issues, use default GaussLegendre{3}() + +See also: [`GaussLegendre`](@ref), [`QuadraturePoint`](@ref) +""" + +# ============================================================================ +# GaussLegendre{1}: 1-point rule (centroid) +# ============================================================================ + +""" + get_quadrature_points(::Type{Triangle}, ::GaussLegendre{1}) + +1-point Gauss-Legendre rule for triangle (centroid). +Exact for linear polynomials (degree 1). +""" +@inline function get_quadrature_points(::Type{Triangle}, ::GaussLegendre{1,V}) where V + return SVector( + QuadraturePoint(Vec{2}(1/3, 1/3), 0.5) + ) +end + +# ============================================================================ +# GaussLegendre{2}: 3-point rules (two variants) +# ============================================================================ + +""" + get_quadrature_points(::Type{Triangle}, ::GaussLegendre{2}) + +3-point Gauss-Legendre rule for triangle (default variant). +Exact for quadratic polynomials (degree 2). +Points located at edge centers of sub-triangles. +""" +@inline function get_quadrature_points(::Type{Triangle}, ::GaussLegendre{2,:default}) + w = 1/6 + return SVector( + QuadraturePoint(Vec{2}(2/3, 1/6), w), + QuadraturePoint(Vec{2}(1/6, 2/3), w), + QuadraturePoint(Vec{2}(1/6, 1/6), w) + ) +end + +""" + get_quadrature_points(::Type{Triangle}, ::GaussLegendre{2,:B}) + +3-point Gauss-Legendre rule for triangle (variant B). +Exact for quadratic polynomials (degree 2). +Points located at edge midpoints. +""" +@inline function get_quadrature_points(::Type{Triangle}, ::GaussLegendre{2,:B}) + w = 1/6 + return SVector( + QuadraturePoint(Vec{2}(0.0, 1/2), w), + QuadraturePoint(Vec{2}(1/2, 0.0), w), + QuadraturePoint(Vec{2}(1/2, 1/2), w) + ) +end + +# ============================================================================ +# GaussLegendre{3}: 4-point rules (two variants) +# ============================================================================ + +""" + get_quadrature_points(::Type{Triangle}, ::GaussLegendre{3}) + +4-point Gauss-Legendre rule for triangle (default variant). +Exact for cubic polynomials (degree 3). +All weights positive. +""" +@inline function get_quadrature_points(::Type{Triangle}, ::GaussLegendre{3,:default}) + return SVector( + QuadraturePoint(Vec{2}(0.15505102572168219, 0.17855872826361642), 0.15902069087198858), + QuadraturePoint(Vec{2}(0.64494897427831781, 0.07503111022260812), 0.09097930912801142), + QuadraturePoint(Vec{2}(0.15505102572168219, 0.66639024601470139), 0.15902069087198858), + QuadraturePoint(Vec{2}(0.64494897427831781, 0.28001991549907407), 0.09097930912801142) + ) +end + +""" + get_quadrature_points(::Type{Triangle}, ::GaussLegendre{3,:B}) + +4-point Gauss-Legendre rule for triangle (variant B). +Exact for cubic polynomials (degree 3). + +**WARNING:** This rule has one negative weight (-27/96), which may cause +numerical issues. Use default variant instead unless specifically needed. +""" +@inline function get_quadrature_points(::Type{Triangle}, ::GaussLegendre{3,:B}) + return SVector( + QuadraturePoint(Vec{2}(1/3, 1/3), -27/96), # Negative weight! + QuadraturePoint(Vec{2}(1/5, 1/5), 25/96), + QuadraturePoint(Vec{2}(1/5, 3/5), 25/96), + QuadraturePoint(Vec{2}(3/5, 1/5), 25/96) + ) +end + +# ============================================================================ +# GaussLegendre{4}: 6-point rule +# ============================================================================ + +""" + get_quadrature_points(::Type{Triangle}, ::GaussLegendre{4}) + +6-point Gauss-Legendre rule for triangle. +Exact for quartic polynomials (degree 4). +""" +@inline function get_quadrature_points(::Type{Triangle}, ::GaussLegendre{4,V}) where V + P1 = 0.11169079483905 + P2 = 0.0549758718227661 + A = 0.445948490915965 + B = 0.091576213509771 + + return SVector( + QuadraturePoint(Vec{2}(B, B), P2), + QuadraturePoint(Vec{2}(1.0 - 2.0*B, B), P2), + QuadraturePoint(Vec{2}(B, 1.0 - 2.0*B), P2), + QuadraturePoint(Vec{2}(A, 1.0 - 2.0*A), P1), + QuadraturePoint(Vec{2}(A, A), P1), + QuadraturePoint(Vec{2}(1.0 - 2.0*A, A), P1) + ) +end + +# ============================================================================ +# GaussLegendre{5}: 7-point rule +# ============================================================================ + +""" + get_quadrature_points(::Type{Triangle}, ::GaussLegendre{5}) + +7-point Gauss-Legendre rule for triangle. +Exact for quintic polynomials (degree 5). +Includes centroid point. +""" +@inline function get_quadrature_points(::Type{Triangle}, ::GaussLegendre{5,V}) where V + A = 0.470142064105115 + B = 0.101286507323456 + P1 = 0.066197076394253 + P2 = 0.062969590272413 + + return SVector( + QuadraturePoint(Vec{2}(1/3, 1/3), 9/80), + QuadraturePoint(Vec{2}(A, A), P1), + QuadraturePoint(Vec{2}(1.0 - 2.0*A, A), P1), + QuadraturePoint(Vec{2}(A, 1.0 - 2.0*A), P1), + QuadraturePoint(Vec{2}(B, B), P2), + QuadraturePoint(Vec{2}(1.0 - 2.0*B, B), P2), + QuadraturePoint(Vec{2}(B, 1.0 - 2.0*B), P2) + ) +end + +# ============================================================================ +# GaussLegendre{6}: 12-point rule +# ============================================================================ + +""" + get_quadrature_points(::Type{Triangle}, ::GaussLegendre{6}) + +12-point Gauss-Legendre rule for triangle. +Exact for 6th degree polynomials. +""" +@inline function get_quadrature_points(::Type{Triangle}, ::GaussLegendre{6,V}) where V + A = 0.063089014491502 + B = 0.249286745170910 + C = 0.310352451033785 + D = 0.053145049844816 + P1 = 0.025422453185103 + P2 = 0.058393137863189 + P3 = 0.041425537809187 + + return SVector( + QuadraturePoint(Vec{2}(A, A), P1), + QuadraturePoint(Vec{2}(1.0 - 2.0*A, A), P1), + QuadraturePoint(Vec{2}(A, 1.0 - 2.0*A), P1), + QuadraturePoint(Vec{2}(B, B), P2), + QuadraturePoint(Vec{2}(1.0 - 2.0*B, B), P2), + QuadraturePoint(Vec{2}(B, 1.0 - 2.0*B), P2), + QuadraturePoint(Vec{2}(C, D), P3), + QuadraturePoint(Vec{2}(D, C), P3), + QuadraturePoint(Vec{2}(1.0 - C - D, C), P3), + QuadraturePoint(Vec{2}(1.0 - C - D, D), P3), + QuadraturePoint(Vec{2}(C, 1.0 - C - D), P3), + QuadraturePoint(Vec{2}(D, 1.0 - C - D), P3) + ) +end + +# ============================================================================ +# Legacy symbol-based API (deprecated, kept for backwards compatibility) +# ============================================================================ + +# These map old :GLTRI symbols to new API +get_quadrature_points(::Type{Val{:GLTRI1}}) = + get_quadrature_points(Triangle, GaussLegendre{1}()) + +get_quadrature_points(::Type{Val{:GLTRI3}}) = + get_quadrature_points(Triangle, GaussLegendre{2}()) + +get_quadrature_points(::Type{Val{:GLTRI3B}}) = + get_quadrature_points(Triangle, GaussLegendre{2,:B}()) + +get_quadrature_points(::Type{Val{:GLTRI4}}) = + get_quadrature_points(Triangle, GaussLegendre{3}()) + +get_quadrature_points(::Type{Val{:GLTRI4B}}) = + get_quadrature_points(Triangle, GaussLegendre{3,:B}()) + +get_quadrature_points(::Type{Val{:GLTRI6}}) = + get_quadrature_points(Triangle, GaussLegendre{4}()) + +get_quadrature_points(::Type{Val{:GLTRI7}}) = + get_quadrature_points(Triangle, GaussLegendre{5}()) + +get_quadrature_points(::Type{Val{:GLTRI12}}) = + get_quadrature_points(Triangle, GaussLegendre{6}()) + +# Order queries (deprecated) +get_order(::Type{Val{:GLTRI1}}) = 1 +get_order(::Type{Val{:GLTRI3}}) = 2 +get_order(::Type{Val{:GLTRI3B}}) = 2 +get_order(::Type{Val{:GLTRI4}}) = 3 +get_order(::Type{Val{:GLTRI4B}}) = 3 +get_order(::Type{Val{:GLTRI6}}) = 4 +get_order(::Type{Val{:GLTRI7}}) = 5 +get_order(::Type{Val{:GLTRI12}}) = 6 diff --git a/src/quadrature/gltri.jl b/src/quadrature/gltri.jl deleted file mode 100644 index 7f659f5..0000000 --- a/src/quadrature/gltri.jl +++ /dev/null @@ -1,131 +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 triangular elements - -""" Gauss-Legendre quadrature, 1 point rule on triangle. """ -function get_quadrature_points(::Type{Val{:GLTRI1}}) - weights = (0.5, ) - points = ((1.0/3.0, 1.0/3.0), ) - return zip(weights, points) -end - -function get_order(::Type{Val{:GLTRI1}}) - return 1 -end - -""" Gauss-Legendre quadrature, 3 point rule on triangle. """ -function get_quadrature_points(::Type{Val{:GLTRI3}}) - weights = (1.0/6.0, 1.0/6.0, 1.0/6.0) - points = ( - (2.0/3.0, 1.0/6.0), - (1.0/6.0, 2.0/3.0), - (1.0/6.0, 1.0/6.0) - ) - return zip(weights, points) -end - -function get_order(::Type{Val{:GLTRI3}}) - return 2 -end - -""" Gauss-Legendre quadrature, 3 point rule on triangle. """ -function get_quadrature_points(::Type{Val{:GLTRI3B}}) - weights = (1.0/6.0, 1.0/6.0, 1.0/6.0) - points = ( - (0.0, 1.0/2.0), - (1.0/2.0, 0.0), - (1.0/2.0, 1.0/2.0) - ) - return zip(weights, points) -end - -function get_order(::Type{Val{:GLTRI3B}}) - return 2 -end - -""" Gauss-Legendre quadrature, 4 point rule on triangle. """ -function get_quadrature_points(::Type{Val{:GLTRI4}}) - weights = ( - 1.5902069087198858469718450103758e-01, - 9.0979309128011415302815498962418e-02, - 1.5902069087198858469718450103758e-01, - 9.0979309128011415302815498962418e-02) - points = ( - (1.5505102572168219018027159252941e-01, 1.7855872826361642311703513337422e-01), - (6.4494897427831780981972840747059e-01, 7.5031110222608118177475598324603e-02), - (1.5505102572168219018027159252941e-01, 6.6639024601470138670269327409637e-01), - (6.4494897427831780981972840747059e-01, 2.8001991549907407200279599420481e-01)) - return zip(weights, points) -end - -function get_order(::Type{Val{:GLTRI4}}) - return 3 -end - -""" Gauss-Legendre quadrature, 4 point rule on triangle. """ -function get_quadrature_points(::Type{Val{:GLTRI4B}}) - weights = (-27.0/96.0, 25.0/96.0, 25.0/96.0, 25.0/96.0) - points = ( - (1.0/3.0, 1.0/3.0), - (1.0/5.0, 1.0/5.0), - (1.0/5.0, 3.0/5.0), - (3.0/5.0, 1.0/5.0)) - return zip(weights, points) -end - -function get_order(::Type{Val{:GLTRI4B}}) - return 3 -end - -""" Gauss-Legendre quadrature, 6 point rule on triangle. """ -function get_quadrature_points(::Type{Val{:GLTRI6}}) - P1 = 0.11169079483905 - P2 = 0.0549758718227661 - A = 0.445948490915965 - B = 0.091576213509771 - weights = (P2, P2, P2, P1, P1, P1) - points = ((B, B), (1.0-2.0*B, B), (B, 1.0-2.0*B), - (A, 1.0-2*A), (A, A), (1.0 - 2.0*A, A)) - return zip(weights, points) -end - -function get_order(::Type{Val{:GLTRI6}}) - return 4 -end - -""" Gauss-Legendre quadrature, 7 point rule on triangle. """ -function get_quadrature_points(::Type{Val{:GLTRI7}}) - A = 0.470142064105115 - B = 0.101286507323456 - P1 = 0.066197076394253 - P2 = 0.062969590272413 - weights = (9/80, P1, P1, P1, P2, P2, P2) - points = ((1/3, 1/3), (A, A), (1.0-2.0*A, A), (A, 1.0-2.0*A), - (B, B), (1.0-2.0*B, B), (B, 1.0-2.0*B)) - return zip(weights, points) -end - -function get_order(::Type{Val{:GLTRI7}}) - return 5 -end - -""" Gauss-Legendre quadrature, 12 point rule on triangle. """ -function get_quadrature_points(::Type{Val{:GLTRI12}}) - A = 0.063089014491502 - B = 0.249286745170910 - C = 0.310352451033785 - D = 0.053145049844816 - P1 = 0.025422453185103 - P2 = 0.058393137863189 - P3 = 0.041425537809187 - weights = (P1, P1, P1, P2, P2, P2, P3, P3, P3, P3, P3, P3) - points = ((A, A), (1.0-2.0*A, A), (A, 1.0-2.0*A), (B, B), (1.0-2.0*B, B), - (B, 1.0-2.0*B), (C, D), (D, C), (1.0-C-D, C), (1.0-C-D, D), - (C, 1.0-C-D), (D, 1.0-C-D)) - return zip(weights, points) -end - -function get_order(::Type{Val{:GLTRI12}}) - return 6 -end