|
11 | 11 | [](https://github.com/JuliaTesting/Aqua.jl) |
12 | 12 |
|
13 | 13 |
|
14 | | -This package implements methods for numerically-computing integrals over geometric polytopes |
15 | | -from [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl) using the following `::IntegrationAlgorithms`: |
| 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: |
16 | 15 | - Gauss-Legendre quadrature rules from [**FastGaussQuadrature.jl**](https://github.com/JuliaApproximation/FastGaussQuadrature.jl): `GaussLegendre(n)` |
17 | 16 | - H-adaptive Gauss-Kronrod quadrature rules from [**QuadGK.jl**](https://github.com/JuliaMath/QuadGK.jl): `GaussKronrod(kwargs...)` |
18 | 17 | - H-adaptive cubature rules from [**HCubature.jl**](https://github.com/JuliaMath/HCubature.jl): `HAdaptiveCubature(kwargs...)` |
19 | 18 |
|
20 | | -Functions available: |
21 | | -- `integral(f, ::Geometry, ::IntegrationAlgorithm)`: integrates a function `f` over a domain defined by `geometry` using a particular `::IntegrationAlgorithm` |
22 | | -- `lineintegral`, `surfaceintegral`, and `volumeintegral` are available as aliases for `integral` that first verify that `geometry` has the appropriate number of parametric dimensions |
| 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. |
23 | 20 |
|
24 | | -# Example Usage |
| 21 | +## Usage |
| 22 | + |
| 23 | +```julia |
| 24 | +integral(f, geometry) |
| 25 | +``` |
| 26 | +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. |
| 27 | + |
| 28 | +```julia |
| 29 | +integral(f, geometry, algorithm, FP=Float64) |
| 30 | +``` |
| 31 | +Performs a numerical integration of some integrand function `f(p::Meshes.Point)` over the domain specified by `geometry` using the specified integration algorithm, e.g. `GaussKronrod()`. |
| 32 | + |
| 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. |
| 34 | + |
| 35 | +```julia |
| 36 | +lineintegral(f, geometry) |
| 37 | +surfaceintegral(f, geometry) |
| 38 | +volumeintegral(f, geometry) |
| 39 | +``` |
| 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 `algorithm` 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 | + |
| 45 | +# Demo |
25 | 46 |
|
26 | 47 | ```julia |
27 | 48 | using Meshes |
28 | 49 | using MeshIntegrals |
| 50 | +using Unitful |
29 | 51 |
|
30 | | -# Define a unit circle on the xy-plane |
31 | | -origin = Point(0,0,0) |
32 | | -ẑ = Vec(0,0,1) |
33 | | -xy_plane = Plane(origin,ẑ) |
34 | | -unit_circle_xy = Circle(xy_plane, 1.0) |
35 | | - |
36 | | -# Approximate unit_circle_xy with a high-order Bezier curve |
37 | | -unit_circle_bz = BezierCurve( |
38 | | - [Point(cos(t), sin(t), 0.0) for t in range(0,2pi,length=361)] |
| 52 | +# Define a path that approximates a sine-wave on the xy-plane |
| 53 | +mypath = BezierCurve( |
| 54 | + [Point(t*u"m", sin(t)*u"m", 0.0u"m") for t in range(-pi, pi, length=361)] |
39 | 55 | ) |
40 | 56 |
|
41 | | -# A Real-valued function |
42 | | -f(x, y, z) = abs(x + y) |
43 | | -f(p) = f(to(p)...) |
| 57 | +# Map f(::Point) -> f(x, y, z) in unitless coordinates |
| 58 | +f(p::Meshes.Point) = f(ustrip(to(p))...) |
44 | 59 |
|
45 | | -integral(f, unit_circle_xy, GaussKronrod()) |
46 | | - # 0.000170 seconds (5.00 k allocations: 213.531 KiB) |
47 | | - # ans == 5.656854249525293 m^2 |
| 60 | +# Integrand function in units of Ohms/meter |
| 61 | +f(x, y, z) = (1 / sqrt(1 + cos(x)^2)) * u"Ω/m" |
48 | 62 |
|
49 | | -integral(f, unit_circle_bz, GaussKronrod()) |
50 | | - # 0.017122 seconds (18.93 k allocations: 78.402 MiB) |
51 | | - # ans = 5.551055333711397 m^2 |
| 63 | +integral(f, mypath) |
| 64 | +# -> Approximately 2*Pi Ω |
52 | 65 | ``` |
0 commit comments