diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 202323e..40ef45c 100644
--- a/src/JuliaFEM.jl
+++ b/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
diff --git a/src/elements/integrate.jl b/src/elements/integrate.jl
index 571da66..1f0957e 100644
--- a/src/elements/integrate.jl
+++ b/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
diff --git a/src/quadrature.jl b/src/quadrature.jl
new file mode 100644
index 0000000..6a1801a
--- /dev/null
+++ b/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
diff --git a/src/quadrature/FEMQuad.jl b/src/quadrature/FEMQuad.jl
new file mode 100644
index 0000000..84f1b57
--- /dev/null
+++ b/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
diff --git a/src/quadrature/glpyr.jl b/src/quadrature/glpyr.jl
new file mode 100644
index 0000000..f49da64
--- /dev/null
+++ b/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
diff --git a/src/quadrature/glquad.jl b/src/quadrature/glquad.jl
new file mode 100644
index 0000000..26addde
--- /dev/null
+++ b/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
+
diff --git a/src/quadrature/gltet.jl b/src/quadrature/gltet.jl
new file mode 100644
index 0000000..96b7c8e
--- /dev/null
+++ b/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
diff --git a/src/quadrature/gltri.jl b/src/quadrature/gltri.jl
new file mode 100644
index 0000000..7f659f5
--- /dev/null
+++ b/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
diff --git a/src/quadrature/glwed.jl b/src/quadrature/glwed.jl
new file mode 100644
index 0000000..2bb4ef0
--- /dev/null
+++ b/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
diff --git a/src/quadrature/quaddata.jl b/src/quadrature/quaddata.jl
new file mode 100644
index 0000000..dec1cf8
--- /dev/null
+++ b/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,)],
+]