Skip to content

Commit 3801370

Browse files
mikeingoldmikeingold
authored
Consolidate onto Jacobian-based integration methods (#18)
* Add a Jacobian function * Add an integral method for Torus using HCubature * Add a torus test * Bugfix * Adjust arg type since hcubature passes SVector * Add placeholder methods for Torus * Bump version and update support matrix * Add generalized internal method for 2D integral using HCubature * Extend usage of internal method * Continue extending generalized method * Add support for integral on CylinderSurface using hcubature * Disable failing CylinderSurface method, rename worker * Implement generalized 2D GaussLegendre * Bugfixes * Disable test for CylinderSurface with hcubature * Extend application of generalized method * Add methods for ParaboloidSurface * Add placeholder for ParaboloidSurface using GaussKronrod * Implement generalized 2D GaussKronrod * Update matrix * Extend generalized GaussKronrod method * Continue extending * Add generalized 1D GaussLegendre * Bugfix * Extending generalized method * Add generalized 1D GaussKronrod * Bugfix * Add generalized 1D hcubature * Add generalized 3D hcubature * Bugfix * Strip old code * Add generalized 3D GaussLegendre * Update docs * Bugfix * Update README plans
1 parent 47af3f0 commit 3801370

File tree

7 files changed

+294
-253
lines changed

7 files changed

+294
-253
lines changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "MeshIntegrals"
22
uuid = "dadec2fd-bbe0-4da4-9dbe-476c782c8e47"
33
authors = ["Mike Ingold <[email protected]>"]
4-
version = "0.10.0"
4+
version = "0.11.0"
55

66
[deps]
77
FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"

‎README.md‎

Lines changed: 11 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -30,7 +30,7 @@ Methods are tested to ensure compatibility with
3030
### Integral
3131
| Geometry | Gauss-Legendre | Gauss-Kronrod | H-Adaptive Cubature |
3232
|----------|----------------|---------------|---------------------|
33-
| `Meshes.Box{Dim,T}` | :x: | :x: | :x: |
33+
| `Meshes.Box{Dim>3,T}` | :x: | :x: | :x: |
3434

3535
### Line Integral
3636
| Geometry | Gauss-Legendre | Gauss-Kronrod | H-Adaptive Cubature |
@@ -51,10 +51,10 @@ Methods are tested to ensure compatibility with
5151
| `Meshes.Box{2,T}` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
5252
| `Meshes.CylinderSurface` | :x: | :white_check_mark: | :x: |
5353
| `Meshes.Disk` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
54-
| `Meshes.ParaboloidSurface` | :x: | :x: | :x: |
54+
| `Meshes.ParaboloidSurface` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
5555
| `Meshes.Sphere{3,T}` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
5656
| `Meshes.Triangle` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
57-
| `Meshes.Torus` | :x: | :x: | :x: |
57+
| `Meshes.Torus` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
5858
| `Meshes.SimpleMesh` | :x: | :x: | :x: |
5959

6060
### Volume Integral
@@ -95,13 +95,15 @@ fr(p) = fr(p.coords...)
9595

9696
# Plans and Work in Progress
9797

98-
- Short term:
99-
- Implement all methods in the support matrix above
98+
- TODO
10099
- Update Example Usage and benchmarks
101100
- Re-implement all tests for Unitful compatibility
102-
103-
- Longer term:
104-
- Once all methods are implemented, determine which can be consolidated with a more abstract/unioned case
105101
- Implement Aqua.jl tests
106-
- Implement Documenter docs
102+
- Improve README documentation in lieu of Documenter
107103
- Implement Monte Carlo integration methods
104+
- Implement BenchmarkTools test suite to quantify performance
105+
- Continue working to generalize and consolidate methods
106+
107+
- Longer term plans
108+
- Assess impact of upcoming Meshes.jl CRS refactor
109+
- Once functionality is established and stable, evaluate transition to JuliaGeometry organization or direct absorption into Meshes.jl

‎src/integral_line.jl‎

Lines changed: 53 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,25 @@ end
3434
# Gauss-Legendre
3535
################################################################################
3636

37+
# Generalized method
38+
function _integral_1d(f, geometry, settings::GaussLegendre)
39+
# Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
40+
xs, ws = gausslegendre(settings.n)
41+
42+
# Change of variables: x [-1,1] ↦ t [0,1]
43+
t(x) = 0.5x + 0.5
44+
point(x) = geometry(t(x))
45+
46+
function paramfactor(x)
47+
J = jacobian(geometry,[t(x)])
48+
return norm(J[1])
49+
end
50+
51+
# Integrate f along the geometry and apply a domain-correction factor for [-1,1] ↦ [0, 1]
52+
integrand((w,x)) = w * f(point(x)) * paramfactor(x)
53+
return 0.5 * sum(integrand, zip(ws, xs))
54+
end
55+
3756
"""
3857
integral(f, curve::Meshes.BezierCurve, ::GaussLegendre; alg=Meshes.Horner())
3958
@@ -73,15 +92,7 @@ function integral(
7392
# A Box{1,T} is definitionally embedded in 1D-space
7493
_validate_integrand(f,1,T)
7594

76-
# Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
77-
xs, ws = gausslegendre(settings.n)
78-
79-
# Change of variables: x [-1,1] ↦ t [0,1]
80-
t(x) = 0.5x + 0.5
81-
point(x) = line(t(x))
82-
83-
# Integrate f along the line and apply a domain-correction factor for [-1,1] ↦ [0, length]
84-
return 0.5 * length(line) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs))
95+
return _integral_1d(f, line, settings)
8596
end
8697

8798
function integral(
@@ -93,16 +104,7 @@ function integral(
93104
# A Circle is definitionally embedded in 3D-space
94105
_validate_integrand(f,3,T)
95106

96-
# Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
97-
xs, ws = gausslegendre(settings.n)
98-
99-
# Change of variables: x [-1,1] ↦ t [0,1]
100-
t(x) = 0.5x + 0.5
101-
point(x) = circle(t(x))
102-
103-
# Integrate f along the circle's rim and apply a domain-correction
104-
# factor for [-1,1] ↦ [0, circumference]
105-
return 0.5 * length(circle) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs))
107+
return _integral_1d(f, circle, settings)
106108
end
107109

108110
function integral(
@@ -135,15 +137,7 @@ function integral(
135137
# Validate the provided integrand function
136138
_validate_integrand(f,Dim,T)
137139

138-
# Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
139-
xs, ws = gausslegendre(settings.n)
140-
141-
# Change of variables: x [-1,1] ↦ t [0,1]
142-
t(x) = 0.5x + 0.5
143-
point(x) = segment(t(x))
144-
145-
# Integrate f along the line and apply a domain-correction factor for [-1,1] ↦ [0, length]
146-
return 0.5 * length(segment) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs))
140+
return _integral_1d(f, segment, settings)
147141
end
148142

149143
function integral(
@@ -155,23 +149,25 @@ function integral(
155149
# A Sphere{2,T} is simply a circle in 2D-space
156150
_validate_integrand(f,2,T)
157151

158-
# Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
159-
xs, ws = gausslegendre(settings.n)
160-
161-
# Change of variables: x [-1,1] ↦ t [0,1]
162-
t(x) = 0.5x + 0.5
163-
point(x) = circle(t(x))
164-
165-
# Integrate f along the circle's rim and apply a domain-correction
166-
# factor for [-1,1] ↦ [0, circumference]
167-
return 0.5 * length(circle) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs))
152+
return _integral_1d(f, circle, settings)
168153
end
169154

170155

171156
################################################################################
172157
# Gauss-Kronrod
173158
################################################################################
174159

160+
# Generalized method
161+
function _integral_1d(f, geometry, settings::GaussKronrod)
162+
function paramfactor(t)
163+
J = jacobian(geometry,[t])
164+
return norm(J[1])
165+
end
166+
167+
integrand(t) = f(geometry(t)) * paramfactor(t)
168+
return QuadGK.quadgk(integrand, 0, 1; settings.kwargs...)[1]
169+
end
170+
175171
"""
176172
integral(f, curve::BezierCurve, ::GaussKronrod; alg=Horner(), kws...)
177173
@@ -205,9 +201,7 @@ function integral(
205201
# A Box is definitionally embedded in 1D-space
206202
_validate_integrand(f,1,T)
207203

208-
len = length(line)
209-
point(t) = line(t)
210-
return QuadGK.quadgk(t -> len * f(point(t)), 0, 1; settings.kwargs...)[1]
204+
return _integral_1d(f, line, settings)
211205
end
212206

213207
function integral(
@@ -219,9 +213,7 @@ function integral(
219213
# A Circle is definitionally embedded in 3D-space
220214
_validate_integrand(f,3,T)
221215

222-
len = length(circle)
223-
point(t) = circle(t)
224-
return QuadGK.quadgk(t -> len * f(point(t)), 0, 1; settings.kwargs...)[1]
216+
return _integral_1d(f, circle, settings)
225217
end
226218

227219
function integral(
@@ -248,9 +240,7 @@ function integral(
248240
# Validate the provided integrand function
249241
_validate_integrand(f,Dim,T)
250242

251-
len = length(segment)
252-
point(t) = segment(t)
253-
return QuadGK.quadgk(t -> len * f(point(t)), 0, 1; settings.kwargs...)[1]
243+
return _integral_1d(f, segment, settings)
254244
end
255245

256246
function integral(
@@ -262,16 +252,25 @@ function integral(
262252
# A Sphere{2,T} is simply a circle in 2D-space
263253
_validate_integrand(f,2,T)
264254

265-
len = length(circle)
266-
point(t) = circle(t)
267-
return QuadGK.quadgk(t -> len * f(point(t)), 0, 1; settings.kwargs...)[1]
255+
return _integral_1d(f, circle, settings)
268256
end
269257

270258

271259
################################################################################
272260
# HCubature
273261
################################################################################
274262

263+
# Generalized method
264+
function _integral_1d(f, geometry, settings::HAdaptiveCubature)
265+
function paramfactor(t)
266+
J = jacobian(geometry,t)
267+
return norm(J[1])
268+
end
269+
270+
integrand(t) = f(geometry(t[1])) * paramfactor(t)
271+
return HCubature.hcubature(integrand, [0], [1]; settings.kwargs...)[1]
272+
end
273+
275274
"""
276275
integral(f, curve::BezierCurve, ::HAdaptiveCubature; alg=Horner(), kws...)
277276
@@ -305,9 +304,7 @@ function integral(
305304
# A Box is definitionally embedded in 1D-space
306305
_validate_integrand(f,1,T)
307306

308-
len = length(line)
309-
point(t) = line(t)
310-
return hcubature(t -> len * f(point(t[1])), [0], [1]; settings.kwargs...)[1]
307+
return _integral_1d(f, line, settings)
311308
end
312309

313310
function integral(
@@ -319,9 +316,7 @@ function integral(
319316
# A Circle is definitionally embedded in 3D-space
320317
_validate_integrand(f,3,T)
321318

322-
len = length(circle)
323-
point(t) = circle(t)
324-
return hcubature(t -> len * f(point(t[1])), [0], [1]; settings.kwargs...)[1]
319+
return _integral_1d(f, circle, settings)
325320
end
326321

327322
function integral(
@@ -352,9 +347,7 @@ function integral(
352347
# Validate the provided integrand function
353348
_validate_integrand(f,Dim,T)
354349

355-
len = length(segment)
356-
point(t) = segment(t)
357-
return hcubature(t -> len * f(point(t[1])), [0], [1]; settings.kwargs...)[1]
350+
return _integral_1d(f, segment, settings)
358351
end
359352

360353
function integral(
@@ -366,7 +359,5 @@ function integral(
366359
# A Sphere{2,T} is simply a circle in 2D-space
367360
_validate_integrand(f,2,T)
368361

369-
len = length(circle)
370-
point(t) = circle(t)
371-
return hcubature(t -> len * f(point(t[1])), [0], [1]; settings.kwargs...)[1]
362+
return _integral_1d(f, circle, settings)
372363
end

0 commit comments

Comments
 (0)