Calculate nodes and weights at GaussLegendre construction time (#149)

* Add fields and constructor

* Use new GaussLegendre fields

* Add Rope and GaussLegendre benchmarks

Author: mikeingold

benchmark/benchmarks.jl

@@ -45,6 +45,7 @@ spec = (
         line = Line(Point(0, 0, 0), Point(1, 1, 1)),
         plane = Plane(Point(0, 0, 0), Vec(0, 0, 1)),
         ray = Ray(Point(0, 0, 0), Vec(0, 0, 1)),
+        rope = Rope([Point(t, t, t) for t in 1:32]...),
         triangle = Triangle(Point(1, 0, 0), Point(0, 1, 0), Point(0, 0, 1)),
         tetrahedron = let
             a = Point(0, 3, 0)
@@ -64,6 +65,7 @@ SUITE["Specializations/Scalar GaussLegendre"] = let s = BenchmarkGroup()
     s["Line"] = @benchmarkable integral($spec.f_exp, $spec.g.line, $spec.rule_gl)
     s["Plane"] = @benchmarkable integral($spec.f_exp, $spec.g.plane, $spec.rule_gl)
     s["Ray"] = @benchmarkable integral($spec.f_exp, $spec.g.ray, $spec.rule_gl)
+    s["Rope"] = @benchmarkable integral($spec.f, $spec.g.rope, $spec.rule_gl)
     s["Triangle"] = @benchmarkable integral($spec.f, $spec.g.triangle, $spec.rule_gl)
     s["Tetrahedron"] = @benchmarkable integral($spec.f, $spec.g.tetrahedron, $spec.rule_gl)
     s
@@ -82,5 +84,14 @@ SUITE["Differentials"] = let s = BenchmarkGroup()
     s
 end
 
+############################################################################################
+#                                   Integration Rules
+###########################################################################################
+
+SUITE["Rules"] = let s = BenchmarkGroup()
+    s["GaussLegendre"] = @benchmarkable GaussLegendre($1000)
+    s
+end
+
 #tune!(SUITE)
 #run(SUITE, verbose=true)

src/integral.jl

@@ -75,17 +75,16 @@ function _integral(
     N = Meshes.paramdim(geometry)
 
     # Get Gauss-Legendre nodes and weights of type FP for a region [-1,1]ᴺ
-    xs, ws = FastGaussQuadrature.gausslegendre(rule.n)
-    xsFP = Iterators.map(FP, xs)
-    wsFP = Iterators.map(FP, ws)
-    weight_grid = Iterators.product(ntuple(Returns(wsFP), N)...)
-    node_grid = Iterators.product(ntuple(Returns(xsFP), N)...)
+    xs = Iterators.map(FP, rule.nodes)
+    ws = Iterators.map(FP, rule.weights)
+    weight_grid = Iterators.product(ntuple(Returns(ws), N)...)
+    node_grid = Iterators.product(ntuple(Returns(xs), N)...)
 
     # Domain transformation: x [-1,1] ↦ t [0,1]
     t(x) = (1 // 2) * x + (1 // 2)
 
     function integrand((weights, nodes))
-        # Transforms nodes/xs, store in an NTuple 
+        # ts = t.(nodes), but non-allocating
         ts = ntuple(i -> t(nodes[i]), length(nodes))
         # Integrand function
         prod(weights) * f(geometry(ts...)) * differential(geometry, ts, diff_method)

src/integration_rules.jl

@@ -33,6 +33,10 @@ e.g. `length(geometry)/λ`.
 """
 struct GaussLegendre <: IntegrationRule
     n::Int64
+    nodes::Vector{Float64}
+    weights::Vector{Float64}
+
+    GaussLegendre(n::Int64) = new(n, FastGaussQuadrature.gausslegendre(n)...)
 end
 
 """