diff --git a/Project.toml b/Project.toml
index c2cdd9450..e5ee4cd52 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,7 +11,9 @@ 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"
+IntegrationInterface = "03ca86a4-b7cd-4b10-a38f-d7a84f40de1e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -38,7 +40,9 @@ CoordRefSystems = "0.19"
 DelaunayTriangulation = "1.0"
 DifferentiationInterface = "0.7"
 Distances = "0.10"
+FastGaussQuadrature = "1.1"
 FiniteDifferences = "0.12"
+IntegrationInterface = "0.1.0"
 LinearAlgebra = "1.10"
 Makie = "0.24"
 NearestNeighbors = "0.4"
diff --git a/src/Meshes.jl b/src/Meshes.jl
index df39176fe..3b2b69bd5 100644
--- a/src/Meshes.jl
+++ b/src/Meshes.jl
@@ -44,6 +44,10 @@ import NearestNeighbors: MinkowskiMetric
 import DifferentiationInterface as DI
 import FiniteDifferences: central_fdm
 
+# Integration API
+import IntegrationInterface as II
+import FastGaussQuadrature: gausslegendre
+
 # Transforms API
 import TransformsBase: Transform, →
 import TransformsBase: isrevertible, isinvertible
@@ -104,6 +108,7 @@ include("predicates.jl")
 
 # calculus
 include("differentation.jl")
+include("integration.jl")
 
 # operations
 include("centroid.jl")
@@ -458,6 +463,8 @@ export
   derivative,
   jacobian,
   differential,
+  integral,
+  localintegral,
 
   # centroids
   centroid,
diff --git a/src/boundary.jl b/src/boundary.jl
index 3fefdf06d..13d0fe776 100644
--- a/src/boundary.jl
+++ b/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)
 
