diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
new file mode 100644
index 00000000..580b7511
--- /dev/null
+++ b/.JuliaFormatter.toml
@@ -0,0 +1 @@
+style = "sciml"
diff --git a/.github/workflows/FormatPR.yml b/.github/workflows/FormatPR.yml
new file mode 100644
index 00000000..8ce8a2d4
--- /dev/null
+++ b/.github/workflows/FormatPR.yml
@@ -0,0 +1,11 @@
+name: Format suggestions
+on:
+  pull_request
+
+jobs:
+  code-style:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: julia-actions/julia-format@v3
+        with:
+          version: '1' # Set `version` to '1.0.54' if you need to use JuliaFormatter.jl v1.0.54 (default: '1')
diff --git a/docs/make.jl b/docs/make.jl
index b4400099..f29f6713 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -2,7 +2,7 @@ using Documenter
 using MeshIntegrals
 
 makedocs(
-    sitename="MeshIntegrals.jl",
+    sitename = "MeshIntegrals.jl",
     pages = [
         "Home" => [
             "About" => "index.md",
@@ -17,5 +17,5 @@ makedocs(
 )
 
 deploydocs(repo = "github.com/mikeingold/MeshIntegrals.jl.git",
-           devbranch = "main",
-           push_preview = true)
+    devbranch = "main",
+    push_preview = true)
diff --git a/src/MeshIntegrals.jl b/src/MeshIntegrals.jl
index 9bedf94c..eefbb991 100644
--- a/src/MeshIntegrals.jl
+++ b/src/MeshIntegrals.jl
@@ -1,35 +1,35 @@
 module MeshIntegrals
-    using CoordRefSystems
-    using LinearAlgebra
-    using Meshes
-    using Unitful
+using CoordRefSystems
+using LinearAlgebra
+using Meshes
+using Unitful
 
-    import FastGaussQuadrature
-    import HCubature
-    import QuadGK
+import FastGaussQuadrature
+import HCubature
+import QuadGK
 
-    include("utils.jl")
-    export jacobian, derivative, unitdirection
+include("utils.jl")
+export jacobian, derivative, unitdirection
 
-    include("integration_rules.jl")
-    export IntegrationRule, GaussKronrod, GaussLegendre, HAdaptiveCubature
+include("integration_rules.jl")
+export IntegrationRule, GaussKronrod, GaussLegendre, HAdaptiveCubature
 
-    include("integral.jl")
-    export integral
+include("integral.jl")
+export integral
 
-    include("integral_aliases.jl")
-    export lineintegral, surfaceintegral, volumeintegral
+include("integral_aliases.jl")
+export lineintegral, surfaceintegral, volumeintegral
 
-    # Integration methods specialized for particular geometries
-    include("specializations/BezierCurve.jl")
-    include("specializations/ConeSurface.jl")
-    include("specializations/CylinderSurface.jl")
-    include("specializations/FrustumSurface.jl")
-    include("specializations/Line.jl")
-    include("specializations/Plane.jl")
-    include("specializations/Ray.jl")
-    include("specializations/Ring.jl")
-    include("specializations/Rope.jl")
-    include("specializations/Tetrahedron.jl")
-    include("specializations/Triangle.jl")
+# Integration methods specialized for particular geometries
+include("specializations/BezierCurve.jl")
+include("specializations/ConeSurface.jl")
+include("specializations/CylinderSurface.jl")
+include("specializations/FrustumSurface.jl")
+include("specializations/Line.jl")
+include("specializations/Plane.jl")
+include("specializations/Ray.jl")
+include("specializations/Ring.jl")
+include("specializations/Rope.jl")
+include("specializations/Tetrahedron.jl")
+include("specializations/Triangle.jl")
 end
diff --git a/src/integral.jl b/src/integral.jl
index 4010988c..b9711481 100644
--- a/src/integral.jl
+++ b/src/integral.jl
@@ -25,9 +25,9 @@ function integral end
 
 # If only f and geometry are specified, select default rule
 function integral(
-    f::F,
-    geometry::G
-) where {F<:Function, G<:Meshes.Geometry}
+        f::F,
+        geometry::G
+) where {F <: Function, G <: Meshes.Geometry}
     N = Meshes.paramdim(geometry)
     rule = (N == 1) ? GaussKronrod() : HAdaptiveCubature()
     _integral(f, geometry, rule)
@@ -35,26 +35,25 @@ end
 
 # with rule and T specified
 function integral(
-    f::F,
-    geometry::G,
-    rule::I,
-    FP::Type{T} = Float64
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationRule, T<:AbstractFloat}
+        f::F,
+        geometry::G,
+        rule::I,
+        FP::Type{T} = Float64
+) where {F <: Function, G <: Meshes.Geometry, I <: IntegrationRule, T <: AbstractFloat}
     _integral(f, geometry, rule, FP)
 end
 
-
 ################################################################################
 #                    Generalized (n-Dimensional) Worker Methods
 ################################################################################
 
 # GaussKronrod
 function _integral(
-    f,
-    geometry,
-    rule::GaussKronrod,
-    FP::Type{T} = Float64
-) where {T<:AbstractFloat}
+        f,
+        geometry,
+        rule::GaussKronrod,
+        FP::Type{T} = Float64
+) where {T <: AbstractFloat}
     # Run the appropriate integral type
     N = Meshes.paramdim(geometry)
     if N == 1
@@ -68,11 +67,11 @@ end
 
 # GaussLegendre
 function _integral(
-    f,
-    geometry,
-    rule::GaussLegendre,
-    FP::Type{T} = Float64
-) where {T<:AbstractFloat}
+        f,
+        geometry,
+        rule::GaussLegendre,
+        FP::Type{T} = Float64
+) where {T <: AbstractFloat}
     N = Meshes.paramdim(geometry)
 
     # Get Gauss-Legendre nodes and weights for a region [-1,1]^N
