From 55d589c2694a3caa9f3f7b3152a3cf67c00fbf8c Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 21 Oct 2024 19:18:25 -0400 Subject: [PATCH 1/6] Bump min version of CliffordNumbers (new feature) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9a2d3922..d8b4f938 100644 --- a/Project.toml +++ b/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" From ea1e4b721b858e6afdd23eea89aef7949b0305d4 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 21 Oct 2024 19:33:53 -0400 Subject: [PATCH 2/6] New KVector and units methods --- src/utils.jl | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index e0117df8..7beef005 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -8,11 +8,6 @@ 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." @@ -20,11 +15,15 @@ function _error_unsupported_combination(geometry, rule) 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) From 721425518e235ff74ca25959a03f1f69d84cf135 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 21 Oct 2024 19:42:27 -0400 Subject: [PATCH 3/6] Use new KVector method to condense differential code --- src/differentiation.jl | 17 ++++------------- 1 file changed, 4 insertions(+), 13 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 9fcb53ac..2783c5ff 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -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, 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 + J = Iterators.map(_KVector, jacobian(geometry, ts)) + return LinearAlgebra.norm(foldl(∧, J)) end From 5de4385aa4e5744b481de1b94db97cfbd4c3cb26 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 21 Oct 2024 19:43:34 -0400 Subject: [PATCH 4/6] Condense redundant type params --- src/differentiation.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index 2783c5ff..a40db5a6 100644 --- a/src/differentiation.jl +++ b/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.")) @@ -89,9 +89,9 @@ 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{M, CRS}, ts::V -) where {M, CRS, G <: Meshes.Geometry{M, CRS}, V <: Union{AbstractVector, Tuple}} +) where {V <: Union{AbstractVector, Tuple}} J = Iterators.map(_KVector, jacobian(geometry, ts)) return LinearAlgebra.norm(foldl(∧, J)) end From 5da201d7e3bcf3c6e47234ee251649095262c551 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Mon, 21 Oct 2024 19:55:03 -0400 Subject: [PATCH 5/6] Remove remnant type params --- src/differentiation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation.jl b/src/differentiation.jl index a40db5a6..94731c90 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -89,7 +89,7 @@ Return the magnitude of the differential element (length, area, volume, etc) of the parametric function for `geometry` at arguments `ts`. """ function differential( - geometry::Meshes.Geometry{M, CRS}, + geometry::Meshes.Geometry, ts::V ) where {V <: Union{AbstractVector, Tuple}} J = Iterators.map(_KVector, jacobian(geometry, ts)) From 50c12c4c310a3fd8a488c25d2e3e671279ee2676 Mon Sep 17 00:00:00 2001 From: Michael Ingold Date: Wed, 30 Oct 2024 07:29:07 -0400 Subject: [PATCH 6/6] Update kvector test --- test/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/utils.jl b/test/utils.jl index 45453a50..5b25c11e 100644 --- a/test/utils.jl +++ b/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")