diff --git a/.typos.toml b/.typos.toml new file mode 100644 index 00000000..0f6b1ff9 --- /dev/null +++ b/.typos.toml @@ -0,0 +1,2 @@ +[default.extend-words] +parametrized = "parametrized" diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index f8f70ee9..b528be0c 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -15,9 +15,9 @@ integrands = ( (name = "Vector", f = p -> fill(norm(to(p)), 3)) ) rules = ( - (name = "GaussLegendre", rule = GaussLegendre(100), N = 100), - (name = "GaussKronrod", rule = GaussKronrod(), N = 100), - (name = "HAdaptiveCubature", rule = HAdaptiveCubature(), N = 500) + (name = "GaussLegendre", rule = GaussLegendre(100)), + (name = "GaussKronrod", rule = GaussKronrod()), + (name = "HAdaptiveCubature", rule = HAdaptiveCubature()) ) geometries = ( (name = "Segment", item = Segment(Point(0, 0, 0), Point(1, 1, 1))), @@ -26,7 +26,8 @@ geometries = ( SUITE["Integrals"] = let s = BenchmarkGroup() for (int, rule, geometry) in Iterators.product(integrands, rules, geometries) - n1, n2, N = geometry.name, "$(int.name) $(rule.name)", rule.N + n1 = geometry.name + n2 = "$(int.name) $(rule.name)" s[n1][n2] = @benchmarkable integral($int.f, $geometry.item, $rule.rule) end s @@ -64,8 +65,7 @@ SUITE["Specializations/Scalar GaussLegendre"] = let s = BenchmarkGroup() 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["Triangle"] = @benchmarkable integral($spec.f, $spec.g.triangle, $spec.rule_gl) - # s["Tetrahedron"] = @benchmarkable integral($spec.f, $spec.g.tetrahedron, $spec.rule) - # TODO re-enable Tetrahedron-GaussLegendre test when supported by main branch + s["Tetrahedron"] = @benchmarkable integral($spec.f, $spec.g.tetrahedron, $spec.rule_gl) s end diff --git a/docs/src/specializations.md b/docs/src/specializations.md index dc873601..ab0da442 100644 --- a/docs/src/specializations.md +++ b/docs/src/specializations.md @@ -11,43 +11,3 @@ There are several notable exceptions to how Meshes.jl defines [parametric functi - `Meshes.Rope` is a composite type that lacks a parametric function, but can be decomposed into `Meshes.Segment`s and then integrated by adding together the individual integrals. - `Meshes.Triangle` has a parametric function that takes coordinates on a 2D barycentric coordinate system. So, for `(::Meshes.Triangle)(t1, t2)`, the coordinates must obey: $t_1, t_2 \in [0,1]$ where $t_1 + t_2 \le 1$. - `Meshes.Tetrahedron` has a parametric function that takes coordinates on a 3D barycentric coordinate system. So, for `(::Meshes.Tetrahedron)(t1, t2)`, the coordinates must obey: $t_1, t_2, t_3 \in [0,1]$ where $t_1 + t_2 + t_3 \le 1$. - -## Triangle - -For a specified `Meshes.Triangle` surface with area $A$, let $u$ and $v$ be Barycentric coordinates that span the surface. -```math -\int_\triangle f(\bar{r}) \, \text{d}A - = \iint_\triangle f\left( \bar{r}(u,v) \right) \, \left( \text{d}u \wedge \text{d}v \right) -``` - -Since the geometric transformation from the originally-arbitrary domain to a Barycentric domain is linear, the magnitude of the surface element $\text{d}u \wedge \text{d}v$ is constant throughout the integration domain. This constant will be equal to twice the magnitude of $A$. -```math -\int_\triangle f(\bar{r}) \, \text{d}A - = 2A \int_0^1 \int_0^{1-v} f\left( \bar{r}(u,v) \right) \, \text{d}u \, \text{d}v -``` - -Since the integral domain is a right-triangle in the Barycentric domain, a nested application of Gauss-Kronrod quadrature rules is capable of computing the result, albeit inefficiently. However, many numerical integration methods that require rectangular bounds can not be directly applied. - -In order to enable integration methods that operate over rectangular bounds, two coordinate system transformations are applied: the first maps from Barycentric coordinates $(u, v)$ to polar coordinates $(r, \phi)$, and the second is a non-linear map from polar coordinates to a new curvilinear basis $(R, \phi)$. - -For the first transformation, let $u = r~\cos\phi$ and $v = r~\sin\phi$ where $\text{d}u~\text{d}v = r~\text{d}r~\text{d}\phi$. The Barycentric triangle's hypotenuse boundary line is described by the function $v(u) = 1 - u$. Substituting in the previous definitions leads to a new boundary line function in polar coordinate space $r(\phi) = 1 / (\sin\phi + \cos\phi)$. -```math -\int_0^1 \int_0^{1-v} f\left( \bar{r}(u,v) \right) \, \text{d}u \, \text{d}v = - \int_0^{\pi/2} \int_0^{1/(\sin\phi+\cos\phi)} f\left( \bar{r}(r,\phi) \right) \, r \, \text{d}r \, \text{d}\phi -``` - -These integral boundaries remain non-rectangular, so an additional transformation will be applied to a curvilinear $(R, \phi)$ space that normalizes all of the hypotenuse boundary line points to $R=1$. To achieve this, a function $R(r,\phi)$ is required such that $R(r_0, \phi) = 1$ where $r_0 = 1 / (\sin\phi + \cos\phi)$ - -To achieve this, let $R(r, \phi) = r~(\sin\phi + \cos\phi)$. Now, substituting some terms leads to -```math -\int_0^{\pi/2} \int_0^{1/(\sin\phi+\cos\phi)} f\left( \bar{r}(r,\phi) \right) \, r \, \text{d}r \, \text{d}\phi - = \int_0^{\pi/2} \int_0^{r_0} f\left( \bar{r}(r,\phi) \right) \, \left(\frac{R}{\sin\phi + \cos\phi}\right) \, \text{d}r \, \text{d}\phi -``` - -Since $\text{d}R/\text{d}r = \sin\phi + \cos\phi$, a change of integral domain leads to -```math -\int_0^{\pi/2} \int_0^{r_0} f\left( \bar{r}(r,\phi) \right) \, \left(\frac{R}{\sin\phi + \cos\phi}\right) \, \text{d}r \, \text{d}\phi - = \int_0^{\pi/2} \int_0^1 f\left( \bar{r}(R,\phi) \right) \, \left(\frac{R}{\left(\sin\phi + \cos\phi\right)^2}\right) \, \text{d}R \, \text{d}\phi -``` - -The second term in this new integrand function serves as a correction factor that corrects for the impact of the non-linear domain transformation. Since all of the integration bounds are now constants, specialized integration methods can be defined for triangles that performs these domain transformations and then solve the new rectangular integration problem using a wider range of solver options. diff --git a/src/MeshIntegrals.jl b/src/MeshIntegrals.jl index 7b270c7b..cbc7260f 100644 --- a/src/MeshIntegrals.jl +++ b/src/MeshIntegrals.jl @@ -10,7 +10,7 @@ import QuadGK import Unitful include("differentiation.jl") -export DifferentiationMethod, Analytical, FiniteDifference, jacobian +export DifferentiationMethod, FiniteDifference, jacobian include("utils.jl") diff --git a/src/differentiation.jl b/src/differentiation.jl index 7ea44db0..95232357 100644 --- a/src/differentiation.jl +++ b/src/differentiation.jl @@ -9,7 +9,7 @@ A category of types used to specify the desired method for calculating derivativ Derivatives are used to form Jacobian matrices when calculating the differential element size throughout the integration region. -See also [`FiniteDifference`](@ref), [`Analytical`](@ref). +See also [`FiniteDifference`](@ref). """ abstract type DifferentiationMethod end @@ -27,22 +27,6 @@ end FiniteDifference{T}() where {T <: AbstractFloat} = FiniteDifference{T}(T(1e-6)) FiniteDifference() = FiniteDifference{Float64}() -""" - Analytical() - -Use to specify use of analytically-derived solutions for calculating derivatives. -These solutions are currently defined only for a subset of geometry types. - -# Supported Geometries: -- `BezierCurve` -- `Line` -- `Plane` -- `Ray` -- `Tetrahedron` -- `Triangle` -""" -struct Analytical <: DifferentiationMethod end - # Future Support: # struct AutoEnzyme <: DifferentiationMethod end # struct AutoZygote <: DifferentiationMethod end diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl index 01b0f24d..a79f5dd3 100644 --- a/src/specializations/BezierCurve.jl +++ b/src/specializations/BezierCurve.jl @@ -39,8 +39,7 @@ function integral( kwargs... ) # Generate a _ParametricGeometry whose parametric function auto-applies the alg kwarg - paramfunction(t) = _parametric(curve, t, alg) - param_curve = _ParametricGeometry(paramfunction, Meshes.paramdim(curve)) + param_curve = _ParametricGeometry(_parametric(curve, alg), Meshes.paramdim(curve)) # Integrate the _ParametricGeometry using the standard methods return _integral(f, param_curve, rule; kwargs...) @@ -50,6 +49,7 @@ end # Parametric ################################################################################ -function _parametric(curve::Meshes.BezierCurve, t, alg::Meshes.BezierEvalMethod) - return curve(t, alg) +# Wrap (::BezierCurve)(t, ::BezierEvalMethod) into f(t) by embedding second argument +function _parametric(curve::Meshes.BezierCurve, alg::Meshes.BezierEvalMethod) + return t -> curve(t, alg) end diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl index 0c7da2c1..e099928c 100644 --- a/src/specializations/Line.jl +++ b/src/specializations/Line.jl @@ -10,81 +10,26 @@ function integral( f, line::Meshes.Line, - rule::GaussLegendre; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Line, diff_method) - - # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] - xs, ws = _gausslegendre(FP, rule.n) - - # Normalize the Line s.t. line(t) is distance t from origin point - line = Meshes.Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 - - # Integrate f along the Line - differential(line, x) = t′(x) * _units(line(0)) - integrand(x) = f(line(t(x))) * differential(line, x) - return sum(w .* integrand(x) for (w, x) in zip(ws, xs)) -end - -function integral( - f, - line::Meshes.Line, - rule::GaussKronrod; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Line, diff_method) - - # Normalize the Line s.t. line(t) is distance t from origin point - line = Meshes.Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) - - # Integrate f along the Line - domainunits = _units(line(0)) - integrand(t) = f(line(t)) * domainunits - return QuadGK.quadgk(integrand, FP(-Inf), FP(Inf); rule.kwargs...)[1] -end - -function integral( - f, - line::Meshes.Line, - rule::HAdaptiveCubature; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Line, diff_method) - - # Normalize the Line s.t. line(t) is distance t from origin point - line = Meshes.Line(line.a, line.a + Meshes.unormalize(line.b - line.a)) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 - - # Integrate f along the Line - differential(line, x) = t′(x) * _units(line(0)) - integrand(xs) = f(line(t(xs[1]))) * differential(line, xs[1]) - - # HCubature doesn't support functions that output Unitful Quantity types - # Establish the units that are output by f - testpoint_parametriccoord = (FP(0.5),) - 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, (-one(FP),), (one(FP),); rule.kwargs...)[1] - - # Reapply units - return value .* integrandunits + rule::IntegrationRule; + kwargs... +) + # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1] + param_line = _ParametricGeometry(_parametric(line), Meshes.paramdim(line)) + + # Integrate the _ParametricGeometry using the standard methods + return _integral(f, param_line, rule; kwargs...) end ################################################################################ -# jacobian +# Parametric ################################################################################ -_has_analytical(::Type{T}) where {T <: Meshes.Line} = true +# Map argument domain from [0, 1] to (-∞, ∞) for (::Line)(t) +function _parametric(line::Meshes.Line) + # [-1, 1] ↦ (-∞, ∞) + f1(t) = t / (1 - t^2) + # [0, 1] ↦ [-1, 1] + f2(t) = 2t - 1 + # Compose the two transforms + return t -> line((f1 ∘ f2)(t)) +end diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl index ee2706cb..df7c4e2d 100644 --- a/src/specializations/Plane.jl +++ b/src/specializations/Plane.jl @@ -11,96 +11,27 @@ function integral( f, plane::Meshes.Plane, - rule::GaussLegendre; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Plane, diff_method) - - # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]² - xs, ws = _gausslegendre(FP, rule.n) - wws = Iterators.product(ws, ws) - xxs = Iterators.product(xs, xs) - - # Normalize the Plane's orthogonal vectors - uu = Meshes.unormalize(plane.u) - uv = Meshes.unormalize(plane.v) - uplane = Meshes.Plane(plane.p, uu, uv) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 - - # Integrate f over the Plane - domainunits = _units(uplane(0, 0)) - function integrand(((wi, wj), (xi, xj))) - wi * wj * f(uplane(t(xi), t(xj))) * t′(xi) * t′(xj) * domainunits^2 - end - return sum(integrand, zip(wws, xxs)) + rule::IntegrationRule; + kwargs... +) + # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1]² + param_plane = _ParametricGeometry(_parametric(plane), Meshes.paramdim(plane)) + + # Integrate the _ParametricGeometry using the standard methods + return _integral(f, param_plane, rule; kwargs...) end -function integral( - f, - plane::Meshes.Plane, - rule::GaussKronrod; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Plane, diff_method) - - # Normalize the Plane's orthogonal vectors - uu = Meshes.unormalize(plane.u) - uv = Meshes.unormalize(plane.v) - uplane = Meshes.Plane(plane.p, uu, uv) - - # Integrate f over the Plane - domainunits = _units(uplane(0, 0))^2 - integrand(u, v) = f(uplane(u, v)) * domainunits - inner∫(v) = QuadGK.quadgk(u -> integrand(u, v), FP(-Inf), FP(Inf); rule.kwargs...)[1] - return QuadGK.quadgk(inner∫, FP(-Inf), FP(Inf); rule.kwargs...)[1] +############################################################################################ +# Parametric +############################################################################################ + +# Map argument domain from [0, 1]² to (-∞, ∞)² for (::Plane)(t1, t2) +function _parametric(plane::Meshes.Plane) + # [-1, 1] ↦ (-∞, ∞) + f1(t) = t / (1 - t^2) + # [0, 1] ↦ [-1, 1] + f2(t) = 2t - 1 + # Compose the two transforms + f = f1 ∘ f2 + return (t1, t2) -> plane(f(t1), f(t2)) end - -function integral( - f, - plane::Meshes.Plane, - rule::HAdaptiveCubature; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Plane, diff_method) - - # Normalize the Plane's orthogonal vectors - uu = Meshes.unormalize(plane.u) - uv = Meshes.unormalize(plane.v) - uplane = Meshes.Plane(plane.p, uu, uv) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ (-∞,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 - - # Integrate f over the Plane - domainunits = _units(uplane(0, 0)) - function integrand(xs) - f(uplane(t(xs[1]), t(xs[2]))) * t′(xs[1]) * t′(xs[2]) * domainunits^2 - end - - # HCubature doesn't support functions that output Unitful Quantity types - # Establish the units that are output by f - testpoint_parametriccoord = zeros(FP, 2) - 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 - a = .-_ones(FP, 2) - b = _ones(FP, 2) - value = HCubature.hcubature(uintegrand, a, b; rule.kwargs...)[1] - - # Reapply units - return value .* integrandunits -end - -################################################################################ -# jacobian -################################################################################ - -_has_analytical(::Type{T}) where {T <: Meshes.Plane} = true diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl index eede812e..8c6e9e98 100644 --- a/src/specializations/Ray.jl +++ b/src/specializations/Ray.jl @@ -11,85 +11,22 @@ function integral( f, ray::Meshes.Ray, - rule::GaussLegendre; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Ray, diff_method) - - # Compute Gauss-Legendre nodes/weights for x in interval [-1,1] - xs, ws = _gausslegendre(FP, rule.n) - - # Normalize the Ray s.t. ray(t) is distance t from origin point - ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v)) - - # Domain transformation: x ∈ [-1,1] ↦ t ∈ [0,∞) - t₁(x) = (1 // 2) * x + (1 // 2) - t₁′(x) = (1 // 2) - t₂(x) = x / (1 - x^2) - t₂′(x) = (1 + x^2) / (1 - x^2)^2 - t = t₂ ∘ t₁ - t′(x) = t₂′(t₁(x)) * t₁′(x) - - # Integrate f along the Ray - domainunits = _units(ray(0)) - integrand(x) = f(ray(t(x))) * t′(x) * domainunits - return sum(w .* integrand(x) for (w, x) in zip(ws, xs)) -end - -function integral( - f, - ray::Meshes.Ray, - rule::GaussKronrod; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Ray, diff_method) - - # Normalize the Ray s.t. ray(t) is distance t from origin point - ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v)) - - # Integrate f along the Ray - domainunits = _units(ray(0)) - integrand(t) = f(ray(t)) * domainunits - return QuadGK.quadgk(integrand, zero(FP), FP(Inf); rule.kwargs...)[1] -end - -function integral( - f, - ray::Meshes.Ray, - rule::HAdaptiveCubature; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Ray, diff_method) - - # Normalize the Ray s.t. ray(t) is distance t from origin point - ray = Meshes.Ray(ray.p, Meshes.unormalize(ray.v)) - - # Domain transformation: x ∈ [0,1] ↦ t ∈ [0,∞) - t(x) = x / (1 - x^2) - t′(x) = (1 + x^2) / (1 - x^2)^2 - - # Integrate f along the Ray - domainunits = _units(ray(0)) - integrand(xs) = f(ray(t(xs[1]))) * t′(xs[1]) * domainunits - - # HCubature doesn't support functions that output Unitful Quantity types - # Establish the units that are output by f - testpoint_parametriccoord = _zeros(FP, 1) - 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, 1), _ones(FP, 1); rule.kwargs...)[1] - - # Reapply units - return value .* integrandunits + rule::IntegrationRule; + kwargs... +) + # Generate a _ParametricGeometry whose parametric function spans the domain [0, 1] + param_ray = _ParametricGeometry(_parametric(ray), Meshes.paramdim(ray)) + + # Integrate the _ParametricGeometry using the standard methods + return _integral(f, param_ray, rule; kwargs...) end ################################################################################ -# jacobian +# Parametric ################################################################################ -_has_analytical(::Type{T}) where {T <: Meshes.Ray} = true +# Map argument domain from [0, 1] to [0, ∞) for (::Ray)(t) +function _parametric(ray::Meshes.Ray) + f(t) = t / (1 - t^2) + return t -> ray(f(t)) +end diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl index 2f6fdece..f074a2fd 100644 --- a/src/specializations/Tetrahedron.jl +++ b/src/specializations/Tetrahedron.jl @@ -15,9 +15,8 @@ function integral( rule::IntegrationRule; kwargs... ) - # Generate a _ParametricGeometry whose parametric function domain spans [0,1]^3 - paramfunction(t1, t2, t3) = _parametric(tetrahedron, t1, t2, t3) - tetra = _ParametricGeometry(paramfunction, Meshes.paramdim(tetrahedron)) + # Generate a _ParametricGeometry whose parametric function domain spans [0,1]³ + tetra = _ParametricGeometry(_parametric(tetrahedron), Meshes.paramdim(tetrahedron)) # Integrate the _ParametricGeometry using the standard methods return _integral(f, tetra, rule; kwargs...) @@ -27,17 +26,21 @@ end # Parametric ################################################################################ -function _parametric(tetrahedron::Meshes.Tetrahedron, t1, t2, t3) - if (t3 < 0) || (t3 > 1) - msg = "tetrahedron(t1, t2, t3) is not defined for t3 outside [0, 1]." - throw(DomainError(t3, msg)) - end +# Map argument domain from [0, 1]³ to Barycentric domain for (::Tetrahedron)(t1, t2, t3) +function _parametric(tetrahedron::Meshes.Tetrahedron) + function f(t1, t2, t3) + if any(Iterators.map(n -> (n < 0) || (n > 1), (t1, t2, t3))) + msg = "tetrahedron(t1, t2, t3) is not defined for (t1, t2, t3) outside [0, 1]³." + throw(DomainError((t1, t2, t3), msg)) + end - # Take a triangular cross-section at t3 - a = tetrahedron(t3, 0, 0) - b = tetrahedron(0, t3, 0) - c = tetrahedron(0, 0, t3) - cross_section = Meshes.Triangle(a, b, c) + # Take a triangular cross-section at t3 + a = tetrahedron(t3, 0, 0) + b = tetrahedron(0, t3, 0) + c = tetrahedron(0, 0, t3) + cross_section = _parametric(Meshes.Triangle(a, b, c)) - return _parametric(cross_section, t1, t2) + return cross_section(t1, t2) + end + return f end diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl index 2c70f400..26f3ea15 100644 --- a/src/specializations/Triangle.jl +++ b/src/specializations/Triangle.jl @@ -12,115 +12,30 @@ function integral( f, triangle::Meshes.Triangle, - rule::GaussLegendre; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Triangle, diff_method) - - # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]^2 - xs, ws = _gausslegendre(FP, rule.n) - wws = Iterators.product(ws, ws) - xxs = Iterators.product(xs, xs) - - # Domain transformations: - # xᵢ [-1,1] ↦ R [0,1] - # xⱼ [-1,1] ↦ φ [0,π/2] - uR(xᵢ) = (1 // 2) * (xᵢ + 1) - uφ(xⱼ) = (T(π) / 4) * (xⱼ + 1) - - # Integrate the Barycentric triangle by transforming it into polar coordinates - # with a modified radius - # R = r ( sinφ + cosφ ) - # s.t. integration bounds become rectangular - # R ∈ [0, 1] and φ ∈ [0, π/2] - function integrand(((wᵢ, wⱼ), (xᵢ, xⱼ))) - R = uR(xᵢ) - φ = uφ(xⱼ) - a, b = sincos(φ) - u = R * (1 - a / (a + b)) - v = R * (1 - b / (a + b)) - return wᵢ * wⱼ * f(triangle(u, v)) * R / (a + b)^2 - end - - # Calculate 2D Gauss-Legendre integral over modified-polar-Barycentric coordinates - # Apply a linear domain-correction factor - return FP(π / 4) * Meshes.area(triangle) .* sum(integrand, zip(wws, xxs)) -end - -function integral( - f, - triangle::Meshes.Triangle, - rule::GaussKronrod; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Triangle, diff_method) - - # Integrate the Barycentric triangle in (u,v)-space: (0,0), (0,1), (1,0) - # i.e. \int_{0}^{1} \int_{0}^{1-u} f(u,v) dv du - ∫u(u) = QuadGK.quadgk(v -> f(triangle(u, v)), zero(FP), FP(1 - u); rule.kwargs...)[1] - ∫ = QuadGK.quadgk(∫u, zero(FP), one(FP); rule.kwargs...)[1] - - # Apply a linear domain-correction factor 0.5 ↦ area(triangle) - return 2 * Meshes.area(triangle) .* ∫ -end - -function integral( - f, - triangle::Meshes.Triangle, - rule::HAdaptiveCubature; - diff_method::DM = Analytical(), - FP::Type{T} = Float64 -) where {DM <: DifferentiationMethod, T <: AbstractFloat} - _guarantee_analytical(Meshes.Triangle, diff_method) - - # Integrate the Barycentric triangle by transforming it into polar coordinates - # with a modified radius - # R = r ( sinφ + cosφ ) - # s.t. integration bounds become rectangular - # R ∈ [0, 1] and φ ∈ [0, π/2] - function integrand(Rφ) - R, φ = Rφ - a, b = sincos(φ) - u = R * (1 - a / (a + b)) - v = R * (1 - b / (a + b)) - return f(triangle(u, v)) * R / (a + b)^2 - end - - # HCubature doesn't support functions that output Unitful Quantity types - # Establish the units that are output by f - testpoint_parametriccoord = _zeros(T, 2) - integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord)) - # Create a wrapper that returns only the value component in those units - uintegrand(ts) = Unitful.ustrip.(integrandunits, integrand(ts)) - # Integrate only the unitless values - ∫ = HCubature.hcubature(uintegrand, _zeros(FP, 2), (T(1), T(π / 2)), rule.kwargs...)[1] - - # Apply a linear domain-correction factor 0.5 ↦ area(triangle) - return 2 * Meshes.area(triangle) * integrandunits .* ∫ + rule::IntegrationRule; + kwargs... +) + # Generate a _ParametricGeometry whose parametric function domain spans [0,1]² + param_triangle = _ParametricGeometry(_parametric(triangle), Meshes.paramdim(triangle)) + + # Integrate the _ParametricGeometry using the standard methods + return _integral(f, param_triangle, rule; kwargs...) end -################################################################################ -# jacobian -################################################################################ - -_has_analytical(::Type{T}) where {T <: Meshes.Triangle} = true - ################################################################################ # Parametric ################################################################################ -function _parametric(triangle::Meshes.Triangle, t1, t2) - if (t1 < 0 || t1 > 1) || (t2 < 0 || t2 > 1) - msg = "triangle(t1, t2) is not defined for (t1, t2) outside [0, 1]^2." - throw(DomainError((t1, t2), msg)) - end - - # Form a line segment between sides - a = triangle(0, t2) - b = triangle(t2, 0) - segment = Meshes.Segment(a, b) +# Map argument domain from [0, 1]² to Barycentric domain for (::Triangle)(t1, t2) +function _parametric(triangle::Meshes.Triangle) + function f(t1, t2) + if any(Iterators.map(n -> (n < 0) || (n > 1), (t1, t2))) + msg = "triangle(t1, t2) is not defined for (t1, t2) outside [0, 1]²." + throw(DomainError((t1, t2), msg)) + end - return segment(t1) + t1t2 = t1 * t2 + return triangle(t1t2, t2 - t1t2) + end + return f end diff --git a/src/specializations/_ParametricGeometry.jl b/src/specializations/_ParametricGeometry.jl index 16904d27..bdf0f6c9 100644 --- a/src/specializations/_ParametricGeometry.jl +++ b/src/specializations/_ParametricGeometry.jl @@ -55,7 +55,7 @@ end Meshes.paramdim(::_ParametricGeometry{M, C, F, Dim}) where {M, C, F, Dim} = Dim """ - _parametric(geometry::G, ts...) where {G <: Meshes.Geometry} + _parametric(geometry::G) where {G <: Meshes.Geometry} -> Function Used in MeshIntegrals.jl for defining parametric functions that transform non-standard geometries into a form that can be integrated over the standard rectangular [0,1]^n limits. diff --git a/src/utils.jl b/src/utils.jl index e69a38ef..1e5292f8 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -26,25 +26,11 @@ _ones(N::Int) = _ones(Float64, N) # DifferentiationMethod ################################################################################ -# Throw an ArgumentError if Analytical() jacobian not defined for this type -function _guarantee_analytical( - ::Type{G}, - ::DifferentiationMethod -) where {G <: Geometry} - throw(ArgumentError("$G geometries require kwarg diff_method = Analytical()")) -end - -_guarantee_analytical(::Type{G}, ::Analytical) where {G <: Geometry} = nothing - -# Return whether a geometry type has jacobian methods defined -_has_analytical(::Type{G}) where {G <: Geometry} = false -_has_analytical(g::G) where {G <: Geometry} = _has_analytical(G) - # Return the default DifferentiationMethod instance for a particular geometry type function _default_method( g::Type{G} ) where {G <: Geometry} - return _has_analytical(G) ? Analytical() : FiniteDifference() + return FiniteDifference() end # Return the default DifferentiationMethod instance for a particular geometry instance diff --git a/test/utils.jl b/test/utils.jl index 69a8fde6..545b2540 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -20,33 +20,16 @@ end @testitem "DifferentiationMethod" setup=[Setup] begin - using MeshIntegrals: _has_analytical, _default_method, _guarantee_analytical + using MeshIntegrals: _default_method - # _has_analytical of instances - triangle = Triangle(Point(0, 0, 0), Point(0, 1, 0), Point(1, 0, 0)) - @test _has_analytical(triangle) == true + # Test geometries sphere = Sphere(Point(0, 0, 0), 1.0) - @test _has_analytical(sphere) == false + triangle = Triangle(Point(0, 0, 0), Point(0, 1, 0), Point(1, 0, 0)) # _default_method - @test _default_method(Meshes.Triangle) isa Analytical - @test _default_method(triangle) isa Analytical @test _default_method(Meshes.Sphere) isa FiniteDifference @test _default_method(sphere) isa FiniteDifference - # _has_analytical of types - @test _has_analytical(Meshes.BezierCurve) == false - @test _has_analytical(Meshes.Line) == true - @test _has_analytical(Meshes.Plane) == true - @test _has_analytical(Meshes.Ray) == true - @test _has_analytical(Meshes.Sphere) == false - @test _has_analytical(Meshes.Tetrahedron) == false - @test _has_analytical(Meshes.Triangle) == true - - # _guarantee_analytical - @test _guarantee_analytical(Meshes.Line, Analytical()) === nothing - @test_throws "Analytical" _guarantee_analytical(Meshes.Line, FiniteDifference()) - # FiniteDifference @test FiniteDifference().ε ≈ 1e-6 end @@ -58,18 +41,14 @@ end pt_w = Point(-7, 0, 0) pt_e = Point(8, 0, 0) ẑ = Vec(0, 0, 1) - triangle = Triangle(pt_n, pt_w, pt_e) - tetrahedron = Tetrahedron(pt_n, pt_w, pt_e, pt_n + ẑ) + triangle = _parametric(Triangle(pt_n, pt_w, pt_e)) + tetrahedron = _parametric(Tetrahedron(pt_n, pt_w, pt_e, pt_n + ẑ)) - # Ensure error is thrown for an out-of-bounds coordinate - for ts in [(-0.5, 0.5), (0.5, -0.5), (1.5, 0.5), (0.5, 1.5)] - @test_throws "not defined" _parametric(triangle, ts...) - # Tetrahedron forwards t1 and t2 to _parametric(::Triangle, ts...) - @test_throws "not defined" _parametric(tetrahedron, ts..., 0) + # Ensure error is thrown for out-of-bounds coordinate + for ts in [(-1, 0), (0, -1), (2, 0), (0, 2)] + @test_throws "not defined" triangle(ts...) end - - # Ensue error is thrown for an out-of-bounds third coordinate - for t3 in [-0.5, 1.5] - @test_throws "not defined" _parametric(tetrahedron, 0, 0, t3) + for ts in [(-1, 0, 0), (0, -1, 0), (0, 0, -1), (2, 0, 0), (0, 2, 0), (0, 0, 2)] + @test_throws "not defined" tetrahedron(ts...) end end