consolidate some of the standard Gauss-Legendre rules a bit

bin/genquaddata.jl

@@ -0,0 +1,15 @@
+# Generates a tablaue of weights and points for Gauss Legendre quadrature
+# on a line
+using FastGaussQuadrature
+
+const MAX_ORDER = 5
+open(joinpath(@__DIR__, "..", "src", "quaddata.jl"), "w") do f
+    println(f, "const QUAD_DATA = [")
+    for g in 1:MAX_ORDER
+        points, weights = gausslegendre(g)
+        println(f, "[(", join(points, ", "), ",),")
+        println(f, "(", join(weights, ", "), ",)],")
+    end
+    println(f, "]")
+end
+

docs/src/api.md

@@ -69,8 +69,8 @@ These rules are get from 1d quadratures by using tensor production.
 FEMQuad.get_quadrature_points(::Type{Val{:GLHEX1}})
 FEMQuad.get_quadrature_points(::Type{Val{:GLHEX8}})
 FEMQuad.get_quadrature_points(::Type{Val{:GLHEX27}})
-FEMQuad.get_quadrature_points(::Type{Val{:GLHEX81}})
-FEMQuad.get_quadrature_points(::Type{Val{:GLHEX243}})
+FEMQuad.get_quadrature_points(::Type{Val{:GLHEX64}})
+FEMQuad.get_quadrature_points(::Type{Val{:GLHEX125}})
 ```
 
 ### Gauss-Legendre rules in prismatic domain

src/FEMQuad.jl

@@ -3,15 +3,12 @@
 
 module FEMQuad
 
-# Gaussian-Legendre quadratures in 1d segments
-include("glseg.jl")
-# Gaussian-Legendre quadratures in 2d quadrangles
+# 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("glhex.jl")
-# Gaussian-Legendre quadratures in 3d tetrahedrons
 include("gltet.jl")
 # Gaussian-Legendre quadratures in 3d wedges
 include("glwed.jl")

src/glquad.jl

@@ -1,90 +1,40 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE
 
-### Gauss quadrature rules for quadrilaterals
+# Tensorial quadrature rules in 1, 2, 3 dimensions
+
+function tensor_product(w::Tuple, p::Tuple, dim::Int)
+    @assert length(w) == length(p)
 
-function tensor_product_2d(w::Tuple, p::Tuple)
     N = length(w)
     weights = Float64[]
-    points = Tuple{Float64, Float64}[]
-    for i in 1:N
-        for j in 1:N
-            push!(weights, w[i]*w[j])
-            push!(points, (p[i], p[j]))
-        end
+    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(weights, points)
-end
-
-""" Gauss-Legendre quadrature, 1 point rule on quadrilateral. """
-function get_quadrature_points(::Type{Val{:GLQUAD1}})
-    w = (2.0, )
-    p = (0.0, )
-    return tensor_product_2d(w, p)
-end
-
-function get_order(::Type{Val{:GLQUAD1}})
-    return 1
-end
-
-""" Gauss-Legendre quadrature, 4 point rule on quadrilateral. """
-function get_quadrature_points(::Type{Val{:GLQUAD4}})
-    w = (1.0, 1.0)
-    p = (-sqrt(1.0/3.0), sqrt(1.0/3.0))
-    return tensor_product_2d(w, p)
-end
-
-function get_order(::Type{Val{:GLQUAD4}})
-    return 3
+    return zip(Tuple(weights), Tuple(points))
 end
 
-""" Gauss-Legendre quadrature, 9 point rule on quadrilateral. """
-function get_quadrature_points(::Type{Val{:GLQUAD9}})
-    w = (5.0/9.0, 8.0/9.0, 5.0/9.0)
-    p = (-sqrt(3.0/5.0), 0.0, sqrt(3.0/5.0))
-    return tensor_product_2d(w, p)
-end
 
-function get_order(::Type{Val{:GLQUAD9}})
-    return 5
-end
+names = [:GLSEG, :GLQUAD, :GLHEX]
+names2 = ["segment", "quadrilateral", "hexahedron"]
 
-""" Gauss-Legendre quadrature, 16 point rule on quadrilateral. """
-function get_quadrature_points(::Type{Val{:GLQUAD16}})
-    w = (
-         1.0/36.0*(18.0-sqrt(30.0)),
-         1.0/36.0*(18.0+sqrt(30.0)),
-         1.0/36.0*(18.0+sqrt(30.0)),
-         1.0/36.0*(18.0-sqrt(30.0)))
-    p = (
-        -sqrt(3.0/7.0 + 2.0/7.0*sqrt(6.0/5.0)),
-        -sqrt(3.0/7.0 - 2.0/7.0*sqrt(6.0/5.0)),
-         sqrt(3.0/7.0 - 2.0/7.0*sqrt(6.0/5.0)),
-         sqrt(3.0/7.0 + 2.0/7.0*sqrt(6.0/5.0)))
-    return tensor_product_2d(w, p)
-end
+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))})
 
-function get_order(::Type{Val{:GLQUAD16}})
-    return 7
-end
-
-""" Gauss-Legendre quadrature, 25 point rule on quadrilateral. """
-function get_quadrature_points(::Type{Val{:GLQUAD25}})
-    w = (
-        1.0/900.0*(322.0-13.0*sqrt(70.0)),
-        1.0/900.0*(322.0+13.0*sqrt(70.0)),
-        128.0/225.0,
-        1.0/900.0*(322.0+13.0*sqrt(70.0)),
-        1.0/900.0*(322.0-13.0*sqrt(70.0)))
-    p = (
-        -1.0/3.0*sqrt(5.0 + 2.0*sqrt(10.0/7.0)),
-        -1.0/3.0*sqrt(5.0 - 2.0*sqrt(10.0/7.0)),
-         0.0,
-         1.0/3.0*sqrt(5.0 - 2.0*sqrt(10.0/7.0)),
-         1.0/3.0*sqrt(5.0 + 2.0*sqrt(10.0/7.0)))
-    return tensor_product_2d(w, p)
+            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
 
-function get_order(::Type{Val{:GLQUAD25}})
-    return 9
-end

src/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,)],
+]

test/runtests.jl

@@ -6,8 +6,7 @@ using Test
 include("../docs/make.jl")
 
 @testset "FEMQuad.jl" begin
-    @testset "test_glseg" begin include("test_glseg.jl") end
-    @testset "test_glhex" begin include("test_glhex.jl") end
+    @testset "test_glquad" begin include("test_glquad.jl") end
     @testset "test_gltri" begin include("test_gltri.jl") end
     @testset "test_gltet" begin include("test_gltet.jl") end
     @testset "test_glwed" begin include("test_glwed.jl") end