Fix incompatibility with recent Meshes changes (#42)

* Add branch TODO and update Meshes compat

* Remove some unnecessary param types

* Restructure such that T is an optional first argument

* Bump version and uodate readme with status

* Replace T of DataType param with FP of Type T

* Remove checking for MethodError with user function

* Remove more MethodError checks

* Continue implementing new structure

* Reorganize arguments (optional FP is now last), various bugfixes

* Update tests

* Bugfixes

* Fix Box parameterization

* Fix Geometry parameterization

* Bugfixes

* Bugfix in Plane vector normalization

* Bugfixes in Line and Ray vector normalization

* Add shim code to handle Unitful types in HCubature methods

* Bugfix missing using statement

* Bugfix typo

* Bugfix for broadcastable f returns

* More broadcasting required

* Update HCubature line and volume integrals

* Implement new HCubature technique for BezierCurve specialization

* Update CylinderSurface to outsource disk integrations

* Update HCubature methods with more generalized parametric coord generators

* Make a generic HCubature method and point 1D integrals to it

* Remove old commented out code and reorganize

* Update integral function signatures, split alias functions into new file

* Convert all _integral_nd methods to point toward new generalized version

* Bugfix to method signatures

* Bugfix add missing include for new file

* Bugfix typos

* Add Ring and Rope specializations for HCubature methods to avoid ambiguities

* Bugfix internal structure of Point type was changed

* Bugfix tests taking exp of unitful quantity

* Implement support for HCubature on CylinderSurface

* Update docs and formatting

* Remove unused type parameter

* Bugfix add Unitful dep to tests

* Bugfix copy paste leftover

* Update Line, Ray, and Plane methods to ensure unitful domain handling via CRS

* Bugfixes

* Bugfixes

* Remove branch-specific notes

Project.toml

@@ -1,19 +1,23 @@
 name = "MeshIntegrals"
 uuid = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47"
 authors = ["Mike Ingold <[email protected]>"]
-version = "0.11.7"
+version = "0.12.0"
 
 [deps]
+CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
+CoordRefSystems = "0.10"
 FastGaussQuadrature = "1"
 HCubature = "1.5"
 julia = "1.6"
 LinearAlgebra = "1"
-Meshes = "0.36 - 0.42"
+Meshes = "0.36 - 0.47"
 QuadGK = "2"
+Unitful = "1"

README.md

@@ -1,23 +1,15 @@
 # MeshIntegrals.jl
 
 This package implements methods for numerically-computing integrals over geometric polytopes
-from [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl) using:
+from [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl) using the following `::IntegrationAlgorithms`:
 - Gauss-Legendre quadrature rules from [**FastGaussQuadrature.jl**](https://github.com/JuliaApproximation/FastGaussQuadrature.jl): `GaussLegendre(n)`
 - H-adaptive Gauss-Kronrod quadrature rules from [**QuadGK.jl**](https://github.com/JuliaMath/QuadGK.jl): `GaussKronrod(kwargs...)`
 - H-adaptive cubature rules from [**HCubature.jl**](https://github.com/JuliaMath/HCubature.jl): `HAdaptiveCubature(kwargs...)`
 
 Functions available:
-- `integral(f, geometry, ::IntegrationAlgorithm)`: integrates a function `f` over a domain defined by `geometry` using a particular
+- `integral(f, ::Geometry, ::IntegrationAlgorithm)`: integrates a function `f` over a domain defined by `geometry` using a particular `::IntegrationAlgorithm`
 - `lineintegral`, `surfaceintegral`, and `volumeintegral` are available as aliases for `integral` that first verify that `geometry` has the appropriate number of parametric dimensions
 
-Integration methods are also compatible with:
-- Meshes.jl geometries with **Unitful.jl** coordinate types, e.g. `Point(1.0u"m", 2.0u"m")`
-- Any `f(::Meshes.Point{Dim,<:Real})` that maps to a value type that **QuadGK.jl** can integrate, including:
-    - Real or complex-valued scalars
-    - Real or complex-valued vectors
-    - Dimensionful scalars or vectors from Unitful.jl
-    - Dimensionful scalars or vectors from DynamicQuantities.jl
-
 # Example Usage
 
 ```julia

src/MeshIntegrals.jl

@@ -1,6 +1,8 @@
 module MeshIntegrals
+    using CoordRefSystems
     using LinearAlgebra
     using Meshes
+    using Unitful
 
     import FastGaussQuadrature
     import HCubature
@@ -10,6 +12,7 @@ module MeshIntegrals
     export jacobian, derivative, unitdirection
 
     include("integral.jl")
+    include("integral_aliases.jl")
     include("integral_line.jl")
     include("integral_surface.jl")
     include("integral_volume.jl")

src/integral.jl

@@ -48,90 +48,103 @@ end
 ################################################################################
 
 """
-    integral(f, geometry, algorithm::IntegrationAlgorithm)
+    integral(f, geometry)
+    integral(f, geometry, algorithm)
+    integral(f, geometry, algorithm, FP)
 
 Numerically integrate a given function `f(::Point)` over the domain defined by
-a `geometry` using a particular `integration algorithm`.
+a `geometry` using a particular `integration algorithm` with floating point
+precision of type `FP`.
 """
+function integral end
+
+# If only f and geometry are specified, default to HAdaptiveCubature
 function integral(
     f::F,
-    geometry::G,
-    settings::I=HAdaptiveCubature()
-) where {F<:Function, Dim, T, G<:Meshes.Geometry{Dim,T}, I<:IntegrationAlgorithm}
-    # Validate that the provided function has an appropriate f(::Point{Dim,T}) method
-    _validate_integrand(f, Dim, T)
-
-    # Run the appropriate integral type
-    dim_param = paramdim(geometry)
-    if dim_param == 1
-        return _integral_1d(f, geometry, settings)
-    elseif dim_param == 2
-        return _integral_2d(f, geometry, settings)
-    elseif dim_param == 3
-        return _integral_3d(f, geometry, settings)
-    end
+    geometry::G
+) where {F<:Function, G<:Meshes.Geometry}
+    _integral(f, geometry, HAdaptiveCubature())
 end
 
-################################################################################
-#                    Line, Surface, and Volume Aliases
-################################################################################
-
-"""
-    lineintegral(f, geometry, algorithm::IntegrationAlgorithm=GaussKronrod)
+# If algorithm is HAdaptiveCubature, use the generalized n-dim solution
+function integral(
+    f::F,
+    geometry::G,
+    settings::HAdaptiveCubature
+) where {F<:Function, G<:Meshes.Geometry}
+    _integral(f, geometry, settings)
+end
 
-Numerically integrate a given function `f(::Point)` along a line-like `geometry`
-using a particular `integration algorithm`.
+# If algorithm is HAdaptiveCubature, and FP specified, use the generalized n-dim solution
+function integral(
+    f::F,
+    geometry::G,
+    settings::HAdaptiveCubature,
+    FP::Type{T}
+) where {F<:Function, G<:Meshes.Geometry, T<:AbstractFloat}
+    _integral(f, geometry, settings, FP)
+end
 
-Algorithm types available:
-- GaussKronrod
-- GaussLegendre
-- HAdaptiveCubature
-"""
-function lineintegral(
+# If algorithm is not HAdaptiveCubature, specialize on number of dimensions
+function integral(
     f::F,
     geometry::G,
-    settings::I=GaussKronrod()
+    settings::I
 ) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm}
-    if (dim = paramdim(geometry)) == 1
-        return integral(f, geometry, settings)
-    else
-        error("Performing a line integral on a geometry with $dim parametric dimensions not supported.")
+    # Run the appropriate integral type
+    Dim = Meshes.paramdim(geometry)
+    if Dim == 1
+        return _integral_1d(f, geometry, settings)
+    elseif Dim == 2
+        return _integral_2d(f, geometry, settings)
+    elseif Dim == 3
+        return _integral_3d(f, geometry, settings)
     end
 end
 
-"""
-    surfaceintegral(f, geometry, algorithm::IntegrationAlgorithm=HAdaptiveCubature)
-
-Numerically integrate a given function `f(::Point)` over a surface `geometry`
-using a particular `integration algorithm`.
-"""
-function surfaceintegral(
+# If algorithm is not HAdaptiveCubature, and FP specified, specialize on number of dimensions
+function integral(
     f::F,
     geometry::G,
-    settings::I=HAdaptiveCubature()
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm}
-    if (dim = paramdim(geometry)) == 2
-        return integral(f, geometry, settings)
-    else
-        error("Performing a surface integral on a geometry with $dim parametric dimensions not supported.")
+    settings::I,
+    FP::Type{T}
+) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm, T<:AbstractFloat}
+    # Run the appropriate integral type
+    Dim = Meshes.paramdim(geometry)
+    if Dim == 1
+        return _integral_1d(f, geometry, settings, FP)
+    elseif Dim == 2
+        return _integral_2d(f, geometry, settings, FP)
+    elseif Dim == 3
+        return _integral_3d(f, geometry, settings, FP)
     end
 end
 
-"""
-    volumeintegral(f, geometry, algorithm::IntegrationAlgorithm=HAdaptiveCubature)
 
-Numerically integrate a given function `f(::Point)` throughout a volumetric
-`geometry` using a particular `integration algorithm`.
+################################################################################
+#                    Generalized (n-Dimensional) Worker Methods
+################################################################################
 
-"""
-function volumeintegral(
-    f::F,
-    geometry::G,
-    settings::I=HAdaptiveCubature()
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm}
-    if (dim = paramdim(geometry)) == 3
-        return integral(f, geometry, settings)
-    else
-        error("Performing a volume integral on a geometry with $dim parametric dimensions not supported.")
-    end
+# General solution for HAdaptiveCubature methods
+function _integral(
+    f,
+    geometry,
+    settings::HAdaptiveCubature,
+    FP::Type{T} = Float64
+) where {T<:AbstractFloat}
+    Dim = Meshes.paramdim(geometry)
+
+    integrand(t) = f(geometry(t...)) * differential(geometry, t)
+
+    # HCubature doesn't support functions that output Unitful Quantity types
+    # Establish the units that are output by f
+    testpoint_parametriccoord = fill(FP(0.5),Dim)
+    integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord))
+    # Create a wrapper that returns only the value component in those units
+    uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
+    # Integrate only the unitless values
+    value = HCubature.hcubature(uintegrand, zeros(FP,Dim), ones(FP,Dim); settings.kwargs...)[1]
+
+    # Reapply units
+    return value .* integrandunits
 end