Use newly-supported KVector Unitful extension to simplify differential

Project.toml

@@ -13,7 +13,7 @@ QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
-CliffordNumbers = "0.1.4"
+CliffordNumbers = "0.1.9"
 CoordRefSystems = "0.12, 0.13, 0.14, 0.15"
 FastGaussQuadrature = "1"
 HCubature = "1.5"

src/differentiation.jl

@@ -14,10 +14,10 @@ finite-difference approximation with step size `ε`.
 - `ε`: step size to use for the finite-difference approximation
 """
 function jacobian(
-        geometry::G,
+        geometry::Meshes.Geometry,
         ts::V;
         ε = 1e-6
-) where {G <: Meshes.Geometry, V <: Union{AbstractVector, Tuple}}
+) where {V <: Union{AbstractVector, Tuple}}
     Dim = Meshes.paramdim(geometry)
     if Dim != length(ts)
         throw(ArgumentError("ts must have same number of dimensions as geometry."))
@@ -85,22 +85,13 @@ end
 """
     differential(geometry, ts)
 
-Calculate the differential element (length, area, volume, etc) of the parametric
-function for `geometry` at arguments `ts`.
+Return the magnitude of the differential element (length, area, volume, etc) of
+the parametric function for `geometry` at arguments `ts`.
 """
 function differential(
-        geometry::G,
+        geometry::Meshes.Geometry,
         ts::V
-) where {M, CRS, G <: Meshes.Geometry{M, CRS}, V <: Union{AbstractVector, Tuple}}
-    # Calculate the Jacobian, convert Vec -> KVector
-    J = jacobian(geometry, ts)
-    J_kvecs = Iterators.map(_kvector, J)
-
-    # Extract units from Geometry type
-    Dim = Meshes.paramdim(geometry)
-    units = _units(geometry)^Dim
-
-    # Return norm of the exterior products
-    element = foldl(∧, J_kvecs)
-    return LinearAlgebra.norm(element) * units
+) where {V <: Union{AbstractVector, Tuple}}
+    J = Iterators.map(_KVector, jacobian(geometry, ts))
+    return LinearAlgebra.norm(foldl(∧, J))
 end

src/utils.jl

@@ -8,23 +8,22 @@ function _gausslegendre(T, n)
     return T.(xs), T.(ws)
 end
 
-# Extract the length units used by the CRS of a Geometry
-function _units(g::Meshes.Geometry{M, CRS}) where {M, CRS}
-    return Unitful.unit(CoordRefSystems.lentype(CRS))
-end
-
 # Common error message structure
 function _error_unsupported_combination(geometry, rule)
     msg = "Integrating a $geometry using a $rule rule not supported."
     throw(ArgumentError(msg))
 end
 
 ################################################################################
-#                        CliffordNumbers Interface
+#                        CliffordNumbers and Units
 ################################################################################
 
-# Meshes.Vec -> ::CliffordNumber.KVector
-function _kvector(v::Meshes.Vec{Dim, T}) where {Dim, T}
+# Meshes.Vec -> Unitful.Quantity{CliffordNumber.KVector}
+function _KVector(v::Meshes.Vec{Dim, T}) where {Dim, T}
     ucoords = Iterators.map(Unitful.ustrip, v.coords)
-    return CliffordNumbers.KVector{1, VGA(Dim)}(ucoords...)
+    return CliffordNumbers.KVector{1, VGA(Dim)}(ucoords...) * _units(v)
 end
+
+# Extract the length units used by the CRS of a Geometry
+_units(::Meshes.Geometry{M, CRS}) where {M, CRS} = Unitful.unit(CoordRefSystems.lentype(CRS))
+_units(::Meshes.Vec{Dim, T}) where {Dim, T} = Unitful.unit(T)

test/utils.jl

@@ -3,7 +3,7 @@
 
     # _kvector
     v = Meshes.Vec(3, 4)
-    @test norm(MeshIntegrals._kvector(v)) ≈ 5.0
+    @test norm(MeshIntegrals._KVector(v)) ≈ 5.0u"m"
 
     # _units
     p = Point(1.0u"cm", 2.0u"mm", 3.0u"m")