Add GaussLegendre quadrature type

src/AbstractTypes.jl

@@ -64,6 +64,25 @@ struct GaussLobattoLegendre <: AbstractQuadratureType
   end
 end
 
+"""
+$(TYPEDEF)
+"""
+struct GaussLegendre <: AbstractQuadratureType
+  cell_degree::Int
+  surf_degree::Int
+
+  function GaussLegendre(degree::Int)
+    @assert degree > 0 "Quadrature degree must be greater than zero"
+    new(degree, degree)
+  end
+
+  function GaussLegendre(cell_degree::Int, surf_degree::Int)
+    @assert cell_degree > 0 "Cell quadrature degree must be greater than zero"
+    @assert surf_degree > 0 "Surface quadrature degree must be greater than zero"
+    new(cell_degree, surf_degree)
+  end
+end
+
 # methods to define for quadrature types
 """
 $(TYPEDSIGNATURES)

src/ReferenceFiniteElements.jl

@@ -80,7 +80,9 @@ export RaviartThomas
 export Serendipity
 
 # quadrature types
+export AbstractQuadratureType
 export GaussLobattoLegendre
+export GaussLegendre
 
 # reference fe
 export ReferenceFE

src/elements/Edge.jl

@@ -141,6 +141,31 @@ num_cell_dofs(::Edge{Lagrange, PD}) where PD = PD + 1
 num_dofs_on_boundary(::Edge{Lagrange, PD}, ::Int) where PD = 1
 num_interior_dofs(::Edge{Lagrange, PD}) where PD = PD < 2 ? 0 : PD - 1
 
+function cell_quadrature_points_and_weights(e::AbstractEdge, q_rule::GaussLegendre)
+    ξs, ws = gausslegendre(cell_quadrature_degree(q_rule))
+
+    if e.shifted
+        ξs .= (ξs .+ 1.) ./ 2.
+        ws .= ws ./ 2.
+    end
+
+    return reshape(ξs, 1, length(ξs)), ws
+end
+
+function surface_quadrature_points_and_weights(e::AbstractEdge, ::GaussLegendre)
+    if e.shifted
+        x_min = 0.
+    else
+        x_min = -1.
+    end
+
+    ξs = zeros(1, 1, 2)
+    ξs[1, 1, 1] = x_min
+    ξs[1, 1, 2] = 1.
+    ws = ones(1, 2)
+    return ξs, ws
+end
+
 function cell_quadrature_points_and_weights(e::AbstractEdge, q_rule::GaussLobattoLegendre)
     ξs, ws = gausslegendre(cell_quadrature_degree(q_rule))
 

src/elements/Hex.jl

@@ -144,6 +144,49 @@ function num_interior_dofs(::Hex{Lagrange, PD}) where PD
     end
 end
 
+function cell_quadrature_points_and_weights(e::AbstractHex, q_rule::GaussLegendre)
+    ξs, ws = cell_quadrature_points_and_weights(boundary_element(boundary_element(e, 0), 0), q_rule)
+    n = length(ws)
+    ξ_return = Matrix{eltype(ξs)}(undef, 3, n * n * n)
+    w_return = Vector{eltype(ξs)}(undef, n * n * n)
+    for (q, ξ) in enumerate(Base.Iterators.product(ξs, ξs, ξs))
+        ξ_return[1, q] = ξ[1]
+        ξ_return[2, q] = ξ[2]
+        ξ_return[3, q] = ξ[3]
+    end
+    for (q, w) in enumerate(Base.Iterators.product(ws, ws, ws))
+        w_return[q] = w[1] * w[2] * w[3]
+    end
+    return ξ_return, w_return
+end
+
+function surface_quadrature_points_and_weights(e::AbstractHex, q_rule::GaussLegendre)
+    ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule)
+
+    ξ_return = zeros(3, length(ws), 6)
+    w_return = zeros(length(ws), 6)
+
+    ξ_return[1:2, :, 1] .= ξs
+    ξ_return[3, :, 1]   .= -1.
+    ξ_return[1, :, 2]   .= 1.
+    ξ_return[2:3, :, 2] .= ξs
+    ξ_return[1:2, :, 3] .= ξs
+    ξ_return[3, :, 3]   .= 1.
+    ξ_return[1, :, 4]   .= -1.
+    ξ_return[2:3, :, 4] .= ξs
+    ξ_return[1, :, 5]   .= ξs[1, :]
+    ξ_return[2, :, 5]   .= -1.
+    ξ_return[3, :, 5]   .= ξs[2, :]
+    ξ_return[1, :, 6]   .= ξs[1, :]
+    ξ_return[2, :, 6]   .= 1.
+    ξ_return[3, :, 6]   .= ξs[2, :]
+
+    for n in 1:6
+        w_return[:, n] .= ws
+    end
+    return ξ_return, w_return
+end
+
 function cell_quadrature_points_and_weights(e::AbstractHex, q_rule::GaussLobattoLegendre)
     ξs, ws = cell_quadrature_points_and_weights(boundary_element(boundary_element(e, 0), 0), q_rule)
     ξ_return = Matrix{eltype(ξs)}(undef, 3, length(ξs) * length(ξs) * length(ξs) * length(ξs))

src/elements/Quad.jl

@@ -75,6 +75,41 @@ function num_interior_dofs(::Quad{Lagrange, PD}) where PD
     end
 end
 
+function cell_quadrature_points_and_weights(e::AbstractQuad, q_rule::GaussLegendre)
+    ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule)
+    ξ_return = Matrix{eltype(ξs)}(undef, 2, length(ws) * length(ws))
+    w_return = Vector{eltype(ξs)}(undef, length(ws) * length(ws))
+    for (q, ξ) in enumerate(Base.Iterators.product(ξs, ξs))
+        ξ_return[1, q] = ξ[1]
+        ξ_return[2, q] = ξ[2]
+    end
+    for (q, w) in enumerate(Base.Iterators.product(ws, ws))
+        w_return[q] = w[1] * w[2]
+    end
+    return ξ_return, w_return
+end
+
+function surface_quadrature_points_and_weights(e::AbstractQuad, q_rule::GaussLegendre)
+    ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule)
+
+    ξ_return = zeros(2, length(ws), 4)
+    w_return = zeros(length(ws), 4)
+
+    ξ_return[1, :, 1] .= ξs[1, :]
+    ξ_return[2, :, 1] .= -1.
+    ξ_return[1, :, 2] .= 1.
+    ξ_return[2, :, 2] .= ξs[1, :]
+    ξ_return[1, :, 3] .= ξs[1, :]
+    ξ_return[2, :, 3] .= 1.
+    ξ_return[1, :, 4] .= -1.
+    ξ_return[2, :, 4] .= ξs[1, :]
+
+    for n in 1:4
+        w_return[:, n] .= ws
+    end
+    return ξ_return, w_return
+end
+
 function cell_quadrature_points_and_weights(e::AbstractQuad, q_rule::GaussLobattoLegendre)
     # ξs, ws = gausslegendre(cell_quadrature_degree(e))
     ξs, ws = cell_quadrature_points_and_weights(boundary_element(e, 0), q_rule)

src/elements/Tet.jl

@@ -178,6 +178,43 @@ end
 num_cell_dofs(::Tet{Lagrange, PD}) where PD = (PD + 1) * (PD + 2) * (PD + 3) ÷ 6
 num_interior_dofs(::Tet{Lagrange, PD}) where PD = PD < 4 ? 0 : (PD - 1) * (PD - 2) * (PD - 3) ÷ 6
 
+function cell_quadrature_points_and_weights(::AbstractTet, q_rule::GaussLegendre)
+    deg = cell_quadrature_degree(q_rule)
+    if deg == 1
+        # 1-point centroid rule (degree 1)
+        ξs = Matrix{Float64}(undef, 3, 1)
+        ξs[:, 1] = [1. / 4., 1. / 4., 1. / 4.]
+        ws = [1. / 6.]
+    elseif deg == 2
+        # 4-point symmetric rule (degree 2)
+        s = sqrt(5.0)
+        a = (5. + 3. * s) / 20.
+        b = (5. - s) / 20.
+        ξs = Matrix{Float64}(undef, 3, 4)
+        ξs[:, 1] = [b, b, b]
+        ξs[:, 2] = [a, b, b]
+        ξs[:, 3] = [b, a, b]
+        ξs[:, 4] = [b, b, a]
+        ws = [1. / 24., 1. / 24., 1. / 24., 1. / 24.]
+    elseif deg == 3
+        # 5-point rule (degree 3)
+        ξs = Matrix{Float64}(undef, 3, 5)
+        ξs[:, 1] = [1. / 4., 1. / 4., 1. / 4.]
+        ξs[:, 2] = [1. / 6., 1. / 6., 1. / 6.]
+        ξs[:, 3] = [1. / 6., 1. / 6., 1. / 2.]
+        ξs[:, 4] = [1. / 6., 1. / 2., 1. / 6.]
+        ξs[:, 5] = [1. / 2., 1. / 6., 1. / 6.]
+        ws = [-2. / 15., 3. / 40., 3. / 40., 3. / 40., 3. / 40.]
+    else
+        @assert false "GaussLegendre degree 1 through 3 supported for Tet."
+    end
+    return ξs, ws
+end
+
+function surface_quadrature_points_and_weights(e::AbstractTet, q_rule::GaussLegendre)
+    return surface_quadrature_points_and_weights(e, GaussLobattoLegendre(cell_quadrature_degree(q_rule), surface_quadrature_degree(q_rule)))
+end
+
 function cell_quadrature_points_and_weights(e::AbstractTet, q_rule::GaussLobattoLegendre)
     if cell_quadrature_degree(q_rule) == 1
         ξs = Matrix{Float64}(undef, 3, 1)

src/elements/Tri.jl

@@ -1,3 +1,35 @@
+function cell_quadrature_points_and_weights(::AbstractTri, q_rule::GaussLegendre)
+    deg = cell_quadrature_degree(q_rule)
+    if deg == 1
+        # 1-point centroid rule (degree 1)
+        ξs = Matrix{Float64}(undef, 2, 1)
+        ξs[:, 1] = [1. / 3., 1. / 3.]
+        ws = [0.5]
+    elseif deg == 2
+        # 3-point rule (degree 2)
+        ξs = Matrix{Float64}(undef, 2, 3)
+        ξs[:, 1] = [1. / 6., 1. / 6.]
+        ξs[:, 2] = [4. / 6., 1. / 6.]
+        ξs[:, 3] = [1. / 6., 4. / 6.]
+        ws = [1. / 6., 1. / 6., 1. / 6.]
+    elseif deg == 3
+        # 4-point rule (degree 3): centroid + 3 edge midpoints
+        ξs = Matrix{Float64}(undef, 2, 4)
+        ξs[:, 1] = [1. / 3., 1. / 3.]
+        ξs[:, 2] = [1. / 5., 3. / 5.]
+        ξs[:, 3] = [3. / 5., 1. / 5.]
+        ξs[:, 4] = [1. / 5., 1. / 5.]
+        ws = [-27. / 96., 25. / 96., 25. / 96., 25. / 96.]
+    else
+        @assert false "GaussLegendre degree 1 through 3 supported for Tri."
+    end
+    return ξs, ws
+end
+
+function surface_quadrature_points_and_weights(e::AbstractTri, q_rule::GaussLegendre)
+    return surface_quadrature_points_and_weights(e, GaussLobattoLegendre(cell_quadrature_degree(q_rule), surface_quadrature_degree(q_rule)))
+end
+
 function cell_quadrature_points_and_weights(::AbstractTri, q_rule::GaussLobattoLegendre)
     if cell_quadrature_degree(q_rule) == 1
       ξs = Matrix{Float64}(undef, 2, 1)