feat: Consolidate FEMQuad.jl into JuliaFEM (quadrature rules)

Consolidated entire FEMQuad.jl package (436 lines) into src/quadrature/:
- quaddata.jl: Quadrature data definitions
- glquad.jl: 2D quadrilateral Gauss-Legendre rules
- gltri.jl: 2D triangle Gauss-Legendre rules (131 lines)
- gltet.jl: 3D tetrahedron Gauss-Legendre rules
- glwed.jl: 3D wedge Gauss-Legendre rules
- glpyr.jl: 3D pyramid Gauss-Legendre rules

Changes:
- Created src/quadrature.jl as main include file
- Removed 'import FEMQuad' from JuliaFEM.jl
- Updated integrate.jl: FEMQuad.get_quadrature_points → get_quadrature_points
- Added export add_element! (was missing)

Result:
- ✅ JuliaFEM loads successfully
- ✅ Integration points work correctly
- ✅ 5 tests still passing (no regression)
- ✅ One less vendor package dependency

Next: Continue consolidating vendor packages

src/JuliaFEM.jl

@@ -117,7 +117,7 @@ using Tensors  # For basis functions (Vec type)
 import Calculus  # For symbolic differentiation in basis generation
 
 import FEMSparse
-import FEMQuad  # Still using vendor FEMQuad for now
+# import FEMQuad  # Consolidated into src/quadrature.jl
 
 # Note: Consolidating FEMBase and FEMBasis into JuliaFEM
 # Previously: @reexport using FEMBase
@@ -141,6 +141,9 @@ include("basis/nurbs_surface.jl")
 include("basis/nurbs_solid.jl")
 include("basis/math.jl")
 
+# Quadrature rules (consolidated from FEMQuad.jl)
+include("quadrature.jl")
+
 # Consolidate FEMBase.jl into src/ (Phase 1 continued)
 # Order matters: fields → types → sparse → elements → integrate → problems → assembly
 include("fields/fields.jl")          # Field system (DCTI, DVTI, etc.)
@@ -238,7 +241,7 @@ export DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment
 export FieldProblem, BoundaryProblem, Problem, Node, Element, Assembly
 export Poi1, Seg2, Seg3, Tri3, Tri6, Tri7, Quad4, Quad8, Quad9,
    Tet4, Tet10, Pyr5, Wedge6, Wedge15, Hex8, Hex20, Hex27
-export update!, add_elements!, get_unknown_field_name, add!,
+export update!, add_element!, add_elements!, get_unknown_field_name, add!,
    is_field_problem, is_boundary_problem, get_gdofs,
    initialize!, get_integration_points, group_by_element_type,
    get_unknown_field_dimension, get_connectivity

src/elements/integrate.jl

@@ -37,13 +37,13 @@ for (E, R) in integration_rule_mapping
         if isequal(i, 1)
             code = quote
                 function get_integration_points(element::$E)
-                    return FEMQuad.get_quadrature_points($P)
+                    return get_quadrature_points($P)
                 end
             end
         else
             code = quote
                 function get_integration_points(element::$E, ::Type{$order})
-                    return FEMQuad.get_quadrature_points($P)
+                    return get_quadrature_points($P)
                 end
             end
         end

src/quadrature.jl

@@ -0,0 +1,62 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE
+#
+# Gaussian-Legendre quadrature rules consolidated from FEMQuad.jl
+
+# Quadrature data
+include("quadrature/quaddata.jl")
+
+# Element-specific quadrature rules
+include("quadrature/glquad.jl")  # 2D quadrilaterals
+include("quadrature/gltri.jl")   # 2D triangles
+include("quadrature/gltet.jl")   # 3D tetrahedrons
+include("quadrature/glwed.jl")   # 3D wedges
+include("quadrature/glpyr.jl")   # 3D pyramids
+
+"""
+    get_rule(order::Int, rules::Symbol...)
+
+Get the first quadrature rule that meets the required order.
+"""
+function get_rule(order::Int, rules::Vararg{Symbol})
+    for rule in rules
+        if get_order(Val{rule}) >= order
+            return rule
+        end
+    end
+    @warn("No accurate rule enough found, picking last.", order, rules)
+    return rules[end]
+end
+
+"""
+    integrate_1d(f::Function, rule::Symbol)
+
+Integrate a 1D function using the specified quadrature rule.
+"""
+function integrate_1d(f::Function, rule::Symbol)
+    points = get_quadrature_points(Val{rule})
+    result = sum(w*f(ip) for (w, ip) in points)
+    return result
+end
+
+"""
+    integrate_2d(f::Function, rule::Symbol)
+
+Integrate a 2D function using the specified quadrature rule.
+"""
+function integrate_2d(f::Function, rule::Symbol)
+    points = get_quadrature_points(Val{rule})
+    result = sum(w*f(ip) for (w, ip) in points)
+    return result
+end
+
+"""
+    integrate_3d(f::Function, rule::Symbol)
+
+Integrate a 3D function using the specified quadrature rule.
+"""
+function integrate_3d(f::Function, rule::Symbol)
+    points = get_quadrature_points(Val{rule})
+    result = sum(w*f(ip) for (w, ip) in points)
+    return result
+end

src/quadrature/FEMQuad.jl

