Add integral and localintegral

Project.toml

@@ -11,6 +11,7 @@ CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
 DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
 DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
@@ -38,6 +39,7 @@ CoordRefSystems = "0.19"
 DelaunayTriangulation = "1.0"
 DifferentiationInterface = "0.7"
 Distances = "0.10"
+FastGaussQuadrature = "1.1"
 FiniteDifferences = "0.12"
 LinearAlgebra = "1.10"
 Makie = "0.24"

src/Meshes.jl

@@ -44,6 +44,9 @@ import NearestNeighbors: MinkowskiMetric
 import DifferentiationInterface as DI
 import FiniteDifferences: central_fdm
 
+# Integration API
+import FastGaussQuadrature: gausslegendre
+
 # Transforms API
 import TransformsBase: Transform, →
 import TransformsBase: isrevertible, isinvertible
@@ -104,6 +107,7 @@ include("predicates.jl")
 
 # calculus
 include("differentation.jl")
+include("integration.jl")
 
 # operations
 include("centroid.jl")
@@ -458,6 +462,8 @@ export
   derivative,
   jacobian,
   differential,
+  integral,
+  localintegral,
 
   # centroids
   centroid,

src/boundary.jl

@@ -122,7 +122,7 @@ boundary(::Circle) = nothing
 
 embedboundary(c::Circle) = c
 
-boundary(c::Cylinder) = CylinderSurface(bottom(c), top(c), radius(c))
+boundary(c::Cylinder) = CylinderSurface(plane(bottom(c)), plane(top(c)), radius(c))
 
 embedboundary(c::Cylinder) = boundary(c)
 

src/geometries/polytopes.jl

@@ -157,20 +157,22 @@ function angles(c::Chain)
 end
 
 function (c::Chain)(t)
-  if t < 0 || t > 1
-    throw(DomainError(t, "c(t) is not defined for t outside [0, 1]."))
-  end
   segs = segments(c)
   sums = cumsum(measure.(segs))
   sums /= last(sums)
-  # find k such that sums[k] ≤ t < sums[k + 1]
-  k = searchsortedfirst(sums, t) - 1
-  # select segment s at index k
-  s, _ = iterate(segs, k)
+  # find first k such that t ≤ sums[k]
+  k = searchsortedfirst(sums, clamp(t, zero(t), one(t)))
   # reparametrization of t within s
-  ∑ₖ = iszero(k) ? zero(eltype(sums)) : sums[k]
-  ∑ₖ₊₁ = sums[k + 1]
-  s((t - ∑ₖ) / (∑ₖ₊₁ - ∑ₖ))
+  if isone(k)
+    s = first(segs)
+    Σₖ = sums[k]
+    s(t / Σₖ)
+  else
+    s = first(Iterators.drop(segs, k - 1))
+    Σₖ = sums[k]
+    Σₖ₋₁ = sums[k - 1]
+    s((t - Σₖ₋₁) / (Σₖ - Σₖ₋₁))
+  end
 end
 
 # implementations of Chain

src/geometries/primitives/cylinder.jl

@@ -56,9 +56,13 @@ end
 
 paramdim(::Type{<:Cylinder}) = 3
 
-bottom(c::Cylinder) = c.bot
+bottom(c::Cylinder) = Disk(c.bot, bottomradius(c))
 
-top(c::Cylinder) = c.top
+top(c::Cylinder) = Disk(c.top, topradius(c))
+
+bottomradius(c::Cylinder) = norm(c(1, 0, 0) - c(0, 0, 0))
+
+topradius(c::Cylinder) = norm(c(1, 0, 1) - c(0, 0, 1))
 
 radius(c::Cylinder) = c.radius
 
