Consolidate differential element calculations (#36)

* Consolidate paramfactor functions into a common differential function

* Fix typo

* Fix bug by vectorizing scalar arg

* Tweak differential and add docstring

* Bump version

Author: mikeingold

Project.toml

@@ -1,7 +1,7 @@
 name = "MeshIntegrals"
 uuid = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47"
 authors = ["Mike Ingold <[email protected]>"]
-version = "0.11.3"
+version = "0.11.4"
 
 [deps]
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"

src/integral_line.jl

@@ -16,13 +16,8 @@ function _integral_1d(
     t(x) = T(1/2) * x + T(1/2)
     point(x) = geometry(t(x))
 
-    function paramfactor(x)
-        J = jacobian(geometry, T[t(x)])
-        return norm(J[1])
-    end
-
     # Integrate f along the geometry and apply a domain-correction factor for [-1,1] ↦ [0, 1]
-    integrand((w,x)) = w * f(point(x)) * paramfactor(x)
+    integrand((w,x)) = w * f(point(x)) * differential(geometry, [t(x)])
     return T(1/2) * sum(integrand, zip(ws, xs))
 end
 
@@ -32,13 +27,7 @@ function _integral_1d(
     settings::GaussKronrod
 )
     T = coordtype(geometry)
-    
-    function paramfactor(t)
-        J = jacobian(geometry, T[t])
-        return norm(J[1])
-    end
-
-    integrand(t) = f(geometry(t)) * paramfactor(t)
+    integrand(t) = f(geometry(t)) * differential(geometry, [t])
     return QuadGK.quadgk(integrand, T(0), T(1); settings.kwargs...)[1]
 end
 
@@ -48,13 +37,7 @@ function _integral_1d(
     settings::HAdaptiveCubature
 )
     T = coordtype(geometry)
-
-    function paramfactor(t)
-        J = jacobian(geometry, t)
-        return norm(J[1])
-    end
-
-    integrand(t) = f(geometry(t[1])) * paramfactor(t)
+    integrand(t) = f(geometry(t[1])) * differential(geometry, t)
     return HCubature.hcubature(integrand, T[0], T[1]; settings.kwargs...)[1]
 end
 

src/integral_surface.jl

@@ -17,13 +17,7 @@ function _integral_2d(
     v(x) = T(1/2)*x + T(1/2)
     point(xi,xj) = geometry2d(u(xi), v(xj))
 
-    function paramfactor(xi, xj)
-        uv = [u(xi), v(xj)]
-        J = jacobian(geometry2d, uv)
-        return norm(J[1] × J[2])
-    end
-
-    integrand(((wi,wj), (xi,xj))) = wi * wj * f(point(xi,xj)) * paramfactor(xi,xj)
+    integrand(((wi,wj), (xi,xj))) = wi * wj * f(point(xi,xj)) * differential(geometry2d, [u(xi), v(xj)])
 
     return T(1/4) .* sum(integrand, zip(wws,xxs))
 end
@@ -33,12 +27,7 @@ function _integral_2d(
     geometry2d::G,
     settings::GaussKronrod
 ) where {Dim, T, G<:Meshes.Geometry{Dim,T}}
-    function paramfactor(uv)
-        J = jacobian(geometry2d, uv)
-        return norm(J[1] × J[2])
-    end
-
-    integrand(u,v) = f(geometry2d(u,v)) * paramfactor([u,v])
+    integrand(u,v) = f(geometry2d(u,v)) * differential(geometry2d, [u,v])
     innerintegral(v) = QuadGK.quadgk(u -> integrand(u,v), T(0), T(1); settings.kwargs...)[1]
     return QuadGK.quadgk(v -> innerintegral(v), T(0), T(1); settings.kwargs...)[1]
 end
@@ -48,12 +37,7 @@ function _integral_2d(
     geometry2d::G,
     settings::HAdaptiveCubature
 ) where {Dim, T, G<:Meshes.Geometry{Dim,T}}
-    function paramfactor(uv)
-        J = jacobian(geometry2d, uv)
-        return norm(J[1] × J[2])
-    end
-
-    integrand(uv) = paramfactor(uv) * f(geometry2d(uv...))
+    integrand(uv) = differential(geometry2d, uv) * f(geometry2d(uv...))
     return hcubature(integrand, T[0,0], T[1,1]; settings.kwargs...)[1]
 end
 

src/integral_volume.jl

@@ -19,14 +19,9 @@ function _integral_3d(
 
     point(stu) = geometry3d(stu[1], stu[2], stu[3])
 
-    function paramfactor(stu)
-        J = jacobian(geometry3d, stu)
-        return abs((J[1] × J[2]) ⋅ J[3])
-    end
-
     function integrand(((wi,wj,wk), (xi,xj,xk)))
         stu = [s(xi),t(xj),u(xk)]
-        wi * wj * wk * f(point(stu)) * paramfactor(stu)
+        wi * wj * wk * f(point(stu)) * differential(geometry3d, stu)
     end
 
     return T(1/8) .* sum(integrand, zip(wws,xxs))
@@ -37,12 +32,7 @@ function _integral_3d(
     geometry3d::G,
     settings::HAdaptiveCubature
 ) where {Dim, T, G<:Meshes.Geometry{Dim,T}}
-    function paramfactor(ts)
-        J = jacobian(geometry3d, ts)
-        return abs((J[1] × J[2]) ⋅ J[3])
-    end
-
-    integrand(ts) = paramfactor(ts) * f(geometry3d(ts...))
+    integrand(ts) = differential(geometry3d, ts) * f(geometry3d(ts...))
     return hcubature(integrand, zeros(T,3), ones(T,3); settings.kwargs...)[1]
 end
 

src/utils.jl

@@ -15,18 +15,18 @@ function jacobian(
 ) where {T<:AbstractFloat}
     # Get the partial derivative along the n'th axis via finite difference approximation
     #   where ts is the current parametric position (εv is a reusable buffer)
-    function ∂r_∂tn!(εv, ts, n)
+    function ∂ₙr!(εv, ts, n)
         if ts[n] < T(0.01)
-            return ∂r_∂tn_right!(εv, ts, n)
+            return ∂ₙr_right!(εv, ts, n)
         elseif T(0.99) < ts[n]
-            return ∂r_∂tn_left!(εv, ts, n)
+            return ∂ₙr_left!(εv, ts, n)
         else
-            return ∂r_∂tn_central!(εv, ts, n)
+            return ∂ₙr_central!(εv, ts, n)
         end
     end
 
     # Central finite difference
-    function ∂r_∂tn_central!(εv, ts, n)
+    function ∂ₙr_central!(εv, ts, n)
         εv .= T(0)
         εv[n] = ε
         a = ts - εv
@@ -35,7 +35,7 @@ function jacobian(
     end
 
     # Left finite difference
-    function ∂r_∂tn_left!(εv, ts, n)
+    function ∂ₙr_left!(εv, ts, n)
         εv .= T(0)
         εv[n] = ε
         a = ts - εv
@@ -44,7 +44,7 @@ function jacobian(
     end
 
     # Right finite difference
-    function ∂r_∂tn_right!(εv, ts, n)
+    function ∂ₙr_right!(εv, ts, n)
         εv .= T(0)
         εv[n] = ε
         a = ts
@@ -55,8 +55,29 @@ function jacobian(
     # Allocate a re-usable ε vector
     εv = similar(ts)
 
-    ∂r_∂tn(n) = ∂r_∂tn!(εv,ts,n)
-    return map(∂r_∂tn, 1:length(ts))
+    ∂ₙr(n) = ∂ₙr!(εv,ts,n)
+    return map(∂ₙr, 1:length(ts))
+end
+
+"""
+    differential(geometry, ts)
+
+Calculate the differential element (length, area, volume, etc) of the parametric
+function for `geometry` at arguments `ts`.
+"""
+function differential(geometry, ts::AbstractVector)
+    J = jacobian(geometry, ts)
+
+    # TODO generalize this with geometric algebra, e.g.: norm(foldl(∧, J))
+    if length(J) == 1
+        return norm(J[1])
+    elseif length(J) == 2
+        return norm(J[1] × J[2])
+    elseif length(J) == 3
+        return abs((J[1] × J[2]) ⋅ J[3])
+    else
+        error("Not implemented yet. Please report this as an Issue on GitHub.")
+    end
 end
 
 """