@@ -81,23 +80,23 @@ function _integral(
     nodes = Iterators.product(ntuple(Returns(xs), N)...)
 
     # Domain transformation: x [-1,1] ↦ t [0,1]
-    t(x) = FP(1//2) * x + FP(1//2)
+    t(x) = FP(1 // 2) * x + FP(1 // 2)
 
     function integrand((weights, nodes))
         ts = t.(nodes)
         prod(weights) * f(geometry(ts...)) * differential(geometry, ts)
     end
 
-    return FP(1//(2^N)) .* sum(integrand, zip(weights, nodes))
+    return FP(1 // (2^N)) .* sum(integrand, zip(weights, nodes))
 end
 
 # HAdaptiveCubature
 function _integral(
-    f,
-    geometry,
-    rule::HAdaptiveCubature,
-    FP::Type{T} = Float64
-) where {T<:AbstractFloat}
+        f,
+        geometry,
+        rule::HAdaptiveCubature,
+        FP::Type{T} = Float64
+) where {T <: AbstractFloat}
     N = Meshes.paramdim(geometry)
 
     integrand(t) = f(geometry(t...)) * differential(geometry, t)
@@ -109,7 +108,7 @@ function _integral(
     # Create a wrapper that returns only the value component in those units
     uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
     # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, zeros(FP,N), ones(FP,N); rule.kwargs...)[1]
+    value = HCubature.hcubature(uintegrand, zeros(FP, N), ones(FP, N); rule.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits
@@ -120,32 +119,32 @@ end
 ################################################################################
 
 function _integral_1d(
-    f,
-    geometry,
-    rule::GaussKronrod,
-    FP::Type{T} = Float64
-) where {T<:AbstractFloat}
+        f,
+        geometry,
+        rule::GaussKronrod,
+        FP::Type{T} = Float64
+) where {T <: AbstractFloat}
     integrand(t) = f(geometry(t)) * differential(geometry, (t))
     return QuadGK.quadgk(integrand, zero(FP), one(FP); rule.kwargs...)[1]
 end
 
 function _integral_2d(
-    f,
-    geometry2d,
-    rule::GaussKronrod,
-    FP::Type{T} = Float64
-) where {T<:AbstractFloat}
-    integrand(u,v) = f(geometry2d(u,v)) * differential(geometry2d, (u,v))
-    ∫₁(v) = QuadGK.quadgk(u -> integrand(u,v), zero(FP), one(FP); rule.kwargs...)[1]
+        f,
+        geometry2d,
+        rule::GaussKronrod,
+        FP::Type{T} = Float64
+) where {T <: AbstractFloat}
+    integrand(u, v) = f(geometry2d(u, v)) * differential(geometry2d, (u, v))
+    ∫₁(v) = QuadGK.quadgk(u -> integrand(u, v), zero(FP), one(FP); rule.kwargs...)[1]
     return QuadGK.quadgk(v -> ∫₁(v), zero(FP), one(FP); rule.kwargs...)[1]
 end
 
 # Integrating volumes with GaussKronrod not supported by default
 function _integral_3d(
-    f,
-    geometry,
-    rule::GaussKronrod,
-    FP::Type{T} = Float64
-) where {T<:AbstractFloat}
+        f,
+        geometry,
+        rule::GaussKronrod,
+        FP::Type{T} = Float64
+) where {T <: AbstractFloat}
     error("Integrating this volume type with GaussKronrod not supported.")
 end
diff --git a/src/integral_aliases.jl b/src/integral_aliases.jl
index fff75031..c919c6c0 100644
--- a/src/integral_aliases.jl
+++ b/src/integral_aliases.jl
@@ -17,9 +17,9 @@ Rule types available:
 - HAdaptiveCubature
 """
 function lineintegral(
-    f::F,
-    geometry::G,
-) where {F<:Function, G<:Meshes.Geometry}
+        f::F,
+        geometry::G
+) where {F <: Function, G <: Meshes.Geometry}
     N = Meshes.paramdim(geometry)
 
     if N == 1
@@ -30,10 +30,10 @@ function lineintegral(
 end
 
 function lineintegral(
-    f::F,
-    geometry::G,
-    rule::I
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationRule}
+        f::F,
+        geometry::G,
+        rule::I
+) where {F <: Function, G <: Meshes.Geometry, I <: IntegrationRule}
     N = Meshes.paramdim(geometry)
 
     if N == 1
@@ -44,11 +44,11 @@ function lineintegral(
 end
 
 function lineintegral(
-    f::F,
-    geometry::G,
-    rule::I,
-    FP::Type{T}
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationRule, T<:AbstractFloat}
+        f::F,
+        geometry::G,
+        rule::I,
+        FP::Type{T}
+) where {F <: Function, G <: Meshes.Geometry, I <: IntegrationRule, T <: AbstractFloat}
     N = Meshes.paramdim(geometry)
 
     if N == 1
@@ -58,7 +58,6 @@ function lineintegral(
     end
 end
 
-
 ################################################################################
 #                              Surface Integral
 ################################################################################
@@ -78,9 +77,9 @@ Algorithm types available:
 - HAdaptiveCubature (default)
 """
 function surfaceintegral(
-    f::F,
-    geometry::G,
-) where {F<:Function, G<:Meshes.Geometry}
+        f::F,
+        geometry::G
+) where {F <: Function, G <: Meshes.Geometry}
     N = Meshes.paramdim(geometry)
 
     if N == 2
@@ -91,10 +90,10 @@ function surfaceintegral(
 end
 
 function surfaceintegral(
-    f::F,
-    geometry::G,
-    rule::I
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationRule}
+        f::F,
+        geometry::G,
+        rule::I
+) where {F <: Function, G <: Meshes.Geometry, I <: IntegrationRule}
     N = Meshes.paramdim(geometry)
 
     if N == 2
@@ -105,11 +104,11 @@ function surfaceintegral(
 end
 
 function surfaceintegral(
-    f::F,
-    geometry::G,
-    rule::I,
-    FP::Type{T}
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationRule, T<:AbstractFloat}
+        f::F,
+        geometry::G,
+        rule::I,
+        FP::Type{T}
+) where {F <: Function, G <: Meshes.Geometry, I <: IntegrationRule, T <: AbstractFloat}
     N = Meshes.paramdim(geometry)
 
     if N == 2
@@ -119,7 +118,6 @@ function surfaceintegral(
     end
 end
 
-
 ################################################################################
 #                              Volume Integral
 ################################################################################
@@ -139,9 +137,9 @@ Algorithm types available:
 - HAdaptiveCubature (default)
 """
 function volumeintegral(
-    f::F,
-    geometry::G,
-) where {F<:Function, G<:Meshes.Geometry}
+        f::F,
+        geometry::G
+) where {F <: Function, G <: Meshes.Geometry}
     N = Meshes.paramdim(geometry)
 
     if N == 3
@@ -152,10 +150,10 @@ function volumeintegral(
 end
 
 function volumeintegral(
-    f::F,
-    geometry::G,
-    rule::I
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationRule}
+        f::F,
+        geometry::G,
+        rule::I
+) where {F <: Function, G <: Meshes.Geometry, I <: IntegrationRule}
     N = Meshes.paramdim(geometry)
 
     if N == 3
@@ -166,11 +164,11 @@ function volumeintegral(
 end
 
 function volumeintegral(
-    f::F,
-    geometry::G,
-    rule::I,
-    FP::Type{T}
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationRule, T<:AbstractFloat}
+        f::F,
+        geometry::G,
+        rule::I,
+        FP::Type{T}
+) where {F <: Function, G <: Meshes.Geometry, I <: IntegrationRule, T <: AbstractFloat}
     N = Meshes.paramdim(geometry)
 
     if N == 3
diff --git a/src/integration_rules.jl b/src/integration_rules.jl
index 1a5e11d6..3ac3544c 100644
--- a/src/integration_rules.jl
+++ b/src/integration_rules.jl
@@ -13,7 +13,7 @@ dimensional geometries; some two- and three-dimensional geometries are additiona
 supported using nested integral solvers with the specified `kwarg` settings.
 """
 struct GaussKronrod <: IntegrationRule
-    kwargs
+    kwargs::Any
     GaussKronrod(; kwargs...) = new(kwargs)
 end
 
@@ -40,6 +40,6 @@ The h-adaptive cubature rule implemented by HCubature.jl. All standard
 `HCubature.hcubature` keyword arguments are supported.
 """
 struct HAdaptiveCubature <: IntegrationRule
-    kwargs
+    kwargs::Any
     HAdaptiveCubature(; kwargs...) = new(kwargs)
 end
diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl
index cc92adf4..426b5f50 100644
--- a/src/specializations/BezierCurve.jl
+++ b/src/specializations/BezierCurve.jl
@@ -3,22 +3,22 @@
 ################################################################################
 
 function lineintegral(
-    f::F,
-    curve::Meshes.BezierCurve,
-    rule::I,
-    FP::Type{T};
-    alg::Meshes.BezierEvalMethod=Meshes.Horner()
-) where {F<:Function, I<:IntegrationRule, T<:AbstractFloat}
-    return integral(f, curve, rule, FP; alg=alg)
+        f::F,
+        curve::Meshes.BezierCurve,
+        rule::I,
+        FP::Type{T};
+        alg::Meshes.BezierEvalMethod = Meshes.Horner()
+) where {F <: Function, I <: IntegrationRule, T <: AbstractFloat}
+    return integral(f, curve, rule, FP; alg = alg)
 end
 
 function lineintegral(
-    f::F,
-    curve::Meshes.BezierCurve,
-    rule::I;
-    alg::Meshes.BezierEvalMethod=Meshes.Horner()
-) where {F<:Function, I<:IntegrationRule}
-    return integral(f, curve, rule; alg=alg)
+        f::F,
+        curve::Meshes.BezierCurve,
+        rule::I;
+        alg::Meshes.BezierEvalMethod = Meshes.Horner()
+) where {F <: Function, I <: IntegrationRule}
+    return integral(f, curve, rule; alg = alg)
 end
 
 """
@@ -32,21 +32,21 @@ can be obtained by specifying the use of DeCasteljau's algorithm instead with
 especially for curves with a large number of control points.
 """
 function integral(
-    f::F,
-    curve::Meshes.BezierCurve,
-    rule::GaussLegendre,
-    FP::Type{T} = Float64;
-    alg::Meshes.BezierEvalMethod=Meshes.Horner()
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        curve::Meshes.BezierCurve,
+        rule::GaussLegendre,
+        FP::Type{T} = Float64;
+        alg::Meshes.BezierEvalMethod = Meshes.Horner()
+) where {F <: Function, T <: AbstractFloat}
     # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
     xs, ws = _gausslegendre(FP, rule.n)
 
     # Change of variables: x [-1,1] ↦ t [0,1]
-    t(x) = FP(1/2) * x + FP(1/2)
+    t(x) = FP(1 / 2) * x + FP(1 / 2)
     point(x) = curve(t(x), alg)
 
     # Integrate f along the line and apply a domain-correction factor for [-1,1] ↦ [0, length]
-    return FP(1/2) * length(curve) * sum(w .* f(point(x)) for (w,x) in zip(ws, xs))
+    return FP(1 / 2) * length(curve) * sum(w .* f(point(x)) for (w, x) in zip(ws, xs))
 end
 
 """
@@ -60,12 +60,12 @@ can be obtained by specifying the use of DeCasteljau's algorithm instead with
 especially for curves with a large number of control points.
 """
 function integral(
-    f::F,
-    curve::Meshes.BezierCurve,
-    rule::GaussKronrod,
-    FP::Type{T} = Float64;
-    alg::Meshes.BezierEvalMethod=Meshes.Horner()
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        curve::Meshes.BezierCurve,
+        rule::GaussKronrod,
+        FP::Type{T} = Float64;
+        alg::Meshes.BezierEvalMethod = Meshes.Horner()
+) where {F <: Function, T <: AbstractFloat}
     len = length(curve)
     point(t) = curve(t, alg)
     return QuadGK.quadgk(t -> len * f(point(t)), FP(0), FP(1); rule.kwargs...)[1]
@@ -82,19 +82,19 @@ can be obtained by specifying the use of DeCasteljau's algorithm instead with
 especially for curves with a large number of control points.
 """
 function integral(
-    f::F,
-    curve::Meshes.BezierCurve,
-    rule::HAdaptiveCubature,
-    FP::Type{T} = Float64;
-    alg::Meshes.BezierEvalMethod=Meshes.Horner()
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        curve::Meshes.BezierCurve,
+        rule::HAdaptiveCubature,
+        FP::Type{T} = Float64;
+        alg::Meshes.BezierEvalMethod = Meshes.Horner()
+) where {F <: Function, T <: AbstractFloat}
     len = length(curve)
     point(t) = curve(t, alg)
     integrand(t) = len * f(point(t[1]))
 
     # HCubature doesn't support functions that output Unitful Quantity types
     # Establish the units that are output by f
-    testpoint_parametriccoord = fill(FP(0.5),3)
+    testpoint_parametriccoord = fill(FP(0.5), 3)
     integrandunits = Unitful.unit.(integrand(testpoint_parametriccoord))
     # Create a wrapper that returns only the value component in those units
     uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
diff --git a/src/specializations/ConeSurface.jl b/src/specializations/ConeSurface.jl
index ad0e3524..78cbe639 100644
--- a/src/specializations/ConeSurface.jl
+++ b/src/specializations/ConeSurface.jl
@@ -3,11 +3,11 @@
 ################################################################################
 
 function integral(
-    f::F,
-    cone::Meshes.ConeSurface,
-    rule::I,
-    FP::Type{T} = Float64
-) where {F<:Function, I<:IntegrationRule, T<:AbstractFloat}
+        f::F,
+        cone::Meshes.ConeSurface,
+        rule::I,
+        FP::Type{T} = Float64
+) where {F <: Function, I <: IntegrationRule, T <: AbstractFloat}
     # The generic method only parameterizes the sides
     sides = _integral(f, cone, rule, FP)
 
diff --git a/src/specializations/CylinderSurface.jl b/src/specializations/CylinderSurface.jl
index f0cbe241..859b4173 100644
--- a/src/specializations/CylinderSurface.jl
+++ b/src/specializations/CylinderSurface.jl
@@ -3,11 +3,11 @@
 ################################################################################
 
 function integral(
-    f::F,
-    cyl::Meshes.CylinderSurface,
-    rule::I,
-    FP::Type{T} = Float64
-) where {F<:Function, I<:IntegrationRule, T<:AbstractFloat}
+        f::F,
+        cyl::Meshes.CylinderSurface,
+        rule::I,
+        FP::Type{T} = Float64
+) where {F <: Function, I <: IntegrationRule, T <: AbstractFloat}
     # The generic method only parameterizes the sides
     sides = _integral(f, cyl, rule, FP)
 
diff --git a/src/specializations/FrustumSurface.jl b/src/specializations/FrustumSurface.jl
index cd6427f4..101b9fce 100644
--- a/src/specializations/FrustumSurface.jl
+++ b/src/specializations/FrustumSurface.jl
@@ -3,11 +3,11 @@
 ################################################################################
 
 function integral(
-    f::F,
-    frust::Meshes.FrustumSurface,
-    rule::I,
-    FP::Type{T} = Float64
-) where  {F<:Function, I<:IntegrationRule, T<:AbstractFloat}
+        f::F,
+        frust::Meshes.FrustumSurface,
+        rule::I,
+        FP::Type{T} = Float64
+) where {F <: Function, I <: IntegrationRule, T <: AbstractFloat}
     # The generic method only parameterizes the sides
     sides = _integral(f, frust, rule, FP)
 
diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl
index ad37b341..b22f8fa5 100644
--- a/src/specializations/Line.jl
+++ b/src/specializations/Line.jl
@@ -3,11 +3,11 @@
 ################################################################################
 
 function integral(
-    f::F,
-    line::Meshes.Line,
-    rule::GaussLegendre,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        line::Meshes.Line,
+        rule::GaussLegendre,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
     # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
     xs, ws = _gausslegendre(FP, rule.n)
 
@@ -21,15 +21,15 @@ function integral(
     # Integrate f along the Line
     domainunits = _units(line(0))
     integrand(x) = f(line(t(x))) * t′(x) * domainunits
-    return sum(w .* integrand(x) for (w,x) in zip(ws, xs))
+    return sum(w .* integrand(x) for (w, x) in zip(ws, xs))
 end
 
 function integral(
-    f::F,
-    line::Meshes.Line,
-    rule::GaussKronrod,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        line::Meshes.Line,
+        rule::GaussKronrod,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
     # Normalize the Line s.t. line(t) is distance t from origin point
     line = Line(line.a, line.a + Meshes.unormalize(line.b - line.a))
 
@@ -40,11 +40,11 @@ function integral(
 end
 
 function integral(
-    f::F,
-    line::Meshes.Line,
-    rule::HAdaptiveCubature,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        line::Meshes.Line,
+        rule::HAdaptiveCubature,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
     # Normalize the Line s.t. line(t) is distance t from origin point
     line = Line(line.a, line.a + Meshes.unormalize(line.b - line.a))
 
diff --git a/src/specializations/Plane.jl b/src/specializations/Plane.jl
index fe5d90a2..d78f1159 100644
--- a/src/specializations/Plane.jl
+++ b/src/specializations/Plane.jl
@@ -3,11 +3,11 @@
 ################################################################################
 
 function integral(
-    f::F,
-    plane::Meshes.Plane,
-    rule::GaussLegendre,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        plane::Meshes.Plane,
+        rule::GaussLegendre,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
     # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]²
     xs, ws = _gausslegendre(FP, rule.n)
     wws = Iterators.product(ws, ws)
@@ -21,32 +21,36 @@ function integral(
     t′(x) = (1 + x^2) / (1 - x^2)^2
 
     # Integrate f over the Plane
-    domainunits = _units(plane(0,0))
-    integrand(((wi,wj), (xi,xj))) = wi * wj * f(plane(t(xi), t(xj))) * t′(xi) * t′(xj) * domainunits^2
-    return sum(integrand, zip(wws,xxs))
+    domainunits = _units(plane(0, 0))
+    function integrand(((wi, wj), (xi, xj)))
+        wi * wj * f(plane(t(xi), t(xj))) * t′(xi) * t′(xj) * domainunits^2
+    end
+    return sum(integrand, zip(wws, xxs))
 end
 
 function integral(
-    f::F,
-    plane::Meshes.Plane,
-    rule::GaussKronrod,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        plane::Meshes.Plane,
+        rule::GaussKronrod,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
     # Normalize the Plane's orthogonal vectors
     plane = Plane(plane.p, Meshes.unormalize(plane.u), Meshes.unormalize(plane.v))
 
     # Integrate f over the Plane
-    domainunits = _units(plane(0,0))
-    inner∫(v) = QuadGK.quadgk(u -> f(plane(u,v)) * domainunits, FP(-Inf), FP(Inf); rule.kwargs...)[1]
+    domainunits = _units(plane(0, 0))
+    function inner∫(v)
+        QuadGK.quadgk(u -> f(plane(u, v)) * domainunits, FP(-Inf), FP(Inf); rule.kwargs...)[1]
+    end
     return QuadGK.quadgk(v -> inner∫(v) * domainunits, FP(-Inf), FP(Inf); rule.kwargs...)[1]
 end
 
 function integral(
-    f::F,
-    plane::Meshes.Plane,
-    rule::HAdaptiveCubature,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        plane::Meshes.Plane,
+        rule::HAdaptiveCubature,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
     # Normalize the Plane's orthogonal vectors
     plane = Plane(plane.p, Meshes.unormalize(plane.u), Meshes.unormalize(plane.v))
 
@@ -55,8 +59,10 @@ function integral(
     t′(x) = (1 + x^2) / (1 - x^2)^2
 
     # Integrate f over the Plane
-    domainunits = _units(plane(0,0))
-    integrand(x::AbstractVector) = f(plane(t(x[1]), t(x[2]))) * t′(x[1]) * t′(x[2]) * domainunits^2
+    domainunits = _units(plane(0, 0))
+    function integrand(x::AbstractVector)
+        f(plane(t(x[1]), t(x[2]))) * t′(x[1]) * t′(x[2]) * domainunits^2
+    end
 
     # HCubature doesn't support functions that output Unitful Quantity types
     # Establish the units that are output by f
@@ -65,7 +71,7 @@ function integral(
     # Create a wrapper that returns only the value component in those units
     uintegrand(uv) = Unitful.ustrip.(integrandunits, integrand(uv))
     # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, FP[-1,-1], FP[1,1]; rule.kwargs...)[1]
+    value = HCubature.hcubature(uintegrand, FP[-1, -1], FP[1, 1]; rule.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits
diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl
index bb4eae90..dee13e0e 100644
--- a/src/specializations/Ray.jl
+++ b/src/specializations/Ray.jl
@@ -3,11 +3,11 @@
 ################################################################################
 
 function integral(
-    f::F,
-    ray::Meshes.Ray,
-    rule::GaussLegendre,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        ray::Meshes.Ray,
+        rule::GaussLegendre,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
     # Compute Gauss-Legendre nodes/weights for x in interval [-1,1]
     xs, ws = _gausslegendre(FP, rule.n)
 
@@ -15,8 +15,8 @@ function integral(
     ray = Ray(ray.p, Meshes.unormalize(ray.v))
 
     # Domain transformation: x ∈ [-1,1] ↦ t ∈ [0,∞)
-    t₁(x) = FP(1/2) * x + FP(1/2)
-    t₁′(x) = FP(1/2)
+    t₁(x) = FP(1 / 2) * x + FP(1 / 2)
+    t₁′(x) = FP(1 / 2)
     t₂(x) = x / (1 - x^2)
     t₂′(x) = (1 + x^2) / (1 - x^2)^2
     t = t₂ ∘ t₁
@@ -25,15 +25,15 @@ function integral(
     # Integrate f along the Ray
     domainunits = _units(ray(0))
     integrand(x) = f(ray(t(x))) * t′(x) * domainunits
-    return sum(w .* integrand(x) for (w,x) in zip(ws, xs))
+    return sum(w .* integrand(x) for (w, x) in zip(ws, xs))
 end
 
 function integral(
-    f::F,
-    ray::Meshes.Ray,
-    rule::GaussKronrod,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        ray::Meshes.Ray,
+        rule::GaussKronrod,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
     # Normalize the Ray s.t. ray(t) is distance t from origin point
     ray = Ray(ray.p, Meshes.unormalize(ray.v))
 
@@ -43,18 +43,18 @@ function integral(
 end
 
 function integral(
-    f::F,
-    ray::Meshes.Ray,
-    rule::HAdaptiveCubature,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        ray::Meshes.Ray,
+        rule::HAdaptiveCubature,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
     # Normalize the Ray s.t. ray(t) is distance t from origin point
     ray = Ray(ray.p, Meshes.unormalize(ray.v))
 
     # Domain transformation: x ∈ [0,1] ↦ t ∈ [0,∞)
     t(x) = x / (1 - x^2)
     t′(x) = (1 + x^2) / (1 - x^2)^2
-    
+
     # Integrate f along the Ray
     domainunits = _units(ray(0))
     integrand(x::AbstractVector) = f(ray(t(x[1]))) * t′(x[1]) * domainunits
diff --git a/src/specializations/Ring.jl b/src/specializations/Ring.jl
index 6adeb3ed..87c08407 100644
--- a/src/specializations/Ring.jl
+++ b/src/specializations/Ring.jl
@@ -3,39 +3,39 @@
 ################################################################################
 
 function integral(
-    f::F,
-    ring::Meshes.Ring,
-    rule::HAdaptiveCubature
-) where {F<:Function}
+        f::F,
+        ring::Meshes.Ring,
+        rule::HAdaptiveCubature
+) where {F <: Function}
     # Convert the Ring into Segments, sum the integrals of those 
     return sum(segment -> _integral(f, segment, rule), segments(ring))
 end
 
 function integral(
-    f::F,
-    ring::Meshes.Ring,
-    rule::HAdaptiveCubature,
-    FP::Type{T}
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        ring::Meshes.Ring,
+        rule::HAdaptiveCubature,
+        FP::Type{T}
+) where {F <: Function, T <: AbstractFloat}
     # Convert the Ring into Segments, sum the integrals of those 
     return sum(segment -> _integral(f, segment, rule, FP), segments(ring))
 end
 
 function integral(
-    f::F,
-    ring::Meshes.Ring,
-    rule::I
-) where {F<:Function, I<:IntegrationRule}
+        f::F,
+        ring::Meshes.Ring,
+        rule::I
+) where {F <: Function, I <: IntegrationRule}
     # Convert the Ring into Segments, sum the integrals of those 
     return sum(segment -> integral(f, segment, rule), segments(ring))
 end
 
 function integral(
-    f::F,
-    ring::Meshes.Ring,
-    rule::I,
-    FP::Type{T}
-) where {F<:Function, I<:IntegrationRule, T<:AbstractFloat}
+        f::F,
+        ring::Meshes.Ring,
+        rule::I,
+        FP::Type{T}
+) where {F <: Function, I <: IntegrationRule, T <: AbstractFloat}
     # Convert the Ring into Segments, sum the integrals of those 
     return sum(segment -> integral(f, segment, rule, FP), segments(ring))
 end
diff --git a/src/specializations/Rope.jl b/src/specializations/Rope.jl
index 50882364..d1fe218e 100644
--- a/src/specializations/Rope.jl
+++ b/src/specializations/Rope.jl
@@ -3,39 +3,39 @@
 ################################################################################
 
 function integral(
-    f::F,
-    rope::Meshes.Rope,
-    rule::HAdaptiveCubature
-) where {F<:Function}
+        f::F,
+        rope::Meshes.Rope,
+        rule::HAdaptiveCubature
+) where {F <: Function}
     # Convert the Rope into Segments, sum the integrals of those 
     return sum(segment -> _integral(f, segment, rule), segments(rope))
 end
 
 function integral(
-    f::F,
-    rope::Meshes.Rope,
-    rule::HAdaptiveCubature,
-    FP::Type{T}
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        rope::Meshes.Rope,
+        rule::HAdaptiveCubature,
+        FP::Type{T}
+) where {F <: Function, T <: AbstractFloat}
     # Convert the Rope into Segments, sum the integrals of those 
     return sum(segment -> _integral(f, segment, rule, FP), segments(rope))
 end
 
 function integral(
-    f::F,
-    rope::Meshes.Rope,
-    rule::I
-) where {F<:Function, I<:IntegrationRule}
+        f::F,
+        rope::Meshes.Rope,
+        rule::I
+) where {F <: Function, I <: IntegrationRule}
     # Convert the Rope into Segments, sum the integrals of those 
     return sum(segment -> integral(f, segment, rule), segments(rope))
 end
 
 function integral(
-    f::F,
-    rope::Meshes.Rope,
-    rule::I,
-    FP::Type{T}
-) where {F<:Function, I<:IntegrationRule, T<:AbstractFloat}
+        f::F,
+        rope::Meshes.Rope,
+        rule::I,
+        FP::Type{T}
+) where {F <: Function, I <: IntegrationRule, T <: AbstractFloat}
     # Convert the Rope into Segments, sum the integrals of those 
     return sum(segment -> integral(f, segment, rule, FP), segments(rope))
 end
diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl
index f5953035..9f74269e 100644
--- a/src/specializations/Tetrahedron.jl
+++ b/src/specializations/Tetrahedron.jl
@@ -3,22 +3,24 @@
 ################################################################################
 
 function integral(
-    f::F,
-    tetrahedron::Meshes.Tetrahedron,
-    rule::GaussLegendre,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        tetrahedron::Meshes.Tetrahedron,
+        rule::GaussLegendre,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
     error("Integrating a Tetrahedron with GaussLegendre not supported.")
 end
 
 function integral(
-    f::F,
-    tetrahedron::Meshes.Tetrahedron,
-    rule::GaussKronrod,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
-    inner∫₂(v,w) = QuadGK.quadgk(u -> f(tetrahedron(u,v,w)), FP(0), FP(1-v-w); rule.kwargs...)[1]
-    inner∫₁(w) = QuadGK.quadgk(v -> inner∫₂(v,w), FP(0), FP(1-w); rule.kwargs...)[1]
+        f::F,
+        tetrahedron::Meshes.Tetrahedron,
+        rule::GaussKronrod,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
+    function inner∫₂(v, w)
+        QuadGK.quadgk(u -> f(tetrahedron(u, v, w)), FP(0), FP(1 - v - w); rule.kwargs...)[1]
+    end
+    inner∫₁(w) = QuadGK.quadgk(v -> inner∫₂(v, w), FP(0), FP(1 - w); rule.kwargs...)[1]
     outer∫ = QuadGK.quadgk(w -> inner∫₁(w), FP(0), FP(1); rule.kwargs...)[1]
 
     # Apply barycentric domain correction (volume: 1/6 → actual)
@@ -26,10 +28,10 @@ function integral(
 end
 
 function integral(
-    f::F,
-    tetrahedron::Meshes.Tetrahedron,
-    rule::HAdaptiveCubature,
-    FP::Type{T} = Float64,
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        tetrahedron::Meshes.Tetrahedron,
+        rule::HAdaptiveCubature,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
     error("Integrating a Tetrahedron with HAdaptiveCubature not supported.")
 end
diff --git a/src/specializations/Triangle.jl b/src/specializations/Triangle.jl
index f92a2cdb..2406cc5e 100644
--- a/src/specializations/Triangle.jl
+++ b/src/specializations/Triangle.jl
@@ -11,11 +11,11 @@ using a Gauss-Legendre quadrature rule along each barycentric dimension of the
 triangle.
 """
 function integral(
-    f::F,
-    triangle::Meshes.Ngon{3},
-    rule::GaussLegendre,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        triangle::Meshes.Ngon{3},
+        rule::GaussLegendre,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
     # Get Gauss-Legendre nodes and weights for a 2D region [-1,1]^2
     xs, ws = _gausslegendre(FP, rule.n)
     wws = Iterators.product(ws, ws)
@@ -24,18 +24,18 @@ function integral(
     # Domain transformations:
     #   xᵢ [-1,1] ↦ R [0,1]
     #   xⱼ [-1,1] ↦ φ [0,π/2]
-    uR(xᵢ) = T(1/2) * (xᵢ + 1)
-    uφ(xⱼ) = T(π/4) * (xⱼ + 1)
+    uR(xᵢ) = T(1 / 2) * (xᵢ + 1)
+    uφ(xⱼ) = T(π / 4) * (xⱼ + 1)
 
     # Integrate the Barycentric triangle by transforming it into polar coordinates
     #   with a modified radius
     #     R = r ( sinφ + cosφ )
     #   s.t. integration bounds become rectangular
     #     R ∈ [0, 1] and φ ∈ [0, π/2]
-    function integrand(((wᵢ,wⱼ), (xᵢ,xⱼ)))
+    function integrand(((wᵢ, wⱼ), (xᵢ, xⱼ)))
         R = uR(xᵢ)
         φ = uφ(xⱼ)
-        a,b = sincos(φ)
+        a, b = sincos(φ)
         u = R * (1 - a / (a + b))
         v = R * (1 - b / (a + b))
         return wᵢ * wⱼ * f(triangle(u, v)) * R / (a + b)^2
@@ -43,7 +43,7 @@ function integral(
 
     # Calculate 2D Gauss-Legendre integral over modified-polar-Barycentric coordinates
     # Apply a linear domain-correction factor
-    return FP(π/4) * area(triangle) .* sum(integrand, zip(wws,xxs))
+    return FP(π / 4) * area(triangle) .* sum(integrand, zip(wws, xxs))
 end
 
 """
@@ -53,14 +53,16 @@ Like [`integral`](@ref) but integrates over the surface of a `triangle` using ne
 Gauss-Kronrod quadrature rules along each barycentric dimension of the triangle.
 """
 function integral(
-    f::F,
-    triangle::Meshes.Ngon{3},
-    rule::GaussKronrod,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        triangle::Meshes.Ngon{3},
+        rule::GaussKronrod,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
     # Integrate the Barycentric triangle in (u,v)-space: (0,0), (0,1), (1,0)
     #   i.e. \int_{0}^{1} \int_{0}^{1-u} f(u,v) dv du
-    inner∫(u) = QuadGK.quadgk(v -> f(triangle(u,v)), zero(FP), FP(1-u); rule.kwargs...)[1]
+    function inner∫(u)
+        QuadGK.quadgk(v -> f(triangle(u, v)), zero(FP), FP(1 - u); rule.kwargs...)[1]
+    end
     outer∫ = QuadGK.quadgk(inner∫, zero(FP), one(FP); rule.kwargs...)[1]
 
     # Apply a linear domain-correction factor 0.5 ↦ area(triangle)
@@ -75,24 +77,24 @@ transforming the triangle into a polar-barycentric coordinate system and using
 an h-adaptive cubature rule.
 """
 function integral(
-    f::F,
-    triangle::Meshes.Ngon{3},
-    rule::HAdaptiveCubature,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        triangle::Meshes.Ngon{3},
+        rule::HAdaptiveCubature,
+        FP::Type{T} = Float64
+) where {F <: Function, T <: AbstractFloat}
     # Integrate the Barycentric triangle by transforming it into polar coordinates
     #   with a modified radius
     #     R = r ( sinφ + cosφ )
     #   s.t. integration bounds become rectangular
     #     R ∈ [0, 1] and φ ∈ [0, π/2]
     function integrand(Rφ)
-        R,φ = Rφ
-        a,b = sincos(φ)
+        R, φ = Rφ
+        a, b = sincos(φ)
         u = R * (1 - a / (a + b))
         v = R * (1 - b / (a + b))
         return f(triangle(u, v)) * R / (a + b)^2
     end
-    intval = HCubature.hcubature(integrand, FP[0, 0], FP[1, π/2], rule.kwargs...)[1]
+    intval = HCubature.hcubature(integrand, FP[0, 0], FP[1, π / 2], rule.kwargs...)[1]
 
     # Apply a linear domain-correction factor 0.5 ↦ area(triangle)
     return 2 * area(triangle) .* intval
diff --git a/src/utils.jl b/src/utils.jl
index 89259b35..0d4c5cac 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -9,9 +9,9 @@ Calculate the Jacobian of a geometry at some parametric point `ts` using a simpl
 central-finite-difference approximation with step size `ε`.
 """
 function jacobian(
-    geometry,
-    ts;
-    ε=1e-6
+        geometry,
+        ts;
+        ε = 1e-6
 )
     T = eltype(ts)
 
@@ -54,7 +54,7 @@ function jacobian(
     # Allocate a re-usable ε vector
     εv = zeros(T, length(ts))
 
-    ∂ₙr(n) = ∂ₙr!(εv,ts,n)
+    ∂ₙr(n) = ∂ₙr!(εv, ts, n)
     return map(∂ₙr, 1:length(ts))
 end
 
@@ -65,8 +65,8 @@ Calculate the differential element (length, area, volume, etc) of the parametric
 function for `geometry` at arguments `ts`.
 """
 function differential(
-    geometry,
-    ts
+        geometry,
+        ts
 )
     J = jacobian(geometry, ts)
 
@@ -89,8 +89,8 @@ Determine the vector derivative of a Bezier curve `b` for the point on the
 curve parameterized by value `t`.
 """
 function derivative(
-    bz::Meshes.BezierCurve,
-    t
+        bz::Meshes.BezierCurve,
+        t
 )
     # Parameter t restricted to domain [0,1] by definition
     if t < 0 || t > 1
@@ -107,12 +107,12 @@ function derivative(
     N <= 1028 || error("This algorithm overflows for curves with ⪆1000 control points.")
 
     # Generator for Bernstein polynomial functions
-    B(i,n) = t -> binomial(Int128(n),i) * t^i * (1-t)^(n-i)
+    B(i, n) = t -> binomial(Int128(n), i) * t^i * (1 - t)^(n - i)
 
     # Derivative = N Σ_{i=0}^{N-1} sigma(i)
     #   P indices adjusted for Julia 1-based array indexing
-    sigma(i) = B(i,N-1)(t) .* (P[(i+1)+1] - P[(i)+1])
-    return N .* sum(sigma, 0:N-1)
+    sigma(i) = B(i, N - 1)(t) .* (P[(i + 1) + 1] - P[(i) + 1])
+    return N .* sum(sigma, 0:(N - 1))
 end
 
 """
@@ -122,8 +122,8 @@ Determine a unit vector pointing in the forward (t+) direction of a Bezier
 curve `b` for a point on the curve parameterized by value `t`.
 """
 function unitdirection(
-    bz::Meshes.BezierCurve, 
-    t
+        bz::Meshes.BezierCurve,
+        t
 )
     # Parameter t restricted to domain [0,1] by definition
     if t < 0 || t > 1
@@ -131,7 +131,7 @@ function unitdirection(
     end
 
     # Normalize the derivative of the curve
-    u = derivative(bz,t)
+    u = derivative(bz, t)
     LinearAlgebra.normalize!(u)
     return u
 end
@@ -147,4 +147,4 @@ function _gausslegendre(T, n)
 end
 
 # Extract the length units used by the CRS of a Point
-_units(pt::Meshes.Point{M,CRS}) where {M,CRS} = first(CoordRefSystems.units(CRS))
+_units(pt::Meshes.Point{M, CRS}) where {M, CRS} = first(CoordRefSystems.units(CRS))
diff --git a/test/aqua.jl b/test/aqua.jl
index a90deaf6..d472eddc 100644
--- a/test/aqua.jl
+++ b/test/aqua.jl
@@ -3,5 +3,5 @@
     # As of v0.11.4:
     # - Ambiguities check disabled since it fails due to upstream findings
     # - Verified that no ambiguities exist within MeshIntegrals.jl
-    Aqua.test_all(MeshIntegrals; ambiguities=false)
+    Aqua.test_all(MeshIntegrals; ambiguities = false)
 end
diff --git a/test/auto_tests.jl b/test/auto_tests.jl
index 27b9ce20..bcf03fed 100644
--- a/test/auto_tests.jl
+++ b/test/auto_tests.jl
@@ -3,7 +3,7 @@
 ################################################################################
 
 @testsnippet AutoTests begin
-    struct SupportItem{T,Dim,CRS,G<:Meshes.Geometry{Meshes.𝔼{Dim},CRS}}
+    struct SupportItem{T, Dim, CRS, G <: Meshes.Geometry{Meshes.𝔼{Dim}, CRS}}
         name::String
         type::Type{T}
         geometry::G
@@ -17,7 +17,8 @@
     end
 
     # Constructor to explicitly convert Ints (0,1) to Bool values
-    SupportItem(name, type, geometry, checkboxes::Vararg{I,7}) where {I<:Integer} = SupportItem(name, type, geometry, Bool.(checkboxes)...)
+    SupportItem(name, type, geometry, checkboxes::Vararg{I, 7}) where {I <: Integer} = SupportItem(
+        name, type, geometry, Bool.(checkboxes)...)
 
     # If method is supported, test it on scalar- and vector-valued functions.
     # Otherwise, test that its use throws a MethodError
@@ -59,14 +60,14 @@
         # For each enabled solver type, run the test suite
         @testset "$(item.name)" begin
             for ((method, methodsupport), (alg, algsupport)) in itemsupport
-                integraltest(method, item.geometry, alg, methodsupport && algsupport, item.type)
+                integraltest(
+                    method, item.geometry, alg, methodsupport && algsupport, item.type)
             end
         end
     end
 end
 
-
-@testitem "Integrals" setup = [Setup, AutoTests] begin
+@testitem "Integrals" setup=[Setup, AutoTests] begin
     # Spatial descriptors
     origin3d(T) = Point(T(0), T(0), T(0))
     origin2d(T) = Point(T(0), T(0))
@@ -83,7 +84,8 @@ end
     # Test Geometries
     ball2d(T) = Ball(origin2d(T), T(2.0))
     ball3d(T) = Ball(origin3d(T), T(2.0))
-    bezier(T) = BezierCurve([Point(cos(t), sin(t), 0) for t in range(T(0), T(2π), length=361)])
+    bezier(T) = BezierCurve([Point(cos(t), sin(t), 0)
+                             for t in range(T(0), T(2π), length = 361)])
     box1d(T) = Box(Point(T(-1)), Point(T(1)))
     box2d(T) = Box(Point(T(-1), T(-1)), Point(T(1), T(1)))
     box3d(T) = Box(Point(T(-1), T(-1), T(-1)), Point(T(1), T(1), T(-1)))
@@ -117,7 +119,7 @@ end
         SupportItem("Sphere{3,$T}", T, sphere3d(T), 1, 0, 1, 0, 1, 1, 1),
         SupportItem("Tetrahedron", T, tetra(T), 1, 0, 0, 1, 0, 1, 0),
         SupportItem("Triangle{$T}", T, triangle(T), 1, 0, 1, 0, 1, 1, 1),
-        SupportItem("Torus{$T}", T, torus(T), 1, 0, 1, 0, 1, 1, 1),
+        SupportItem("Torus{$T}", T, torus(T), 1, 0, 1, 0, 1, 1, 1)
     ]
 
     @testset "Float64 Geometries" begin
diff --git a/test/combinations.jl b/test/combinations.jl
index 0707bac6..6aedb812 100644
--- a/test/combinations.jl
+++ b/test/combinations.jl
@@ -50,14 +50,14 @@ end
 
     # Scalar integrand
     sol = _area_cone_rightcircular(h, r)
-    @test integral(f, cone, GaussLegendre(100)) ≈ sol   rtol=1e-6
-    @test integral(f, cone, GaussKronrod()) ≈ sol       rtol=1e-6
+    @test integral(f, cone, GaussLegendre(100))≈sol rtol=1e-6
+    @test integral(f, cone, GaussKronrod())≈sol rtol=1e-6
     @test integral(f, cone, HAdaptiveCubature()) ≈ sol
 
     # Vector integrand
     vsol = fill(sol, 3)
-    @test integral(fv, cone, GaussLegendre(100)) ≈ vsol   rtol=1e-6
-    @test integral(fv, cone, GaussKronrod()) ≈ vsol       rtol=1e-6
+    @test integral(fv, cone, GaussLegendre(100))≈vsol rtol=1e-6
+    @test integral(fv, cone, GaussKronrod())≈vsol rtol=1e-6
     @test integral(fv, cone, HAdaptiveCubature()) ≈ vsol
 
     # Integral aliases
@@ -94,14 +94,14 @@ end
         area_walls = area_walls_projected - area_walls_missing
         area_walls + _area_base(r_top) + _area_base(r_bot)
     end
-    @test integral(f, frustum, GaussLegendre(100)) ≈ sol    rtol=1e-6
-    @test integral(f, frustum, GaussKronrod()) ≈ sol        rtol=1e-6
+    @test integral(f, frustum, GaussLegendre(100))≈sol rtol=1e-6
+    @test integral(f, frustum, GaussKronrod())≈sol rtol=1e-6
     @test integral(f, frustum, HAdaptiveCubature()) ≈ sol
 
     # Vector integrand
     vsol = fill(sol, 3)
-    @test integral(fv, frustum, GaussLegendre(100)) ≈ vsol    rtol=1e-6
-    @test integral(fv, frustum, GaussKronrod()) ≈ vsol        rtol=1e-6
+    @test integral(fv, frustum, GaussLegendre(100))≈vsol rtol=1e-6
+    @test integral(fv, frustum, GaussKronrod())≈vsol rtol=1e-6
     @test integral(fv, frustum, HAdaptiveCubature()) ≈ vsol
 end
 
@@ -110,7 +110,7 @@ end
     b = Point(1.0u"m", 1.0u"m", 1.0u"m")
     line = Line(a, b)
 
-    function f(p::P) where {P<:Meshes.Point}
+    function f(p::P) where {P <: Meshes.Point}
         ur = hypot(p.coords.x, p.coords.y, p.coords.z)
         r = ustrip(u"m", ur)
         exp(-r^2)
@@ -140,7 +140,7 @@ end
     v = Vec(0.0u"m", 0.0u"m", 1.0u"m")
     plane = Plane(p, v)
 
-    function f(p::P) where {P<:Meshes.Point}
+    function f(p::P) where {P <: Meshes.Point}
         ur = hypot(p.coords.x, p.coords.y, p.coords.z)
         r = ustrip(u"m", ur)
         exp(-r^2)
@@ -170,7 +170,7 @@ end
     v = Vec(1.0u"m", 1.0u"m", 1.0u"m")
     ray = Ray(a, v)
 
-    function f(p::P) where {P<:Meshes.Point}
+    function f(p::P) where {P <: Meshes.Point}
         ur = hypot(p.coords.x, p.coords.y, p.coords.z)
         r = ustrip(u"m", ur)
         exp(-r^2)
@@ -202,7 +202,7 @@ end
     pt_d = Point(1.0u"m", 1.0u"m", 1.0u"m")
     rope = Ring(pt_a, pt_b, pt_c, pt_d, pt_c, pt_b)
 
-    function f(p::P) where {P<:Meshes.Point}
+    function f(p::P) where {P <: Meshes.Point}
         x, y, z = (p.coords.x, p.coords.y, p.coords.z)
         (x + 2y + 3z) * u"A/m^2"
     end
@@ -233,7 +233,7 @@ end
     pt_d = Point(1.0u"m", 1.0u"m", 1.0u"m")
     rope = Rope(pt_a, pt_b, pt_c, pt_d)
 
-    function f(p::P) where {P<:Meshes.Point}
+    function f(p::P) where {P <: Meshes.Point}
         x, y, z = (p.coords.x, p.coords.y, p.coords.z)
         (x + 2y + 3z) * u"A/m^2"
     end
@@ -261,12 +261,13 @@ end
     # Connect a line segment from the origin to an arbitrary point on the unit sphere
     phi, theta = (7pi / 6, pi / 3)  # Arbitrary spherical angles
     pt_a = Point(0.0u"m", 0.0u"m", 0.0u"m")
-    pt_b = Point(sin(theta) * cos(phi) * u"m", sin(theta) * sin(phi) * u"m", cos(theta) * u"m")
+    pt_b = Point(
+        sin(theta) * cos(phi) * u"m", sin(theta) * sin(phi) * u"m", cos(theta) * u"m")
     segment = Segment(pt_a, pt_b)
 
     a, b = (7.1, 4.6)  # arbitrary constants > 0
 
-    function f(p::P) where {P<:Meshes.Point}
+    function f(p::P) where {P <: Meshes.Point}
         ur = hypot(p.coords.x, p.coords.y, p.coords.z)
         r = ustrip(u"m", ur)
         exp(r * log(a) + (1 - r) * log(b))
diff --git a/test/floating_point_types.jl b/test/floating_point_types.jl
index 083b8766..06174c00 100644
--- a/test/floating_point_types.jl
+++ b/test/floating_point_types.jl
@@ -9,7 +9,7 @@
     )
 end
 
-@testitem "Alternate floating types" setup = [Setup, BaseAtol] begin
+@testitem "Alternate floating types" setup=[Setup, BaseAtol] begin
     @testset "$FP" for FP in (Float32, BigFloat)
         # typeof @test's are currently broken for Float32, see GitHub Issue 74
 
@@ -21,37 +21,41 @@ end
 
         # Check HCubature integrals (same method invoked for all dimensions)
         int_HC = integral(f, box1d, HAdaptiveCubature(), FP)
-        @test int_HC ≈ one(FP) * u"m" atol = baseatol[FP] * u"m"
-        @test typeof(int_HC.val) == FP broken = (FP == Float32)
+        @test int_HC≈one(FP) * u"m" atol=baseatol[FP] * u"m"
+        @test typeof(int_HC.val)==FP broken=(FP == Float32)
 
         # Check Gauss-Legendre integral in 1D
         int_GL_1D = integral(f, box1d, GaussLegendre(100), FP)
-        @test int_GL_1D ≈ one(FP) * u"m" atol = baseatol[FP] * u"m"
-        @test typeof(int_GL_1D.val) == FP broken = (FP == Float32)
+        @test int_GL_1D≈one(FP) * u"m" atol=baseatol[FP] * u"m"
+        @test typeof(int_GL_1D.val)==FP broken=(FP == Float32)
 
         # Check Gauss-Legendre integral in 2D
         int_GL_2D = integral(f, box2d, GaussLegendre(100), FP)
-        @test int_GL_2D ≈ one(FP) * u"m^2" atol = 2baseatol[FP] * u"m^2"
-        @test typeof(int_GL_2D.val) == FP broken = (FP == Float32)
+        @test int_GL_2D≈one(FP) * u"m^2" atol=2baseatol[FP] * u"m^2"
+        @test typeof(int_GL_2D.val)==FP broken=(FP == Float32)
 
         # Check Gauss-Legendre integral in 3D
         int_GL_3D = integral(f, box3d, GaussLegendre(100), FP)
-        @test int_GL_3D ≈ one(FP) * u"m^3" atol = 3baseatol[FP] * u"m^3"
-        @test typeof(int_GL_3D.val) == FP broken = (FP == Float32)
+        @test int_GL_3D≈one(FP) * u"m^3" atol=3baseatol[FP] * u"m^3"
+        @test typeof(int_GL_3D.val)==FP broken=(FP == Float32)
     end
 end
 
-@testitem "Integral Aliases" setup = [Setup] begin
+@testitem "Integral Aliases" setup=[Setup] begin
     f = p -> one(Float32)
-    box1d = Box(Point(fill(zero(Float32) * u"m", 1)...), Point(fill(one(Float32) * u"m", 1)...))
-    box2d = Box(Point(fill(zero(Float32) * u"m", 2)...), Point(fill(one(Float32) * u"m", 2)...))
-    box3d = Box(Point(fill(zero(Float32) * u"m", 3)...), Point(fill(one(Float32) * u"m", 3)...))
-    box4d = Box(Point(fill(zero(Float32) * u"m", 4)...), Point(fill(one(Float32) * u"m", 4)...))
+    box1d = Box(
+        Point(fill(zero(Float32) * u"m", 1)...), Point(fill(one(Float32) * u"m", 1)...))
+    box2d = Box(
+        Point(fill(zero(Float32) * u"m", 2)...), Point(fill(one(Float32) * u"m", 2)...))
+    box3d = Box(
+        Point(fill(zero(Float32) * u"m", 3)...), Point(fill(one(Float32) * u"m", 3)...))
+    box4d = Box(
+        Point(fill(zero(Float32) * u"m", 4)...), Point(fill(one(Float32) * u"m", 4)...))
 
     # Check alias functions for accuracy
-    @test lineintegral(f, box1d, GaussLegendre(100), Float32) ≈ 1.0f0u"m" atol = 0.01f0u"m"
-    @test surfaceintegral(f, box2d, GaussLegendre(100), Float32) ≈ 1.0f0u"m^2" atol = 0.02f0u"m^2"
-    @test volumeintegral(f, box3d, GaussLegendre(100), Float32) ≈ 1.0f0u"m^3" atol = 0.03f0u"m^3"
+    @test lineintegral(f, box1d, GaussLegendre(100), Float32)≈1.0f0u"m" atol=0.01f0u"m"
+    @test surfaceintegral(f, box2d, GaussLegendre(100), Float32)≈1.0f0u"m^2" atol=0.02f0u"m^2"
+    @test volumeintegral(f, box3d, GaussLegendre(100), Float32)≈1.0f0u"m^3" atol=0.03f0u"m^3"
 
     # Check for unsupported use of alias functions
     @test_throws "not supported" lineintegral(f, box4d, HAdaptiveCubature(), Float32)
diff --git a/test/runtests.jl b/test/runtests.jl
index 02a1c812..9689f0b0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,7 +1,7 @@
 using TestItemRunner
 using TestItems
 
-@run_package_tests verbose=true
+@run_package_tests verbose = true
 
 @testsnippet Setup begin
     using Meshes