@@ -76,10 +80,12 @@ Base.isapprox(c₁::Cylinder, c₂::Cylinder; atol=atol(lentype(c₁)), kwargs..
 function (c::Cylinder)(ρ, φ, z)
   ℒ = lentype(c)
   T = numtype(ℒ)
-  t = top(c)
-  b = bottom(c)
-  r = radius(c)
-  a = axis(c)
+  C = crs(c)
+  D = datum(C)
+  b = c.bot
+  t = c.top
+  r = c.radius
+  a = Line(b(0, 0), t(0, 0))
   d = a(T(1)) - a(T(0))
 
   # rotation to align z axis with cylinder axis
@@ -91,8 +97,8 @@ function (c::Cylinder)(ρ, φ, z)
   # project a parametric segment between the top and bottom planes
   ρ′ = T(ρ) * r
   φ′ = T(φ) * 2 * T(π) * u"rad"
-  p₁ = Point(convert(crs(c), Cylindrical(ρ′, φ′, zero(ℒ))))
-  p₂ = Point(convert(crs(c), Cylindrical(ρ′, φ′, norm(d))))
+  p₁ = Point(convert(C, Cylindrical{D}(ρ′, φ′, zero(ℒ))))
+  p₂ = Point(convert(C, Cylindrical{D}(ρ′, φ′, norm(d))))
   l = Line(p₁, p₂) |> Affine(R, o)
   s = Segment(l ∩ b, l ∩ t)
   s(T(z))

src/geometries/primitives/cylindersurface.jl

@@ -56,35 +56,38 @@ end
 
 paramdim(::Type{<:CylinderSurface}) = 2
 
-bottom(c::CylinderSurface) = c.bot
+bottom(c::CylinderSurface) = Disk(c.bot, bottomradius(c))
 
-top(c::CylinderSurface) = c.top
+top(c::CylinderSurface) = Disk(c.top, topradius(c))
+
+bottomradius(c::CylinderSurface) = bottomradius(Cylinder(c.bot, c.top, c.radius))
+
+topradius(c::CylinderSurface) = topradius(Cylinder(c.bot, c.top, c.radius))
 
 radius(c::CylinderSurface) = c.radius
 
-axis(c::CylinderSurface) = Line(bottom(c)(0, 0), top(c)(0, 0))
+axis(c::CylinderSurface) = Line(c.bot(0, 0), c.top(0, 0))
 
 # cylinder is right if axis is aligned with plane normals
 function isright(c::CylinderSurface)
   T = numtype(lentype(c))
   a = axis(c)
   d = a(T(1)) - a(T(0))
-  u = normal(bottom(c))
-  v = normal(top(c))
+  u = normal(c.bot)
+  v = normal(c.top)
   isapproxzero(norm(d × u)) && isapproxzero(norm(d × v))
 end
 
 function hasintersectingplanes(c::CylinderSurface)
-  l = bottom(c) ∩ top(c)
-  !isnothing(l) && evaluate(Euclidean(), axis(c), l) < radius(c)
+  l = c.bot ∩ c.top
+  !isnothing(l) && evaluate(Euclidean(), axis(c), l) < c.radius
 end
 
-==(c₁::CylinderSurface, c₂::CylinderSurface) =
-  bottom(c₁) == bottom(c₂) && top(c₁) == top(c₂) && radius(c₁) == radius(c₂)
+==(c₁::CylinderSurface, c₂::CylinderSurface) = c₁.bot == c₂.bot && c₁.top == c₂.top && c₁.radius == c₂.radius
 
 Base.isapprox(c₁::CylinderSurface, c₂::CylinderSurface; atol=atol(lentype(c₁)), kwargs...) =
-  isapprox(bottom(c₁), bottom(c₂); atol, kwargs...) &&
-  isapprox(top(c₁), top(c₂); atol, kwargs...) &&
-  isapprox(radius(c₁), radius(c₂); atol, kwargs...)
+  isapprox(c₁.bot, c₂.bot; atol, kwargs...) &&
+  isapprox(c₁.top, c₂.top; atol, kwargs...) &&
+  isapprox(c₁.radius, c₂.radius; atol, kwargs...)
 
-(c::CylinderSurface)(φ, z) = Cylinder(bottom(c), top(c), radius(c))(1, φ, z)
+(c::CylinderSurface)(φ, z) = Cylinder(c.bot, c.top, c.radius)(1, φ, z)

src/integration.jl

@@ -0,0 +1,94 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    integral(fun, geom; n=3)
+
+Calculate the integral over the `geom`etry of the `fun`ction that maps
+[`Point`](@ref)s to values in a linear space.
+
+    integral(fun, dom; n=3)
+
+Alternatively, calculate the integral over the `dom`ain (e.g., mesh) by
+summing the integrals for each constituent geometry.
+
+Polynomials of degree up to `2n-1` are integrated exactly.
+
+See also [`localintegral`](@ref).
+"""
+integral(fun, geom::Geometry; n=3) = localintegral(fun ∘ geom, geom; n)
+
+integral(fun, cylsurf::CylinderSurface; n=3) =
+  localintegral(fun ∘ cylsurf, cylsurf; n) + integral(fun, top(cylsurf); n) + integral(fun, bottom(cylsurf); n)
+
+integral(fun, conesurf::ConeSurface; n=3) =
+  localintegral(fun ∘ conesurf, conesurf; n) + integral(fun, base(conesurf); n)
+
+integral(fun, frustumsurf::FrustumSurface; n=3) =
+  localintegral(fun ∘ frustumsurf, frustumsurf; n) +
+  integral(fun, top(frustumsurf); n) +
+  integral(fun, bottom(frustumsurf); n)
+
+integral(fun, multi::Multi; n=3) = sum(integral(fun, geom; n) for geom in parent(multi))
+
+integral(fun, dom::Domain; n=3) = sum(integral(fun, geom; n) for geom in dom)
+
+"""
+    localintegral(fun, geom; n=3)
+
+Calculate the integral over the `geom`etry of the `fun`ction that maps
+parametric coordinates `uvw` to values in a linear space.
+
+Polynomials of degree up to `2n-1` are integrated exactly.
+
+See also [`integral`](@ref).
+"""
+localintegral(fun, geom::Geometry; n=3) = _uvwintegral(fun, geom, n)
+
+# ray is parametrized over [0, ∞] interval
+localintegral(fun, ray::Ray; n=3) = _uvwintegral(fun, ray, n, trans=t -> @. t / (1 - t))
+
+# line is parametrized over [-∞, ∞] interval
+localintegral(fun, line::Line; n=3) = _uvwintegral(fun, line, n, trans=t -> @. log(t) - log(1 - t))
+
+# plane is parametrized over [-∞, ∞] interval
+localintegral(fun, plane::Plane; n=3) = _uvwintegral(fun, plane, n, trans=t -> @. log(t) - log(1 - t))
+
+# triangle is parametrized with barycentric coordinates
+localintegral(fun, tri::Triangle; n=3) = _uvwintegral(fun, tri, n, trans=t -> (t[1] * t[2], t[2] - t[1] * t[2]))
+
+# specialize quadrangle for performance
+localintegral(fun, quad::Quadrangle; n=3) = _uvwintegral(fun, quad, n)
+
+# tetrahedron is parametrized with barycentric coordinates
+localintegral(fun, tetra::Tetrahedron; n=3) =
+  _uvwintegral(fun, tetra, n, trans=t -> (t[2] * t[3] - t[1] * t[2] * t[3], t[1] * t[2] * t[3], t[3] - t[2] * t[3]))
+
+# -----------------
+# HELPER FUNCTIONS
+# -----------------
+
+function _uvwintegral(fun, geom, n; trans=identity)
+  # parametric dimension and number type
+  N = paramdim(geom)
+  T = numtype(lentype(geom))
+
+  # Gauss-Legendre quadrature points and weights
+  ts, ws = gausslegendre(n)
+  tgrid = Iterators.product(ntuple(_ -> T.(ts), N)...)
+  wgrid = Iterators.product(ntuple(_ -> T.(ws), N)...)
+
+  # map quadrature points in [-1, 1] to parametric coordinates in [0, 1],
+  # and then map parametric coordinates in [0, 1] to uvw parametrization
+  g = trans ∘ (t -> @. (t + 1) / 2)
+
+  # compute integral with change of variable and differential element
+  Σwᵢfᵢ = sum(zip(tgrid, wgrid)) do (t, w)
+    uvw = g(t)
+    prod(w) * fun(uvw...) * differential(geom, uvw)
+  end
+
+  # adjust for change of variable
+  Σwᵢfᵢ / 2^N
+end