Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement more _ParametricGeometry transforms #139

Merged
merged 34 commits into from
Dec 2, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
34 commits
Select commit Hold shift + click to select a range
93b1971
Implement _ParametricGeometry for Line
mikeingold Nov 30, 2024
8f3502b
Implement _ParametricGeometry for Plane
mikeingold Nov 30, 2024
201ae62
Implement _ParametricGeometry for Ray
mikeingold Nov 30, 2024
cdac6ec
Update and reorganize analytical tests
mikeingold Nov 30, 2024
f4518c8
Add comments
mikeingold Nov 30, 2024
3ce948c
Use paramdim for clarity
mikeingold Nov 30, 2024
7ce264e
Test _parametric -> anonymous function
mikeingold Nov 30, 2024
8e9c6c8
Update _parametric API to return an anonymous function
mikeingold Nov 30, 2024
77f39f4
Remove obsolete N, enable Tetrahedron test
mikeingold Nov 30, 2024
0b3739a
Remove hanging comma
mikeingold Nov 30, 2024
2e06577
Update field name
mikeingold Nov 30, 2024
5f404f3
Add analytical derivative in comment
mikeingold Nov 30, 2024
dce7c3a
Update _parametric for Triangle and Tetrahedron
mikeingold Dec 1, 2024
63e85fe
Update tests
mikeingold Dec 1, 2024
a2b4953
Merge branch 'main' into line
JoshuaLampert Dec 1, 2024
b5c6886
put parametrized to .typos.toml
JoshuaLampert Dec 1, 2024
4217d11
Implement _ParametricGeometry for Triangle
mikeingold Dec 2, 2024
8b8ff28
Temp disable some Analytical tests
mikeingold Dec 2, 2024
1cac830
Temp disable more tests
mikeingold Dec 2, 2024
044a20d
Test disabling bound check
mikeingold Dec 2, 2024
15710a6
Restore bounds check, try new formulation
mikeingold Dec 2, 2024
7d3c83e
Remove old commented-out code
mikeingold Dec 2, 2024
dbdffed
Add comments
mikeingold Dec 2, 2024
bc4e591
Remove obsolete tests
mikeingold Dec 2, 2024
3178e8f
Remove obsolete utils and update default diff_method
mikeingold Dec 2, 2024
9616901
Remove obsolete imports
mikeingold Dec 2, 2024
c21b433
Remove obsolete Analytical
mikeingold Dec 2, 2024
8a8f892
Remove obsolete export
mikeingold Dec 2, 2024
6743a04
Formatting
mikeingold Dec 2, 2024
de2c784
Remove obsolete derivation
mikeingold Dec 2, 2024
4b3ad9d
Update src/specializations/Tetrahedron.jl
mikeingold Dec 2, 2024
c3637cd
Remove docstring ref to Analytical
mikeingold Dec 2, 2024
5b313f5
Update src/specializations/Triangle.jl
mikeingold Dec 2, 2024
29346fa
Update src/specializations/Triangle.jl
mikeingold Dec 2, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions .typos.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
[default.extend-words]
parametrized = "parametrized"
12 changes: 6 additions & 6 deletions benchmark/benchmarks.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))),
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand Down
40 changes: 0 additions & 40 deletions docs/src/specializations.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.
2 changes: 1 addition & 1 deletion src/MeshIntegrals.jl
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ import QuadGK
import Unitful

include("differentiation.jl")
export DifferentiationMethod, Analytical, FiniteDifference, jacobian
export DifferentiationMethod, FiniteDifference, jacobian

include("utils.jl")

Expand Down
18 changes: 1 addition & 17 deletions src/differentiation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand Down
8 changes: 4 additions & 4 deletions src/specializations/BezierCurve.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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...)
Expand All @@ -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
91 changes: 18 additions & 73 deletions src/specializations/Line.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
111 changes: 21 additions & 90 deletions src/specializations/Plane.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Loading
Loading