Skip to content

Commit 49aaf57

Browse files
mikeingoldgithub-actions[bot]JoshuaLampert
authored
Use a new DifferentiationMethod type for specifying derivative methods (#122)
* Add JacobianMethod types * Change JacobianMethod to DifferentiationMethod * Add a third argument to jacobian for DiffMethod for dispatch * Remove dead line * Remove obsolete analytic jacobian for BezierCurve * Add default step size if unspecified * Convert differential to three-argument with default FiniteDiff method * Bigfix * Add two-arg Jacobian pointing to FD * Apply formatting suggestion Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Apply suggestion Co-authored-by: Joshua Lampert <[email protected]> * Add a dt argument * Update docstring API * Add dt kwargs * Add missing params * Add dt kwargs (with TODOs) * Consolidate kwarg naming * Update src/differentiation.jl Co-authored-by: Joshua Lampert <[email protected]> * Add exports * Tweak suggested naming [skip CI] * Update Line methods to use an internal differential method * Bump version to v0.16.0-DEV * Add block for DifferentiationMethods * Drop abstract type * Add diff_method arguments * Change dt to diff_method * Prepare to revert jacobian of BezierCurve (pending new benchmarks) * Bump CI * Add BézierCurve to benchmarks (#124) (#126) * Add Bézier entry * Missing comma, added newline * Make item implicitly Unitful * Make test item explicit to avoid reference errors * Formatting * Reformat range * Apply format suggestion --------- Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Reduce scald * Activate specialized Jacobian * Define an Analytic method and update jacobian of Bezier to it * Typo [skip CI] * Typo * Update src/specializations/BezierCurve.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Implement has_analytic * Bugfix * Enable analytical for BezierCurve [skip CI] * Implement default_method * Fix arg definition * Update arg type * Add methods on instances and types * Add tests for DifferentiationMethods * Use default_method in _integral methods * Use default_method in Bezier methods, improve docstring * Improve docstrings * Explicit namespace for test of non-exported function * Fix typo * Fix copy/paste error * Add tests for default_method * Type vs instance * Export Analytical * Typo * Typo * Parameterize Bezier type * Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Redefine Line, Plane, and Ray specializations as Analytical * Redefine Triangle and Tetrahedron specializations as Analytical * Try new solution (cant dispatch on kwarg types) * Add jacobian code sections * Add guarantees * Move new utilities to utils and rename with prefix underscore * Typo * Update tests * Update name * Update function name * Update comment * Update name * Add better docstrings * Bugfix - error not actually thrown * Split DifferentiationMethod into separate testitem * Tweak include order * Formatting changes * Remove redundant docstrings * Add docstrings for Ring and Rope, update for Bezier * Apply suggestions from code review Co-authored-by: Joshua Lampert <[email protected]> * (WIP) Add a new kwarg section under Usage * Clarify Usage section * Formatting * Add a WIP differential forms explainer * Rename page and add some structure * Add new page to directory * Update index page with README contents * Update first two sections * Finished version 0 of How It Works * Move How it Works into new Dev Notes section * Properly enforce FiniteDiff FP rules * Bugfix - restore G param * Bugfix - tuple to ntuple * Remove Unicode from math blocks * New default constructors for FiniteDifference * Implement two-arg _default_method with FP type * Reorder kwargs and switch FP to T * Revert to single-arg _default_method * Change jacobian parameterization * Remove ununsed type param * Remove unneeded type param * Remove remnant two-arg default call * Clarify _guarantee_analytical arg types * Test externalizing partial derivative function * Re-internalize partial derivative function * Remove errant character * Apply format suggestions [skip ci] Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Apply suggestions from code review [skip ci] Co-authored-by: Joshua Lampert <[email protected]> * Apply suggestions from code review Co-authored-by: Joshua Lampert <[email protected]> * Try new parameterization * Remove obsolete comment * Apply suggestions from code review Co-authored-by: Joshua Lampert <[email protected]> * Consolidate information on specializations * Apply suggestion * Apply suggestion Co-authored-by: Joshua Lampert <[email protected]> --------- Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: Joshua Lampert <[email protected]>
1 parent c5d559c commit 49aaf57

20 files changed

+509
-147
lines changed

Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "MeshIntegrals"
22
uuid = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47"
3-
version = "0.15.2"
3+
version = "0.16.0-DEV"
44

55
[deps]
66
CliffordNumbers = "3998ac73-6bd4-4031-8035-f167dd3ed523"

README.md

Lines changed: 10 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -20,27 +20,31 @@ These solvers have support for integrand functions that produce scalars, vectors
2020

2121
## Usage
2222

23+
### Basic
24+
2325
```julia
2426
integral(f, geometry)
2527
```
2628
Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry`. A default integration method will be automatically selected according to the geometry: `GaussKronrod()` for 1D, and `HAdaptiveCubature()` for all others.
2729

2830
```julia
29-
integral(f, geometry, rule, FP=Float64)
31+
integral(f, geometry, rule)
3032
```
3133
Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry` using the specified integration rule, e.g. `GaussKronrod()`.
3234

33-
Optionally, a fourth argument can be provided to specify the floating point precision level desired. This setting can be manipulated if your integrand function produces outputs with alternate floating point precision (e.g. `Float16`, `BigFloat`, etc) AND you'd prefer to avoid implicit type promotions.
35+
Additionally, several optional keyword arguments are defined in [the API](https://juliageometry.github.io/MeshIntegrals.jl/stable/api/) to provide additional control over the integration mechanics.
36+
37+
### Aliases
3438

3539
```julia
3640
lineintegral(f, geometry)
3741
surfaceintegral(f, geometry)
3842
volumeintegral(f, geometry)
3943
```
40-
Alias functions are provided for convenience. These are simply wrappers for `integral` that first validate that the provided `geometry` has the expected number of parametric/manifold dimensions. Like with `integral` in the examples above, the `rule` can also be specified as a third-argument.
41-
- `lineintegral` for curve-like geometries or polytopes (e.g. `Segment`, `Ray`, `BezierCurve`, `Rope`, etc)
42-
- `surfaceintegral` for surfaces (e.g. `Disk`, `Sphere`, `CylinderSurface`, etc)
43-
- `volumeintegral` for (3D) volumes (e.g. `Ball`, `Cone`, `Torus`, etc)
44+
Alias functions are provided for convenience. These are simply wrappers for `integral` that also validate that the provided `geometry` has the expected number of parametric dimensions. Like with `integral`, a `rule` can also optionally be specified as a third argument.
45+
- `lineintegral` is used for curve-like geometries or polytopes (e.g. `Segment`, `Ray`, `BezierCurve`, `Rope`, etc)
46+
- `surfaceintegral` is used for surfaces (e.g. `Disk`, `Sphere`, `CylinderSurface`, etc)
47+
- `volumeintegral` is used for (3D) volumes (e.g. `Ball`, `Cone`, `Torus`, etc)
4448

4549
# Demo
4650

benchmark/benchmarks.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,7 @@ geometries = (
2828
SUITE["Integrals"] = let s = BenchmarkGroup()
2929
for (int, rule, geometry) in Iterators.product(integrands, rules, geometries)
3030
n1, n2, N = geometry.name, "$(int.name) $(rule.name)", rule.N
31-
s[n1][n2] = @benchmarkable integral($int.f, $geometry.item, $rule.rule) evals=N
31+
s[n1][n2] = @benchmarkable integral($int.f, $geometry.item, $rule.rule)
3232
end
3333
s
3434
end
@@ -41,8 +41,8 @@ sphere = Sphere(Point(0, 0, 0), 1.0)
4141
differential = MeshIntegrals.differential
4242

4343
SUITE["Differentials"] = let s = BenchmarkGroup()
44-
s["Jacobian"] = @benchmarkable jacobian($sphere, $(0.5, 0.5)) evals=1000
45-
s["Differential"] = @benchmarkable differential($sphere, $(0.5, 0.5)) evals=1000
44+
s["Jacobian"] = @benchmarkable jacobian($sphere, $(0.5, 0.5)) evals=10
45+
s["Differential"] = @benchmarkable differential($sphere, $(0.5, 0.5)) evals=10
4646
s
4747
end
4848

docs/make.jl

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,7 @@ using MeshIntegrals
44
# Dynamically set all files in subdirectories of the source directory to include all files in these subdirectories.
55
# This way they don't need to be listed explicitly.
66
SPECIALIZATIONS_FILES = joinpath.(Ref("specializations"),
7-
readdir(joinpath(dirname(@__DIR__), "src",
8-
"specializations")))
7+
readdir(joinpath(dirname(@__DIR__), "src", "specializations")))
98

109
makedocs(
1110
sitename = "MeshIntegrals.jl",
@@ -15,8 +14,9 @@ makedocs(
1514
"Support Matrix" => "supportmatrix.md",
1615
"Example Usage" => "usage.md"
1716
],
18-
"Derivations" => [
19-
"Integrating a Triangle" => "triangle.md"
17+
"Developer Notes" => [
18+
"How it Works" => "how_it_works.md",
19+
"Specializations" => "specializations.md"
2020
],
2121
"Public API" => "api.md"
2222
]

docs/src/how_it_works.md

Lines changed: 76 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,76 @@
1+
# How it Works (By Example)
2+
3+
## Example Problem
4+
5+
Let $f$ be a function of position $\bar{r}$ in some space.
6+
```julia
7+
function f(r̄::Meshes.Point)
8+
x, y, z = to(r̄)
9+
...
10+
end
11+
```
12+
13+
Let the integration domain be the space (a ball) enclosed by a sphere centered on the origin with a radius of 5 meters.
14+
```julia
15+
center = Meshes.Point(0u"m", 0u"m", 0u"m")
16+
radius = 5.0u"m"
17+
ball = Meshes.Ball(center, radius)
18+
```
19+
20+
This integral is often expressed abstractly as simply the following, where the triple integral signs and $\text{d}V$ indicate that the integration domain is some three-dimensional volume.
21+
```math
22+
\iiint f(\bar{r}) ~ \text{d}V
23+
```
24+
25+
Integrals like this are often solved manually by selecting an appropriate coordinate system and limits that neatly represent the integration domain, e.g.
26+
```math
27+
\int_0^{\pi} \int_0^{2\pi} \int_0^{5} f(\bar{r}) ~ \text{d}\rho~\text{d}\theta~\text{d}\phi
28+
```
29+
30+
This works great for simple geometries, but requires integration code that is geometry-specific. This package leverages parametric functions defined in Meshes.jl and differential forms to define integral methods that are general solutions for all geometries.
31+
32+
## [Parametric Functions](@id how-parametric)
33+
34+
Every supported `Meshes.Geometry` type is defined as having a parametric function that maps from a local parametric coordinate system to every point on the geometry. Curve-like geometries will have a single parametric dimension, surfaces will have two dimensions, and volumes will have three dimensions; this can be checked for a particular geometry via `Meshes.paramdim(geometry)`.
35+
36+
For consistency across geometry types, with [some notable exceptions](@ref specializations), these parametric functions are defined to take coordinates inside a normalized range $[0,1]$. In the example case of `ball`, Meshes.jl defines a parametric function mapped in normalized spherical coordinates $(t_\rho, ~t_\theta, ~t_\phi)$. We find, then:
37+
```julia
38+
Meshes.paramdim(ball) == 3 # a volume
39+
40+
ball(tρ, tθ, tφ) # for args in range [0, 1], maps to a corresponding Meshes.Point
41+
42+
ball(0, tθ, tφ) == center
43+
```
44+
45+
In effect, we can now use the geometry itself as a function that maps from three normalized ($0 \le t \le 1$) arguments to every point on the geometry. For the sake of generalization, let this parametric function be called $g$.
46+
```math
47+
\text{g}: (t_1,~t_2,~t_3) ~\mapsto~ \text{Point}\big[ x, ~y, ~z \big]
48+
```
49+
50+
## Differential Forms
51+
52+
Using differential forms, the general solution for integrating a geometry with three parametric dimensions ($t_1$, $t_2$, and $t_3$) is
53+
```math
54+
\iiint f(r̄) ~ \text{d}V = \iiint f(\bar{r}) ~ \bar{\text{d}t_1} \wedge \bar{\text{d}t_2} \wedge \bar{\text{d}t_3}
55+
```
56+
57+
This resultant differential (volume) element is formed at each point in the integration domain by taking [the Jacobian](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant) of the parametric function.
58+
```math
59+
\mathbf{J}_f = \begin{bmatrix} \bar{\text{d}t_1} & \bar{\text{d}t_2} & \bar{\text{d}t_3} \end{bmatrix}
60+
```
61+
where
62+
```math
63+
\bar{\text{d}t_n} = \frac{\partial}{\partial t_n} ~ \text{g}(t_1,~t_2,~t_3)
64+
```
65+
66+
Each of these partial derivatives is a vector representing the direction that changing each parametric function argument will move the resultant point. The differential element ($E$) size is then calculated using geometric algebra as the magnitude of the exterior product ($\wedge$) of these three vectors.
67+
```math
68+
E(t_1,~t_2,~t_3) = \left\| \bar{\text{d}t_1} \wedge \bar{\text{d}t_2} \wedge \bar{\text{d}t_3} \right\|
69+
```
70+
71+
Finally, we use the parametric function itself, $g$, as a map to all points $\bar{r}$ in the integration domain. Since `Meshes.Geometry` parametric functions all operate on normalized domains, we can now solve any volume integral as simply
72+
```math
73+
\int_0^1 \int_0^1 \int_0^1 f\Big(\text{g}\big(t_1,~t_2,~t_3\big)\Big) ~ E(t_1,~t_2,~t_3) ~ \text{d}t_1 ~ \text{d}t_2 ~ \text{d}t_3
74+
```
75+
76+
This form of integral can be trivially generalized to support $n$-dimensional geometries in a form that enables the use of a wide range of numerical integration libraries.

docs/src/index.md

Lines changed: 35 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -10,5 +10,38 @@
1010
[![Coveralls](https://coveralls.io/repos/github/JuliaGeometry/MeshIntegrals.jl/badge.svg?branch=main)](https://coveralls.io/github/JuliaGeometry/MeshIntegrals.jl?branch=main)
1111
[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
1212

13-
**MeshIntegrals.jl** is a Julia library that leverages differential forms to implement fast
14-
and easy numerical integration of field equations over geometric domains.
13+
14+
**MeshIntegrals.jl** uses differential forms to enable fast and easy numerical integration of arbitrary integrand functions over domains defined via [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl) geometries. This is achieved using:
15+
- Gauss-Legendre quadrature rules from [**FastGaussQuadrature.jl**](https://github.com/JuliaApproximation/FastGaussQuadrature.jl): `GaussLegendre(n)`
16+
- H-adaptive Gauss-Kronrod quadrature rules from [**QuadGK.jl**](https://github.com/JuliaMath/QuadGK.jl): `GaussKronrod(kwargs...)`
17+
- H-adaptive cubature rules from [**HCubature.jl**](https://github.com/JuliaMath/HCubature.jl): `HAdaptiveCubature(kwargs...)`
18+
19+
These solvers have support for integrand functions that produce scalars, vectors, and [**Unitful.jl**](https://github.com/PainterQubits/Unitful.jl) `Quantity` types. While HCubature.jl does not natively support `Quantity` type integrands, this package provides a compatibility layer to enable this feature.
20+
21+
## Usage
22+
23+
### Basic
24+
25+
```julia
26+
integral(f, geometry)
27+
```
28+
Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry`. A default integration method will be automatically selected according to the geometry: `GaussKronrod()` for 1D, and `HAdaptiveCubature()` for all others.
29+
30+
```julia
31+
integral(f, geometry, rule)
32+
```
33+
Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry` using the specified integration rule, e.g. `GaussKronrod()`.
34+
35+
Additionally, several optional keyword arguments are defined in [the API](https://juliageometry.github.io/MeshIntegrals.jl/stable/api/) to provide additional control over the integration mechanics.
36+
37+
### Aliases
38+
39+
```julia
40+
lineintegral(f, geometry)
41+
surfaceintegral(f, geometry)
42+
volumeintegral(f, geometry)
43+
```
44+
Alias functions are provided for convenience. These are simply wrappers for `integral` that also validate that the provided `geometry` has the expected number of parametric dimensions. Like with `integral`, a `rule` can also optionally be specified as a third argument.
45+
- `lineintegral` is used for curve-like geometries or polytopes (e.g. `Segment`, `Ray`, `BezierCurve`, `Rope`, etc)
46+
- `surfaceintegral` is used for surfaces (e.g. `Disk`, `Sphere`, `CylinderSurface`, etc)
47+
- `volumeintegral` is used for (3D) volumes (e.g. `Ball`, `Cone`, `Torus`, etc)

docs/src/triangle.md renamed to docs/src/specializations.md

Lines changed: 15 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,18 @@
1-
# Integrating a Triangle
1+
# [Specializations](@id specializations)
2+
3+
There are several notable exceptions to how Meshes.jl defines [parametric functions](@ref how-parametric).
4+
- `Meshes.ConeSurface` is essentially a composite type and has a parametric function that only maps the conical portion of the geometry, so the `Meshes.Disk` base element has to be integrated separately.
5+
- `Meshes.CylinderSurface` is essentially a composite type and has a parametric function that only maps the cylindrical portion of the geometry, so the `Meshes.Disk` element has to be integrated separately.
6+
- `Meshes.FrustumSurface` is essentially a composite type and has a parametric function that only maps the cylindrical portion of the geometry, so the top and bottom `Meshes.Disk` elements have to be integrated separately.
7+
- `Meshes.Line` represents a line of infinite length that passes through two points, and it has a parametric function that is valid on the domain $(-\infty, \infty)$.
8+
- `Meshes.Plane` represents a plane of infinite extent, and it has a parametric function that is valid on the domain $(-\infty, \infty)^2$.
9+
- `Meshes.Ray` represents a line that begins at a point and extends in a particular direction with infinite length, and it has a parametric function that is valid on the domain $[0, \infty)$.
10+
- `Meshes.Ring` 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.
11+
- `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.
12+
- `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$.
13+
- `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$.
14+
15+
## Triangle
216

317
For a specified `Meshes.Triangle` surface with area $A$, let $u$ and $v$ be Barycentric coordinates that span the surface.
418
```math

src/MeshIntegrals.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -9,10 +9,10 @@ import LinearAlgebra
99
import QuadGK
1010
import Unitful
1111

12-
include("utils.jl")
13-
1412
include("differentiation.jl")
15-
export jacobian
13+
export DifferentiationMethod, Analytical, FiniteDifference, jacobian
14+
15+
include("utils.jl")
1616

1717
include("integration_rules.jl")
1818
export IntegrationRule, GaussKronrod, GaussLegendre, HAdaptiveCubature

0 commit comments

Comments
 (0)