diff --git a/src/geometries/polytopes.jl b/src/geometries/polytopes.jl
index 222c826a2..327aeed72 100644
--- a/src/geometries/polytopes.jl
+++ b/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
diff --git a/src/geometries/primitives/cylinder.jl b/src/geometries/primitives/cylinder.jl
index 17a884cdc..105a88bd3 100644
--- a/src/geometries/primitives/cylinder.jl
+++ b/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))
diff --git a/src/geometries/primitives/cylindersurface.jl b/src/geometries/primitives/cylindersurface.jl
index eb0536845..8789c76ef 100644
--- a/src/geometries/primitives/cylindersurface.jl
+++ b/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)
diff --git a/src/integration.jl b/src/integration.jl
new file mode 100644
index 000000000..4ccb0cd66
--- /dev/null
+++ b/src/integration.jl
@@ -0,0 +1,100 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+# default integration method
+const GAUSSLEGENDRE = II.Backend.Quadrature(gausslegendre(3))
+
+"""
+    integral(fun, geom[, method]) 
+
+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[, method])
+
+Alternatively, calculate the integral over the `dom`ain (e.g., mesh) by
+summing the integrals for each constituent geometry.
+
+By default, use Gauss-Legendre quadrature rule with `n` nodes so that
+polynomials of degree up to `2n-1` are integrated exactly.
+
+See also [`localintegral`](@ref).
+"""
+integral(fun, geom::Geometry, method=GAUSSLEGENDRE) = localintegral(fun ∘ geom, geom, method)
+
+# cylinder surface is the union of the curved surface and the top and bottom disks
+integral(fun, cylsurf::CylinderSurface, method=GAUSSLEGENDRE) =
+  localintegral(fun ∘ cylsurf, cylsurf, method) +
+  integral(fun, top(cylsurf), method) +
+  integral(fun, bottom(cylsurf), method)
+
+# cone surface is the union of the curved surface and the base disk
+integral(fun, conesurf::ConeSurface, method=GAUSSLEGENDRE) =
+  localintegral(fun ∘ conesurf, conesurf, method) + integral(fun, base(conesurf), method)
+
+# frustum surface is the union of the curved surface and the top and bottom disks
+integral(fun, frustumsurf::FrustumSurface, method=GAUSSLEGENDRE) =
+  localintegral(fun ∘ frustumsurf, frustumsurf, method) +
+  integral(fun, top(frustumsurf), method) +
+  integral(fun, bottom(frustumsurf), method)
+
+# multi-geometry is the union of its constituent geometries
+integral(fun, multi::Multi, method=GAUSSLEGENDRE) = sum(integral(fun, geom, method) for geom in parent(multi))
+
+# domain is the union of its constituent geometries
+integral(fun, dom::Domain, method=GAUSSLEGENDRE) = sum(integral(fun, geom, method) for geom in dom)
+
+"""
+    localintegral(fun, geom[, method])
+
+Calculate the integral over the `geom`etry of the `fun`ction that maps
+parametric coordinates `uvw` to values in a linear space.
+
+By default, use Gauss-Legendre quadrature rule with `n` nodes so that
+polynomials of degree up to `2n-1` are integrated exactly.
+
+See also [`integral`](@ref).
+"""
+function localintegral(fun, geom::Geometry, method=GAUSSLEGENDRE)
+  # integrand is equal to function times differential element
+  integrand = uvw -> fun(uvw...) * differential(geom, uvw)
+
+  # domain of integration varies depending on geometry
+  domain = ∫dom(geom)
+
+  # integral of function times differential element
+  I = II.integral(integrand, domain; backend=method)
+
+  # perform numerical integration
+  I()
+end
+
+function ∫dom(geom::Geometry)
+  N = paramdim(geom)
+  T = numtype(lentype(geom))
+  a = ntuple(_ -> zero(T), N)
+  b = ntuple(_ -> one(T), N)
+  II.Domain.Box(a, b)
+end
+
+function ∫dom(ray::Ray)
+  T = numtype(lentype(ray))
+  a = (zero(T),)
+  b = (II.Infinity(one(T)),)
+  II.Domain.Box(a, b)
+end
+
+function ∫dom(line::Line)
+  T = numtype(lentype(line))
+  a = (-II.Infinity(one(T)),)
+  b = (II.Infinity(one(T)),)
+  II.Domain.Box(a, b)
+end
+
+function ∫dom(plane::Plane)
+  T = numtype(lentype(plane))
+  a = (-II.Infinity(one(T)), -II.Infinity(one(T)))
+  b = (II.Infinity(one(T)), II.Infinity(one(T)))
+  II.Domain.Box(a, b)
+end
diff --git a/test/integration.jl b/test/integration.jl
new file mode 100644
index 000000000..83902cc5d
--- /dev/null
+++ b/test/integration.jl
@@ -0,0 +1,339 @@
+@testitem "integral" setup = [Setup] begin
+  # Ray
+  a = cart(0, 0, 0)
+  v = vector(1, 1, 1)
+  ray = Ray(a, v)
+  function funray(p)
+    r = ustrip(u"m", norm(to(p)))
+    exp(-r^2) * u"A"
+  end
+  solution = sqrt(T(π)) / 2 * u"A*m"
+  @test_broken integral(funray, ray) ≈ solution
+
+  # Line
+  a = cart(0, 0, 0)
+  b = cart(1, 1, 1)
+  line = Line(a, b)
+  function funline(p)
+    r = ustrip(u"m", norm(to(p)))
+    exp(-r^2) * u"A"
+  end
+  solution = sqrt(T(π)) * u"A*m"
+  @test_broken integral(funline, line) ≈ solution
+
+  # Bezier Curve
+  bezier = BezierCurve([cart(t, sin(t), 0) for t in range(-π, π, length=361)])
+  function funbezier(p)
+    ux = ustrip(coords(p).x)
+    (1 / sqrt(1 + cos(ux)^2)) * u"Ω"
+  end
+  solution = T(2π) * u"Ω*m"
+  @test integral(funbezier, bezier, n=10) ≈ solution rtol = 1e-2
+
+  # Plane
+  p = cart(0, 0, 0)
+  v = vector(0, 0, 1)
+  plane = Plane(p, v)
+  function funplane(p)
+    r = ustrip(u"m", norm(to(p)))
+    exp(-r^2) * u"A"
+  end
+  solution = T(π) * u"A*m^2"
+  @test_broken integral(funplane, plane) ≈ solution
+
+  # Box 1D
+  a = T(π)
+  box = Box(cart(0), cart(a))
+  function funbox1(p)
+    x₁ = only(ustrip.(to(p)))
+    √(a^2 - x₁^2) * u"A"
+  end
+  solution = T(π) * a^2 / 4 * u"A*m"
+  @test integral(funbox1, box, n=10) ≈ solution rtol = 1e-3
+
+  # Box 2D
+  a = T(π)
+  box = Box(cart(0, 0), cart(a, a))
+  function funbox2(p)
+    x₁, x₂ = ustrip.(to(p))
+    (√(a^2 - x₁^2) + √(a^2 - x₂^2)) * u"A"
+  end
+  solution = 2a * (T(π) * a^2 / 4) * u"A*m^2"
+  @test integral(funbox2, box, n=10) ≈ solution rtol = 1e-3
+
+  # Box 3D
+  a = T(π)
+  box = Box(cart(0, 0, 0), cart(a, a, a))
+  function funbox3(p)
+    x₁, x₂, x₃ = ustrip.(to(p))
+    (√(a^2 - x₁^2) + √(a^2 - x₂^2) + √(a^2 - x₃^2)) * u"A"
+  end
+  solution = 3a^2 * (T(π) * a^2 / 4) * u"A*m^3"
+  @test integral(funbox3, box, n=10) ≈ solution rtol = 1e-3
+
+  # Ball 2D
+  origin = cart(0, 0)
+  radius = T(2.8)
+  ball = Ball(origin, radius)
+  function funball2(p)
+    r = ustrip(u"m", norm(to(p)))
+    exp(-r^2) * u"A"
+  end
+  solution = (T(π) - T(π) * exp(-radius^2)) * u"A*m^2"
+  @test integral(funball2, ball, n=10) ≈ solution rtol = 1e-3
+
+  # Ellipsoid
+  origin = cart(0, 0, 0)
+  R = r₁ = r₂ = r₃ = T(4.1)
+  ellipsoid = Ellipsoid((r₁, r₂, r₃), origin)
+  function funellips(p)
+    x, y, z = ustrip.(u"m", to(p))
+    (z^2) * u"A"
+  end
+  solution = (T(4π) * R^4 / 3) * u"A*m^2"
+  @test integral(funellips, ellipsoid, n=10) ≈ solution rtol = 1e-3
+
+  # Disk
+  center = cart(1, 2, 3)
+  normal = vector(1 / 2, 1 / 2, sqrt(2) / 2)
+  plane = Plane(center, normal)
+  radius = T(2.5)
+  disk = Disk(plane, radius)
+  function fundisk(p)
+    offset = p - center
+    r = ustrip(u"m", norm(offset))
+    exp(-r^2) * u"A"
+  end
+  solution = (T(π) - T(π) * exp(-radius^2)) * u"A*m^2"
+  @test integral(fundisk, disk, n=10) ≈ solution rtol = 1e-3
+
+  # Circle
+  center = cart(1, 2, 3)
+  normal = vector(1 / 2, 1 / 2, sqrt(2) / 2)
+  plane = Plane(center, normal)
+  radius = T(4.4)
+  circle = Circle(plane, radius)
+  function funcircle(p)
+    offset = p - center
+    r = ustrip(u"m", norm(offset))
+    exp(-r^2) * u"A"
+  end
+  solution = T(2π) * radius * exp(-radius^2) * u"A*m"
+  @test integral(funcircle, circle, n=10) ≈ solution rtol = 1e-3
+
+  # Cylinder
+  h = T(8.5)u"m"
+  ρ = T(1.3)u"m"
+  a = cart(0, 0, 0)
+  b = cart(0u"m", 0u"m", h)
+  cyl = Cylinder(a, b, ρ)
+  function funcylinder(p)
+    c = convert(Cylindrical, coords(p))
+    ρ = c.ρ
+    ϕ = c.ϕ
+    z = c.z
+    ρ^(-1) * (ρ + ϕ * u"m" + z) * u"A"
+  end
+  solution = ((T(π) * h * ρ^2) + (T(π) * h^2 * ρ) + (T(2π) * T(π) * u"m" * h * ρ)) * u"A"
+  @test integral(funcylinder, cyl, n=10) ≈ solution
+
+  # CylinderSurface
+  h = T(8.5)u"m"
+  ρ = T(1.3)u"m"
+  a = cart(0, 0, 0)
+  b = cart(0u"m", 0u"m", h)
+  cylsurf = CylinderSurface(a, b, ρ)
+  function funcylsurf(p)
+    c = convert(Cylindrical, coords(p))
+    ρ = c.ρ
+    ϕ = c.ϕ
+    z = c.z
+    ρ^(-1) * (ρ + ϕ * u"m" + z) * u"A"
+  end
+  solution = let
+    A1 = (T(2π) * h * ρ) + (T(π) * ρ^2) + (T(π) * u"m" * ρ * T(2π))
+    A2 = (T(π) * ρ^2) + (T(π) * u"m" * ρ * T(2π))
+    A3 = (T(2π) * h * ρ) + (2T(π)^2 * u"m" * h) + (T(π) * h^2)
+    (A1 + A2 + A3) * u"A"
+  end
+  @test integral(funcylsurf, cylsurf, n=10) ≈ solution
+
+  # Cone
+  r = T(2.5)u"m"
+  h = T(3.5)u"m"
+  origin = cart(0, 0, 0)
+  plane = Plane(origin, vector(0, 0, 1))
+  base = Disk(plane, r)
+  apex = cart(0u"m", 0u"m", h)
+  cone = Cone(base, apex)
+  funcone(p) = T(1.0)u"A"
+  solution = (T(π) * r^2 * h / 3) * u"A"
+  @test integral(funcone, cone) ≈ solution
+
+  # ConeSurface
+  r = T(2.5)u"m"
+  h = T(3.5)u"m"
+  origin = cart(0, 0, 0)
+  plane = Plane(origin, vector(0, 0, 1))
+  base = Disk(plane, r)
+  apex = cart(0u"m", 0u"m", h)
+  conesurf = ConeSurface(base, apex)
+  funconesurf(p) = T(1.0)u"A"
+  solution = ((T(π) * r^2) + (T(π) * r * hypot(h, r))) * u"A"
+  @test integral(funconesurf, conesurf) ≈ solution
+
+  # Frustum
+  r = T(2.5)u"m"
+  h = T(3.5)u"m"
+  origin = cart(0, 0, 0)
+  normal = vector(0, 0, 1)
+  midpoint = cart(0.0u"m", 0.0u"m", h / 2)
+  base = Disk(Plane(origin, normal), r)
+  disk = Disk(Plane(midpoint, normal), r / 2)
+  frustum = Frustum(base, disk)
+  funfrustum(p) = T(1.0)u"A"
+  solution = (T(7) / T(8)) * (T(π) * r^2 * h / T(3)) * u"A"
+  @test integral(funfrustum, frustum) ≈ solution
+
+  # FrustumSurface
+  rbot = T(2.5)u"m"
+  rtop = T(1.25)u"m"
+  height = T(2π) * u"m"
+  origin = cart(0, 0, 0)
+  normal = vector(0, 0, 1)
+  planebot = Plane(origin, normal)
+  diskbot = Disk(planebot, rbot)
+  centertop = cart(0.0u"m", 0.0u"m", height / 2)
+  planetop = Plane(centertop, normal)
+  disktop = Disk(planetop, rtop)
+  frustumsurf = FrustumSurface(diskbot, disktop)
+  funfrustumsurf(p) = T(1.0)u"A"
+  solution = let
+    A1 = T(π) * rbot * hypot(height, rbot)
+    A2 = T(π) * rtop * hypot(height / 2, rtop)
+    A3 = T(π) * rtop^2
+    A4 = T(π) * rbot^2
+    (A1 - A2 + A3 + A4) * u"A"
+  end
+  @test integral(funfrustumsurf, frustumsurf) ≈ solution
+
+  # Segment
+  ϕ = 7T(pi) / 6
+  θ = T(pi) / 3
+  a = cart(0, 0, 0)
+  b = cart(sin(θ) * cos(ϕ), sin(θ) * sin(ϕ), cos(θ))
+  seg = Segment(a, b)
+  ka = T(7.1)
+  kb = T(4.6)
+  function funseg(p)
+    r = ustrip(u"m", norm(to(p)))
+    exp(r * log(ka) + (1 - r) * log(kb)) * u"A"
+  end
+  solution = ((ka - kb) / (log(ka) - log(kb))) * u"A*m"
+  @test integral(funseg, seg) ≈ solution
+
+  # Rope
+  a = cart(0, 0, 0)
+  b = cart(1, 0, 0)
+  c = cart(1, 1, 0)
+  d = cart(1, 1, 1)
+  rope = Rope(a, b, c, d)
+  function funrope(p)
+    x, y, z = ustrip.(to(p))
+    (x + 2y + 3z) * u"A"
+  end
+  solution = T(7.0)u"A*m"
+  @test integral(funrope, rope, n=100) ≈ solution rtol = 1e-2
+
+  # Ring
+  a = cart(0, 0, 0)
+  b = cart(1, 0, 0)
+  c = cart(1, 1, 0)
+  d = cart(1, 1, 1)
+  ring = Ring(a, b, c, d, c, b)
+  function funring(p)
+    x, y, z = ustrip.(to(p))
+    (x + 2y + 3z) * u"A"
+  end
+  solution = T(14.0)u"A*m"
+  @test integral(funring, ring, n=100) ≈ solution rtol = 1e-2
+
+  # PolyArea
+  a, b, c, z = T(0.4), T(0.6), T(1.0), T(0.0)
+  outer = [(z, z), (c, z), (c, c), (z, c)]
+  hole = [(a, a), (a, b), (b, b), (b, a)]
+  poly = PolyArea([outer, hole])
+  function funpoly(p)
+    x, y = ustrip.(u"m", to(p))
+    2x * u"A"
+  end
+  solution = (c^2 - (b - a) * (b^2 - a^2)) * u"A*m^2"
+  @test_broken integral(funpoly, poly) ≈ solution
+
+  # Triangle
+  a = cart(0, 0, 0)
+  b = cart(1, 0, 0)
+  c = cart(0, 1, 0)
+  tri = Triangle(a, b, c)
+  function funtri(p)
+    x, y, z = ustrip.(u"m", to(p))
+    (x + 2y + 3z) * u"A"
+  end
+  solution = T(0.5) * u"A*m^2"
+  @test_broken integral(funtri, tri) ≈ solution
+
+  # Quadrangle
+  quad = Quadrangle(cart(-1.0, 0.0), cart(-1.0, 1.0), cart(1.0, 1.0), cart(1.0, 0.0))
+  function funquad(p)
+    r = ustrip(u"m", norm(to(p)))
+    exp(-r^2) * u"A"
+  end
+  solution = T(π) * T(0.8427007929497149)^2 / 2 * u"A*m^2" # erf(1) = 0.8427007929497149
+  @test integral(funquad, quad, n=10) ≈ solution rtol = 1e-3
+
+  # Tetrahedron
+  a = cart(0, 0, 0)
+  b = cart(1, 0, 0)
+  c = cart(0, 1, 0)
+  d = cart(0, 0, 1)
+  tetra = Tetrahedron(a, b, c, d)
+  function funtetra(p)
+    x, y, z = ustrip.(u"m", to(p))
+    (x + 2y + 3z) * u"A"
+  end
+  solution = T(0.25) * u"A*m^3"
+  @test_broken integral(funtetra, tetra) ≈ solution
+
+  # Hexahedron
+  a = π
+  box = Box(cart(0, 0, 0), cart(a, a, a))
+  hexa = convert(Hexahedron, box)
+  function funhexa(p)
+    x₁, x₂, x₃ = ustrip.(to(p))
+    (√(a^2 - x₁^2) + √(a^2 - x₂^2) + √(a^2 - x₃^2)) * u"A"
+  end
+  solution = 3a^2 * (π * a^2 / 4) * u"A*m^3"
+  @test integral(funhexa, hexa, n=10) ≈ solution rtol = 1e-3
+
+  # Multi
+  box = Box(cart(0, 0), cart(1, 1))
+  ball = Ball(cart(5, 5), T(1))
+  multi = Multi([box, ball])
+  funmulti(p) = sum(to(p))
+  @test integral(funmulti, multi) ≈ integral(funmulti, box) + integral(funmulti, ball)
+
+  # GeometrySet
+  box = Box(cart(0, 0), cart(1, 1))
+  ball = Ball(cart(5, 5), T(1))
+  gset = GeometrySet([box, ball])
+  fungset(p) = sum(to(p))
+  @test integral(fungset, gset) ≈ integral(fungset, box) + integral(fungset, ball)
+
+  # SimpleMesh
+  # TODO:
+
+  # Grid
+  grid = cartgrid(10, 10)
+  fungrid(p) = T(1) * u"A"
+  @test integral(fungrid, grid) ≈ 100 * integral(fungrid, first(grid))
+end
diff --git a/test/primitives.jl b/test/primitives.jl
index b8a95e053..6f965d45b 100644
--- a/test/primitives.jl
+++ b/test/primitives.jl
@@ -863,8 +863,8 @@ end
   @test crs(c) <: Cartesian{NoDatum}
   @test Meshes.lentype(c) == ℳ
   @test radius(c) == T(5) * u"m"