@@ -0,0 +1,48 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE
+
+module FEMQuad
+
+# Gaussian-Legendre quadratures
+include("quaddata.jl")
+include("glquad.jl")
+# Gaussian-Legendre quadratures in 2d triangles
+include("gltri.jl")
+# Gaussian-Legendre quadratures in 3d hexahedrons
+include("gltet.jl")
+# Gaussian-Legendre quadratures in 3d wedges
+include("glwed.jl")
+# Gaussian-Legendre quadratures in 3d pyramid
+include("glpyr.jl")
+
+function get_rule(order::Int, rules::Vararg{Symbol})
+    for rule in rules
+        if get_order(Val{rule}) >= order
+            return rule
+        end
+    end
+    @warn("No accurate rule enough found, picking last.", order, rules)
+    return rules[end]
+end
+
+function integrate_1d(f::Function, rule::Symbol)
+    points = get_quadrature_points(Val{rule})
+    result = sum(w*f(ip) for (w, ip) in points)
+    return result
+end
+
+function integrate_2d(f::Function, rule::Symbol)
+    points = get_quadrature_points(Val{rule})
+    result = sum(w*f(ip) for (w, ip) in points)
+    return result
+end
+
+function integrate_3d(f::Function, rule::Symbol)
+    points = get_quadrature_points(Val{rule})
+    result = sum(w*f(ip) for (w, ip) in points)
+    return result
+end
+
+export get_quadrature_points
+
+end

src/quadrature/glpyr.jl

@@ -0,0 +1,45 @@
+# 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

src/quadrature/glquad.jl

@@ -0,0 +1,40 @@
+# 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
+

src/quadrature/gltet.jl

@@ -0,0 +1,67 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE
+
+### Gauss quadrature rules for tetrahedrons
+
+""" Gauss-Legendre quadrature, 1 point rule on tetrahedron. """
+function get_quadrature_points(::Type{Val{:GLTET1}})
+    weights = (1.0/6.0, )
+    points = ((1.0/4.0, 1.0/4.0, 1.0/4.0), )
+    return zip(weights, points)
+end
+
+function get_order(::Type{Val{:GLTET1}})
+    return 1
+end
+
+""" Gauss-Legendre quadrature, 4 point rule on tetrahedron. """
+function get_quadrature_points(::Type{Val{:GLTET4}})
+    a = (5.0+3.0*sqrt(5.0))/20.0
+    b = (5.0-sqrt(5.0))/20.0
+    w = 1.0/24.0
+    weights = (w, w, w, w)
+    points = ((a, b, b), (b, a, b), (b, b, a), (b, b, b))
+    return zip(weights, points)
+end
+
+function get_order(::Type{Val{:GLTET4}})
+    return 2
+end
+
+""" Gauss-Legendre quadrature, 5 point rule on tetrahedron. """
+function get_quadrature_points(::Type{Val{:GLTET5}})
+    a = 1.0/4.0
+    b = 1.0/6.0
+    c = 1.0/2.0
+    weights = (-2.0/15.0, 3.0/40.0, 3.0/40.0, 3.0/40.0, 3.0/40.0)
+    points = ((a, a, a), (b, b, b), (b, b, c), (b, c, b), (c, b, b))
+    return zip(weights, points)
+end
+
+function get_order(::Type{Val{:GLTET5}})
+    return 3
+end
+
+""" Gauss-Legendre quadrature, 15 point rule on tetrahedron. """
+function get_quadrature_points(::Type{Val{:GLTET15}})
+    a = 1.0/4.0
+    b1 = 1.0/34.0*(7.0 + sqrt(15.0))
+    b2 = 1.0/34.0*(7.0 - sqrt(15.0))
+    c1 = 1.0/34.0*(13.0 - 3.0*sqrt(15.0))
+    c2 = 1.0/34.0*(13.0 + 3.0*sqrt(15.0))
+    d = 1.0/20.0*(5.0 - sqrt(15.0))
+    f = 1.0/20.0*(5.0 + sqrt(15.0))
+    w1 = 8.0/405.0
+    w2 = (2665.0 - 14.0*sqrt(15.0))/226800.0
+    w3 = (2665.0 + 14.0*sqrt(15.0))/226800.0
+    w4 = 5.0/567.0
+    weights = (w1, w2, w2, w2, w2, w3, w3, w3, w3, w4, w4, w4, w4, w4, w4)
+    points = ((a, a, a), (b1, b1, b1), (b1, b1, c1), (b1, c1, b1), (c1, b1, b1),
+              (b2, b2, b2), (b2, b2, c2), (b2, c2, b2), (c2, b2, b2), (d, d, f),
+              (d, f, d), (f, d, d), (d, f, f), (f, d, f), (f, f, d))
+    return zip(weights, points)
+end
+
+function get_order(::Type{Val{:GLTET15}})
+    return 4
+end

src/quadrature/gltri.jl

@@ -0,0 +1,131 @@
+# 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

src/quadrature/glwed.jl

@@ -0,0 +1,93 @@
+# 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

src/quadrature/quaddata.jl

@@ -0,0 +1,12 @@
+const QUAD_DATA = [
+[(0.0,),
+(2.0,)],
+[(-0.5773502691896258, 0.5773502691896258,),
+(1.0, 1.0,)],
+[(-0.7745966692414834, 0.0, 0.7745966692414834,),
+(0.5555555555555556, 0.8888888888888888, 0.5555555555555556,)],
+[(-0.8611363115940526, -0.3399810435848563, 0.3399810435848563, 0.8611363115940526,),
+(0.34785484513745385, 0.6521451548625462, 0.6521451548625462, 0.34785484513745385,)],
+[(-0.906179845938664, -0.5384693101056831, 0.0, 0.5384693101056831, 0.906179845938664,),
+(0.23692688505618908, 0.47862867049936647, 0.5688888888888889, 0.47862867049936647, 0.23692688505618908,)],
+]