Skip to content

Commit feb56fd

Browse files
mikeingoldmikeingold
authored
Updated conic integration methods (#51)
* Try out new generic method for integrating CylinderSurface * Fix ambiguity bug, add methods for ConeSurface and FrustumSurface * Fix errors in FrustumSurface integrals * Update Support Matrix for Cone and ConeSurface * Add tests for Cone * Update Meshes to latest version * Push CRS and Meshes to newest versions * Attempted bugfix for currently unsatisfiable requirements * Bugfix in test Project - bump Meshes version * Add custom tests for Cone without Meshes.measure * Update Support Matrix * Bugfix - wrong name * Switch Cone GaussKronrod tests to test_throws since not supported * Create a generic integration method for volumes with GaussKronrod that throws not supported error * Add ColPrac and move license to front * Add custom tests for ConeSurface * Add custom tests for FrustumSurface * Bugfix * Bugfix * Remove unused type spec * Bump versions * Bugfix * Remove explicit support for FrustumSurface * Remove redundant methods (never called)
1 parent 0fa6e13 commit feb56fd

File tree

6 files changed

+187
-119
lines changed

6 files changed

+187
-119
lines changed

‎Project.toml‎

Lines changed: 3 additions & 3 deletions
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.12.2"
4+
version = "0.13.0"
55

66
[deps]
77
CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
@@ -13,11 +13,11 @@ QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
1313
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
1414

1515
[compat]
16-
CoordRefSystems = "0.10, 0.11"
16+
CoordRefSystems = "0.12, 0.13"
1717
FastGaussQuadrature = "1"
1818
HCubature = "1.5"
1919
LinearAlgebra = "1"
20-
Meshes = "0.47 - 0.49"
20+
Meshes = "0.50"
2121
QuadGK = "2"
2222
Unitful = "1"
2323
julia = "1.6"

‎README.md‎

Lines changed: 10 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,13 @@
11
# MeshIntegrals.jl
22

3+
[![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT)
4+
[![ColPrac](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet?style=flat-square)](https://github.com/SciML/ColPrac)
5+
36
[![Build Status](https://github.com/mikeingold/MeshIntegrals.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/mikeingold/MeshIntegrals.jl/actions/workflows/CI.yml?query=branch%3Amain)
47
[![codecov](https://codecov.io/gh/mikeingold/MeshIntegrals.jl/graph/badge.svg)](https://codecov.io/gh/mikeingold/MeshIntegrals.jl)
58
[![Coveralls](https://coveralls.io/repos/github/mikeingold/MeshIntegrals.jl/badge.svg?branch=main)](https://coveralls.io/github/mikeingold/MeshIntegrals.jl?branch=main)
69
[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
7-
[![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT)
10+
811

912
This package implements methods for numerically-computing integrals over geometric polytopes
1013
from [**Meshes.jl**](https://github.com/JuliaGeometry/Meshes.jl) using the following `::IntegrationAlgorithms`:
@@ -51,8 +54,7 @@ integral(f, unit_circle_bz, GaussKronrod())
5154
| Symbol | Meaning |
5255
|--------|---------|
5356
| :white_check_mark: | Implemented, passes tests |
54-
| :x: | Planned but not yet implemented |
55-
| :warning: | Unable to implement: parameterization not available (see [Issue #28](https://github.com/mikeingold/MeshIntegrals.jl/issues/28)) |
57+
| :x: | Not yet supported |
5658
| :stop_sign: | Not supported |
5759

5860
### Integral
@@ -64,15 +66,14 @@ integral(f, unit_circle_bz, GaussKronrod())
6466
| `Box` in `𝔼{1}` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
6567
| `Box` in `𝔼{2}` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
6668
| `Box` in `𝔼{3}` | :white_check_mark: | :stop_sign: | :white_check_mark: |
67-
| `Box` in `𝔼{N}` | :x: | :x: | :x: |
6869
| `Circle` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
69-
| `Cone` | :warning: | :warning: | :warning: |
70-
| `ConeSurface` | :x: | :x: | :x: |
70+
| `Cone` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
71+
| `ConeSurface` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
7172
| `Cylinder` | :white_check_mark: | :stop_sign: | :white_check_mark: |
72-
| `CylinderSurface` | :x: | :white_check_mark: | :white_check_mark: |
73+
| `CylinderSurface` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
7374
| `Disk` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
74-
| `Frustum` | :warning: | :warning: | :warning: |
75-
| `FrustumSurface` | :x: | :x: | :x: |
75+
| `Frustum` | :stop_sign: | :stop_sign: | :stop_sign: |
76+
| `FrustumSurface` | :stop_sign: | :stop_sign: | :stop_sign: |
7677
| `Line` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
7778
| `ParaboloidSurface` | :white_check_mark: | :white_check_mark: | :white_check_mark: |
7879
| `Plane` | :white_check_mark: | :white_check_mark: | :white_check_mark: |

‎src/integral_surface.jl‎

Lines changed: 71 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -33,91 +33,109 @@ function _integral_2d(
3333
return QuadGK.quadgk(v -> ∫₁(v), FP(0), FP(1); settings.kwargs...)[1]
3434
end
3535

36-
function _integral_2d(
37-
f,
38-
geometry,
39-
settings::HAdaptiveCubature,
40-
)
41-
return _integral(f, geometry, settings)
42-
end
43-
44-
function _integral_2d(
45-
f,
46-
geometry,
47-
settings::HAdaptiveCubature,
48-
FP::Type{T}
49-
) where {T<:AbstractFloat}
50-
return _integral(f, geometry, settings, FP)
51-
end
52-
5336

5437
################################################################################
55-
# Specialized Methods for CylinderSurface
38+
# Specialized Methods for ConeSurface, CylinderSurface, and FrustumSurface
5639
################################################################################
5740

5841
function integral(
5942
f::F,
6043
cyl::Meshes.CylinderSurface,
61-
settings::GaussLegendre,
44+
settings::I,
6245
FP::Type{T} = Float64
63-
) where {F<:Function, T<:AbstractFloat}
64-
error("Integrating a CylinderSurface with GaussLegendre not supported.")
65-
# TODO Planned to support in the future
46+
) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat}
47+
# The generic method only parameterizes the sides
48+
sides = _integral_2d(f, cyl, settings, FP)
49+
50+
# Integrate the Disk at the top
51+
disk_top = Meshes.Disk(cyl.top, cyl.radius)
52+
top = _integral_2d(f, disk_top, settings, FP)
53+
54+
# Integrate the Disk at the bottom
55+
disk_bottom = Meshes.Disk(cyl.bot, cyl.radius)
56+
bottom = _integral_2d(f, disk_bottom, settings, FP)
57+
58+
return sides + top + bottom
6659
end
6760

6861
function integral(
6962
f::F,
7063
cyl::Meshes.CylinderSurface,
71-
settings::GaussKronrod,
64+
settings::HAdaptiveCubature,
7265
FP::Type{T} = Float64
7366
) where {F<:Function, T<:AbstractFloat}
74-
# Integrate the curved sides of the CylinderSurface
75-
# \int ( \int f(r̄) dz ) dφ
76-
function sides_inner∫(φ)
77-
sidelength = norm(cyl(φ,FP(1)) - cyl(φ,FP(0)))
78-
return sidelength * QuadGK.quadgk(z -> f(cyl(φ,z)), FP(0), FP(1); settings.kwargs...)[1]
79-
end
80-
sides = (FP(2π) * cyl.radius) .* QuadGK.quadgk(φ -> sides_inner∫(φ), FP(0), FP(1); settings.kwargs...)[1]
67+
# The generic method only parameterizes the sides
68+
sides = _integral(f, cyl, settings, FP)
8169

82-
# Integrate the Disk at the top of the CylinderSurface
70+
# Integrate the Disk at the top
8371
disk_top = Meshes.Disk(cyl.top, cyl.radius)
84-
top = _integral_2d(f, disk_top, settings, FP)
72+
top = _integral(f, disk_top, settings, FP)
8573

86-
# Integrate the Disk at the bottom of the CylinderSurface
74+
# Integrate the Disk at the bottom
8775
disk_bottom = Meshes.Disk(cyl.bot, cyl.radius)
88-
bottom = _integral_2d(f, disk_bottom, settings, FP)
76+
bottom = _integral(f, disk_bottom, settings, FP)
8977

9078
return sides + top + bottom
9179
end
9280

9381
function integral(
9482
f::F,
95-
cyl::Meshes.CylinderSurface,
83+
cone::Meshes.ConeSurface,
84+
settings::I,
85+
FP::Type{T} = Float64
86+
) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat}
87+
# The generic method only parameterizes the sides
88+
sides = _integral_2d(f, cone, settings, FP)
89+
90+
# Integrate the Disk at the base
91+
base = _integral_2d(f, cone.base, settings, FP)
92+
93+
return sides + base
94+
end
95+
96+
function integral(
97+
f::F,
98+
cone::Meshes.ConeSurface,
9699
settings::HAdaptiveCubature,
97100
FP::Type{T} = Float64
98101
) where {F<:Function, T<:AbstractFloat}
99-
Dim = Meshes.paramdim(cyl)
102+
# The generic method only parameterizes the sides
103+
sides = _integral(f, cone, settings, FP)
100104

101-
integrand(t) = f(cyl(t...)) * differential(cyl, t)
105+
# Integrate the Disk at the base
106+
base = _integral(f, cone.base, settings, FP)
102107

103-
# HCubature doesn't support functions that output Unitful Quantity types
104-
# Establish the units that are output by f
105-
testpoint_parametriccoord = fill(FP(0.5),Dim)
106-
integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord))
107-
# Create a wrapper that returns only the value component in those units
108-
uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
109-
# Integrate only the unitless values
110-
value = HCubature.hcubature(uintegrand, zeros(FP,Dim), ones(FP,Dim); settings.kwargs...)[1]
111-
# Reapply units
112-
sides = value .* integrandunits
108+
return sides + base
109+
end
113110

114-
# Integrate the Disk at the top of the CylinderSurface
115-
disk_top = Meshes.Disk(cyl.top, cyl.radius)
116-
top = _integral_2d(f, disk_top, settings, FP)
111+
function integral(
112+
f::F,
113+
frust::Meshes.FrustumSurface,
114+
settings::I,
115+
FP::Type{T} = Float64
116+
) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat}
117+
# The generic method only parameterizes the sides
118+
sides = _integral_2d(f, frust, settings, FP)
117119

118-
# Integrate the Disk at the bottom of the CylinderSurface
119-
disk_bottom = Meshes.Disk(cyl.bot, cyl.radius)
120-
bottom = _integral_2d(f, disk_bottom, settings, FP)
120+
# Integrate the Disks at the top and bottom
121+
top = _integral_2d(f, frust.top, settings, FP)
122+
bottom = _integral_2d(f, frust.bot, settings, FP)
123+
124+
return sides + top + bottom
125+
end
126+
127+
function integral(
128+
f::F,
129+
frust::Meshes.FrustumSurface,
130+
settings::HAdaptiveCubature,
131+
FP::Type{T} = Float64
132+
) where {F<:Function, T<:AbstractFloat}
133+
# The generic method only parameterizes the sides
134+
sides = _integral(f, frust, settings, FP)
135+
136+
# Integrate the Disks at the top and bottom
137+
top = _integral(f, frust.top, settings, FP)
138+
bottom = _integral(f, frust.bot, settings, FP)
121139

122140
return sides + top + bottom
123141
end

‎src/integral_volume.jl‎

Lines changed: 4 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -24,21 +24,14 @@ function _integral_3d(
2424
return FP(1/8) .* sum(integrand, zip(wws,xxs))
2525
end
2626

27+
# Integrating volumes with GaussKronrod not supported by default
2728
function _integral_3d(
2829
f,
2930
geometry,
30-
settings::HAdaptiveCubature,
31-
)
32-
return _integral(f, geometry, settings)
33-
end
34-
35-
function _integral_3d(
36-
f,
37-
geometry,
38-
settings::HAdaptiveCubature,
39-
FP::Type{T}
31+
settings::GaussKronrod,
32+
FP::Type{T} = Float64
4033
) where {T<:AbstractFloat}
41-
return _integral(f, geometry, settings, FP)
34+
error("Integrating this volume type with GaussKronrod not supported.")
4235
end
4336

4437

@@ -77,35 +70,3 @@ function integral(
7770
) where {F<:Function, T<:AbstractFloat}
7871
error("Integrating a Tetrahedron with HAdaptiveCubature not supported.")
7972
end
80-
81-
82-
################################################################################
83-
# Unsupported Placeholders
84-
################################################################################
85-
86-
function integral(
87-
f::F,
88-
ball::Meshes.Ball{Meshes.𝔼{3},CRS,ℒ},
89-
settings::GaussKronrod,
90-
FP::Type{T} = Float64
91-
) where {F<:Function, CRS, ℒ, T<:AbstractFloat}
92-
error("Integrating a Ball in 𝔼{3} with GaussKronrod not supported.")
93-
end
94-
95-
function integral(
96-
f::F,
97-
box::Meshes.Box{Meshes.𝔼{3},CRS},
98-
settings::GaussKronrod,
99-
FP::Type{T} = Float64,
100-
) where {F<:Function, CRS, T<:AbstractFloat}
101-
error("Integrating a Box in 𝔼{3} with GaussKronrod not supported.")
102-
end
103-
104-
function integral(
105-
f::F,
106-
box::Meshes.Cylinder,
107-
settings::GaussKronrod,
108-
FP::Type{T} = Float64
109-
) where {F<:Function, T<:AbstractFloat}
110-
error("Integrating a Cylinder with GaussKronrod not supported.")
111-
end

‎test/Project.toml‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
88

99
[compat]
1010
Aqua = "0.8"
11-
Meshes = "0.49"
12-
MeshIntegrals = "0.12"
11+
Meshes = "0.50"
12+
MeshIntegrals = "0.13"
1313
QuadGK = "2"
1414
Unitful = "1"

0 commit comments

Comments
 (0)