-  @test bottom(c) == Plane(cart(1, 2, 3), vector(0, 0, 1))
-  @test top(c) == Plane(cart(4, 5, 6), vector(0, 0, 1))
+  @test bottom(c) == Disk(Plane(cart(1, 2, 3), vector(0, 0, 1)), Meshes.bottomradius(c))
+  @test top(c) == Disk(Plane(cart(4, 5, 6), vector(0, 0, 1)), Meshes.topradius(c))
   @test axis(c) == Line(cart(1, 2, 3), cart(4, 5, 6))
   @test !isright(c)
   @test measure(c) == volume(c) ≈ T(5)^2 * pi * T(3) * sqrt(T(3)) * u"m^3"
@@ -897,8 +897,8 @@ end
 
   c = Cylinder(cart(0, 0, 0), cart(0, 0, 1), T(1))
   @test radius(c) == T(1) * u"m"
-  @test bottom(c) == Plane(cart(0, 0, 0), vector(0, 0, 1))
-  @test top(c) == Plane(cart(0, 0, 1), vector(0, 0, 1))
+  @test bottom(c) == Disk(Plane(cart(0, 0, 0), vector(0, 0, 1)), T(1))
+  @test top(c) == Disk(Plane(cart(0, 0, 1), vector(0, 0, 1)), T(1))
   @test centroid(c) == cart(0.0, 0.0, 0.5)
   @test axis(c) == Line(cart(0, 0, 0), cart(0, 0, 1))
   @test isright(c)
@@ -943,8 +943,8 @@ end
   @test crs(c) <: Cartesian{NoDatum}
   @test Meshes.lentype(c) == ℳ
   @test radius(c) == T(2) * u"m"
-  @test bottom(c) == Plane(cart(0, 0, 0), vector(0, 0, 1))
-  @test top(c) == Plane(cart(0, 0, 1), vector(0, 0, 1))
+  @test bottom(c) == Disk(Plane(cart(0, 0, 0), vector(0, 0, 1)), T(2))
+  @test top(c) == Disk(Plane(cart(0, 0, 1), vector(0, 0, 1)), T(2))
   @test centroid(c) == cart(0.0, 0.0, 0.5)
   @test axis(c) == Line(cart(0, 0, 0), cart(0, 0, 1))
   @test isright(c)