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..6785a9e7 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -2,20 +2,14 @@ using Documenter
 using MeshIntegrals
 
 makedocs(
-    sitename="MeshIntegrals.jl",
+    sitename = "MeshIntegrals.jl",
     pages = [
-        "Home" => [
-            "About" => "index.md",
-            "Support Matrix" => "supportmatrix.md",
-            "Example Usage" => "usage.md"
-        ],
-        "Derivations" => [
-            "Integrating a Triangle" => "triangle.md"
-        ],
+        "Home" => ["About" => "index.md", "Support Matrix" => "supportmatrix.md",
+            "Example Usage" => "usage.md"],
+        "Derivations" => ["Integrating a Triangle" => "triangle.md"],
         "Public API" => "api.md"
     ]
 )
 
 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 48825a71..959c69f3 100644
--- a/src/MeshIntegrals.jl
+++ b/src/MeshIntegrals.jl
@@ -1,33 +1,33 @@
 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_algorithms.jl")
-    export GaussKronrod, GaussLegendre, HAdaptiveCubature
+include("integration_algorithms.jl")
+export GaussKronrod, GaussLegendre, HAdaptiveCubature
 
-    include("integral.jl")
-    include("integral_aliases.jl")
-    export integral, lineintegral, surfaceintegral, volumeintegral
+include("integral.jl")
+include("integral_aliases.jl")
+export integral, 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 f170164f..606b7d3b 100644
--- a/src/integral.jl
+++ b/src/integral.jl
@@ -24,10 +24,7 @@ contrast, increasing `FP` to e.g. `BigFloat` will typically increase precision
 function integral end
 
 # If only f and geometry are specified, select default algorithm
-function integral(
-    f::F,
-    geometry::G
-) where {F<:Function, G<:Meshes.Geometry}
+function integral(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 +32,22 @@ end
 
 # with algorithm and T specified
 function integral(
-    f::F,
-    geometry::G,
-    settings::I,
-    FP::Type{T} = Float64
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm, T<:AbstractFloat}
+        f::F,
+        geometry::G,
+        settings::I,
+        FP::Type{T} = Float64
+) where {
+        F <: Function, G <: Meshes.Geometry, I <: IntegrationAlgorithm, T <: AbstractFloat}
     _integral(f, geometry, settings, FP)
 end
 
-
 ################################################################################
 #                    Generalized (n-Dimensional) Worker Methods
 ################################################################################
 
 # GaussKronrod
-function _integral(
-    f,
-    geometry,
-    settings::GaussKronrod,
-    FP::Type{T} = Float64
-) where {T<:AbstractFloat}
+function _integral(f, geometry, settings::GaussKronrod,
+        FP::Type{T} = Float64) where {T <: AbstractFloat}
     # Run the appropriate integral type
     Dim = Meshes.paramdim(geometry)
     if Dim == 1
@@ -67,12 +60,8 @@ function _integral(
 end
 
 # GaussLegendre
-function _integral(
-    f,
-    geometry,
-    settings::GaussLegendre,
-    FP::Type{T} = Float64
-) where {T<:AbstractFloat}
+function _integral(f, geometry, settings::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,35 +70,32 @@ 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,
-    settings::HAdaptiveCubature,
-    FP::Type{T} = Float64
-) where {T<:AbstractFloat}
+function _integral(f, geometry, settings::HAdaptiveCubature,
+        FP::Type{T} = Float64) where {T <: AbstractFloat}
     Dim = Meshes.paramdim(geometry)
 
     integrand(t) = f(geometry(t...)) * differential(geometry, t)
 
     # HCubature doesn't support functions that output Unitful Quantity types
     # Establish the units that are output by f
-    testpoint_parametriccoord = fill(FP(0.5),Dim)
+    testpoint_parametriccoord = fill(FP(0.5), Dim)
     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))
     # Integrate only the unitless values
-    value = HCubature.hcubature(uintegrand, zeros(FP,Dim), ones(FP,Dim); settings.kwargs...)[1]
+    value = HCubature.hcubature(
+        uintegrand, zeros(FP, Dim), ones(FP, Dim); settings.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits
@@ -119,33 +105,21 @@ end
 #                    Specialized GaussKronrod Methods
 ################################################################################
 
-function _integral_1d(
-    f,
-    geometry,
-    settings::GaussKronrod,
-    FP::Type{T} = Float64
-) where {T<:AbstractFloat}
+function _integral_1d(f, geometry, settings::GaussKronrod,
+        FP::Type{T} = Float64) where {T <: AbstractFloat}
     integrand(t) = f(geometry(t)) * differential(geometry, [t])
     return QuadGK.quadgk(integrand, FP(0), FP(1); settings.kwargs...)[1]
 end
 
-function _integral_2d(
-    f,
-    geometry2d,
-    settings::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), FP(0), FP(1); settings.kwargs...)[1]
+function _integral_2d(f, geometry2d, settings::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), FP(0), FP(1); settings.kwargs...)[1]
     return QuadGK.quadgk(v -> ∫₁(v), FP(0), FP(1); settings.kwargs...)[1]
 end
 
 # Integrating volumes with GaussKronrod not supported by default
-function _integral_3d(
-    f,
-    geometry,
-    settings::GaussKronrod,
-    FP::Type{T} = Float64
-) where {T<:AbstractFloat}
+function _integral_3d(f, geometry, settings::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 92da3349..2f0288e9 100644
--- a/src/integral_aliases.jl
+++ b/src/integral_aliases.jl
@@ -16,10 +16,7 @@ Algorithm types available:
 - GaussLegendre
 - HAdaptiveCubature
 """
-function lineintegral(
-    f::F,
-    geometry::G,
-) where {F<:Function, G<:Meshes.Geometry}
+function lineintegral(f::F, geometry::G) where {F <: Function, G <: Meshes.Geometry}
     Dim = Meshes.paramdim(geometry)
 
     if Dim == 1
@@ -29,11 +26,8 @@ function lineintegral(
     end
 end
 
-function lineintegral(
-    f::F,
-    geometry::G,
-    settings::I
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm}
+function lineintegral(f::F, geometry::G,
+        settings::I) where {F <: Function, G <: Meshes.Geometry, I <: IntegrationAlgorithm}
     Dim = Meshes.paramdim(geometry)
 
     if Dim == 1
@@ -44,11 +38,12 @@ function lineintegral(
 end
 
 function lineintegral(
-    f::F,
-    geometry::G,
-    settings::I,
-    FP::Type{T}
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm, T<:AbstractFloat}
+        f::F,
+        geometry::G,
+        settings::I,
+        FP::Type{T}
+) where {
+        F <: Function, G <: Meshes.Geometry, I <: IntegrationAlgorithm, T <: AbstractFloat}
     Dim = Meshes.paramdim(geometry)
 
     if Dim == 1
@@ -58,7 +53,6 @@ function lineintegral(
     end
 end
 
-
 ################################################################################
 #                              Surface Integral
 ################################################################################
@@ -77,10 +71,7 @@ Algorithm types available:
 - GaussLegendre
 - HAdaptiveCubature (default)
 """
-function surfaceintegral(
-    f::F,
-    geometry::G,
-) where {F<:Function, G<:Meshes.Geometry}
+function surfaceintegral(f::F, geometry::G) where {F <: Function, G <: Meshes.Geometry}
     Dim = Meshes.paramdim(geometry)
 
     if Dim == 2
@@ -90,11 +81,8 @@ function surfaceintegral(
     end
 end
 
-function surfaceintegral(
-    f::F,
-    geometry::G,
-    settings::I
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm}
+function surfaceintegral(f::F, geometry::G,
+        settings::I) where {F <: Function, G <: Meshes.Geometry, I <: IntegrationAlgorithm}
     Dim = Meshes.paramdim(geometry)
 
     if Dim == 2
@@ -105,11 +93,12 @@ function surfaceintegral(
 end
 
 function surfaceintegral(
-    f::F,
-    geometry::G,
-    settings::I,
-    FP::Type{T}
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm, T<:AbstractFloat}
+        f::F,
+        geometry::G,
+        settings::I,
+        FP::Type{T}
+) where {
+        F <: Function, G <: Meshes.Geometry, I <: IntegrationAlgorithm, T <: AbstractFloat}
     Dim = Meshes.paramdim(geometry)
 
     if Dim == 2
@@ -119,7 +108,6 @@ function surfaceintegral(
     end
 end
 
-
 ################################################################################
 #                              Volume Integral
 ################################################################################
@@ -138,10 +126,7 @@ Algorithm types available:
 - GaussLegendre
 - HAdaptiveCubature (default)
 """
-function volumeintegral(
-    f::F,
-    geometry::G,
-) where {F<:Function, G<:Meshes.Geometry}
+function volumeintegral(f::F, geometry::G) where {F <: Function, G <: Meshes.Geometry}
     Dim = Meshes.paramdim(geometry)
 
     if Dim == 3
@@ -151,11 +136,8 @@ function volumeintegral(
     end
 end
 
-function volumeintegral(
-    f::F,
-    geometry::G,
-    settings::I
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm}
+function volumeintegral(f::F, geometry::G,
+        settings::I) where {F <: Function, G <: Meshes.Geometry, I <: IntegrationAlgorithm}
     Dim = Meshes.paramdim(geometry)
 
     if Dim == 3
@@ -166,11 +148,12 @@ function volumeintegral(
 end
 
 function volumeintegral(
-    f::F,
-    geometry::G,
-    settings::I,
-    FP::Type{T}
-) where {F<:Function, G<:Meshes.Geometry, I<:IntegrationAlgorithm, T<:AbstractFloat}
+        f::F,
+        geometry::G,
+        settings::I,
+        FP::Type{T}
+) where {
+        F <: Function, G <: Meshes.Geometry, I <: IntegrationAlgorithm, T <: AbstractFloat}
     Dim = Meshes.paramdim(geometry)
 
     if Dim == 3
diff --git a/src/integration_algorithms.jl b/src/integration_algorithms.jl
index 8a8eb0b5..39569909 100644
--- a/src/integration_algorithms.jl
+++ b/src/integration_algorithms.jl
@@ -11,7 +11,7 @@ Numerically integrate using the h-adaptive Gauss-Kronrod quadrature rule impleme
 by QuadGK.jl. All standard `QuadGK.quadgk` keyword arguments are supported.
 """
 struct GaussKronrod <: IntegrationAlgorithm
-    kwargs
+    kwargs::Any
     GaussKronrod(; kwargs...) = new(kwargs)
 end
 
@@ -39,6 +39,6 @@ implemented by HCubature.jl. All standard `HCubature.hcubature` keyword
 arguments are supported.
 """
 struct HAdaptiveCubature <: IntegrationAlgorithm
-    kwargs
+    kwargs::Any
     HAdaptiveCubature(; kwargs...) = new(kwargs)
 end
diff --git a/src/specializations/BezierCurve.jl b/src/specializations/BezierCurve.jl
index f85b9116..947855a9 100644
--- a/src/specializations/BezierCurve.jl
+++ b/src/specializations/BezierCurve.jl
@@ -3,22 +3,22 @@
 ################################################################################
 
 function lineintegral(
-    f::F,
-    curve::Meshes.BezierCurve,
-    settings::I,
-    FP::Type{T};
-    alg::Meshes.BezierEvalMethod=Meshes.Horner()
-) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat}
-    return integral(f, curve, settings, FP; alg=alg)
+        f::F,
+        curve::Meshes.BezierCurve,
+        settings::I,
+        FP::Type{T};
+        alg::Meshes.BezierEvalMethod = Meshes.Horner()
+) where {F <: Function, I <: IntegrationAlgorithm, T <: AbstractFloat}
+    return integral(f, curve, settings, FP; alg = alg)
 end
 
 function lineintegral(
-    f::F,
-    curve::Meshes.BezierCurve,
-    settings::I;
-    alg::Meshes.BezierEvalMethod=Meshes.Horner()
-) where {F<:Function, I<:IntegrationAlgorithm}
-    return integral(f, curve, settings; alg=alg)
+        f::F,
+        curve::Meshes.BezierCurve,
+        settings::I;
+        alg::Meshes.BezierEvalMethod = Meshes.Horner()
+) where {F <: Function, I <: IntegrationAlgorithm}
+    return integral(f, curve, settings; 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,
-    settings::GaussLegendre,
-    FP::Type{T} = Float64;
-    alg::Meshes.BezierEvalMethod=Meshes.Horner()
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        curve::Meshes.BezierCurve,
+        settings::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, settings.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,
-    settings::GaussKronrod,
-    FP::Type{T} = Float64;
-    alg::Meshes.BezierEvalMethod=Meshes.Horner()
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        curve::Meshes.BezierCurve,
+        settings::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); settings.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,
-    settings::HAdaptiveCubature,
-    FP::Type{T} = Float64;
-    alg::Meshes.BezierEvalMethod=Meshes.Horner()
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        curve::Meshes.BezierCurve,
+        settings::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 67cbc84f..714ea894 100644
--- a/src/specializations/ConeSurface.jl
+++ b/src/specializations/ConeSurface.jl
@@ -3,11 +3,11 @@
 ################################################################################
 
 function integral(
-    f::F,
-    cone::Meshes.ConeSurface,
-    settings::I,
-    FP::Type{T} = Float64
-) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat}
+        f::F,
+        cone::Meshes.ConeSurface,
+        settings::I,
+        FP::Type{T} = Float64
+) where {F <: Function, I <: IntegrationAlgorithm, T <: AbstractFloat}
     # The generic method only parameterizes the sides
     sides = _integral(f, cone, settings, FP)
 
diff --git a/src/specializations/CylinderSurface.jl b/src/specializations/CylinderSurface.jl
index 7af8f213..3bc225da 100644
--- a/src/specializations/CylinderSurface.jl
+++ b/src/specializations/CylinderSurface.jl
@@ -3,11 +3,11 @@
 ################################################################################
 
 function integral(
-    f::F,
-    cyl::Meshes.CylinderSurface,
-    settings::I,
-    FP::Type{T} = Float64
-) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat}
+        f::F,
+        cyl::Meshes.CylinderSurface,
+        settings::I,
+        FP::Type{T} = Float64
+) where {F <: Function, I <: IntegrationAlgorithm, T <: AbstractFloat}
     # The generic method only parameterizes the sides
     sides = _integral(f, cyl, settings, FP)
 
diff --git a/src/specializations/FrustumSurface.jl b/src/specializations/FrustumSurface.jl
index 96bea658..7300e554 100644
--- a/src/specializations/FrustumSurface.jl
+++ b/src/specializations/FrustumSurface.jl
@@ -3,11 +3,11 @@
 ################################################################################
 
 function integral(
-    f::F,
-    frust::Meshes.FrustumSurface,
-    settings::I,
-    FP::Type{T} = Float64
-) where  {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat}
+        f::F,
+        frust::Meshes.FrustumSurface,
+        settings::I,
+        FP::Type{T} = Float64
+) where {F <: Function, I <: IntegrationAlgorithm, T <: AbstractFloat}
     # The generic method only parameterizes the sides
     sides = _integral(f, frust, settings, FP)
 
diff --git a/src/specializations/Line.jl b/src/specializations/Line.jl
index 0c4d431a..9aef0f07 100644
--- a/src/specializations/Line.jl
+++ b/src/specializations/Line.jl
@@ -3,11 +3,11 @@
 ################################################################################
 
 function integral(
-    f::F,
-    line::Meshes.Line,
-    settings::GaussLegendre,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        line::Meshes.Line,
+        settings::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, settings.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,
-    settings::GaussKronrod,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        line::Meshes.Line,
+        settings::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,
-    settings::HAdaptiveCubature,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        line::Meshes.Line,
+        settings::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 e074c5c4..95233f03 100644
--- a/src/specializations/Plane.jl
+++ b/src/specializations/Plane.jl
@@ -3,11 +3,11 @@
 ################################################################################
 
 function integral(
-    f::F,
-    plane::Meshes.Plane,
-    settings::GaussLegendre,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        plane::Meshes.Plane,
+        settings::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, settings.n)
     wws = Iterators.product(ws, ws)
@@ -21,32 +21,38 @@ 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,
-    settings::GaussKronrod,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        plane::Meshes.Plane,
+        settings::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); settings.kwargs...)[1]
-    return QuadGK.quadgk(v -> inner∫(v) * domainunits, FP(-Inf), FP(Inf); settings.kwargs...)[1]
+    domainunits = _units(plane(0, 0))
+    function inner∫(v)
+        QuadGK.quadgk(
+            u -> f(plane(u, v)) * domainunits, FP(-Inf), FP(Inf); settings.kwargs...)[1]
+    end
+    return QuadGK.quadgk(
+        v -> inner∫(v) * domainunits, FP(-Inf), FP(Inf); settings.kwargs...)[1]
 end
 
 function integral(
-    f::F,
-    plane::Meshes.Plane,
-    settings::HAdaptiveCubature,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        plane::Meshes.Plane,
+        settings::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 +61,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 +73,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]; settings.kwargs...)[1]
+    value = HCubature.hcubature(uintegrand, FP[-1, -1], FP[1, 1]; settings.kwargs...)[1]
 
     # Reapply units
     return value .* integrandunits
diff --git a/src/specializations/Ray.jl b/src/specializations/Ray.jl
index b7d3f58f..c0042a3b 100644
--- a/src/specializations/Ray.jl
+++ b/src/specializations/Ray.jl
@@ -3,11 +3,11 @@
 ################################################################################
 
 function integral(
-    f::F,
-    ray::Meshes.Ray,
-    settings::GaussLegendre,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        ray::Meshes.Ray,
+        settings::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, settings.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,
-    settings::GaussKronrod,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        ray::Meshes.Ray,
+        settings::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,
-    settings::HAdaptiveCubature,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        ray::Meshes.Ray,
+        settings::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 4876a768..40f67363 100644
--- a/src/specializations/Ring.jl
+++ b/src/specializations/Ring.jl
@@ -3,39 +3,33 @@
 ################################################################################
 
 function integral(
-    f::F,
-    ring::Meshes.Ring,
-    settings::HAdaptiveCubature
-) where {F<:Function}
+        f::F, ring::Meshes.Ring, settings::HAdaptiveCubature) where {F <: Function}
     # Convert the Ring into Segments, sum the integrals of those 
     return sum(segment -> _integral(f, segment, settings), segments(ring))
 end
 
 function integral(
-    f::F,
-    ring::Meshes.Ring,
-    settings::HAdaptiveCubature,
-    FP::Type{T}
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        ring::Meshes.Ring,
+        settings::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, settings, FP), segments(ring))
 end
 
-function integral(
-    f::F,
-    ring::Meshes.Ring,
-    settings::I
-) where {F<:Function, I<:IntegrationAlgorithm}
+function integral(f::F, ring::Meshes.Ring,
+        settings::I) where {F <: Function, I <: IntegrationAlgorithm}
     # Convert the Ring into Segments, sum the integrals of those 
     return sum(segment -> integral(f, segment, settings), segments(ring))
 end
 
 function integral(
-    f::F,
-    ring::Meshes.Ring,
-    settings::I,
-    FP::Type{T}
-) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat}
+        f::F,
+        ring::Meshes.Ring,
+        settings::I,
+        FP::Type{T}
+) where {F <: Function, I <: IntegrationAlgorithm, T <: AbstractFloat}
     # Convert the Ring into Segments, sum the integrals of those 
     return sum(segment -> integral(f, segment, settings, FP), segments(ring))
 end
diff --git a/src/specializations/Rope.jl b/src/specializations/Rope.jl
index 03917618..f4c0471f 100644
--- a/src/specializations/Rope.jl
+++ b/src/specializations/Rope.jl
@@ -3,39 +3,33 @@
 ################################################################################
 
 function integral(
-    f::F,
-    rope::Meshes.Rope,
-    settings::HAdaptiveCubature
-) where {F<:Function}
+        f::F, rope::Meshes.Rope, settings::HAdaptiveCubature) where {F <: Function}
     # Convert the Rope into Segments, sum the integrals of those 
     return sum(segment -> _integral(f, segment, settings), segments(rope))
 end
 
 function integral(
-    f::F,
-    rope::Meshes.Rope,
-    settings::HAdaptiveCubature,
-    FP::Type{T}
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        rope::Meshes.Rope,
+        settings::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, settings, FP), segments(rope))
 end
 
-function integral(
-    f::F,
-    rope::Meshes.Rope,
-    settings::I
-) where {F<:Function, I<:IntegrationAlgorithm}
+function integral(f::F, rope::Meshes.Rope,
+        settings::I) where {F <: Function, I <: IntegrationAlgorithm}
     # Convert the Rope into Segments, sum the integrals of those 
     return sum(segment -> integral(f, segment, settings), segments(rope))
 end
 
 function integral(
-    f::F,
-    rope::Meshes.Rope,
-    settings::I,
-    FP::Type{T}
-) where {F<:Function, I<:IntegrationAlgorithm, T<:AbstractFloat}
+        f::F,
+        rope::Meshes.Rope,
+        settings::I,
+        FP::Type{T}
+) where {F <: Function, I <: IntegrationAlgorithm, T <: AbstractFloat}
     # Convert the Rope into Segments, sum the integrals of those 
     return sum(segment -> integral(f, segment, settings, FP), segments(rope))
 end
diff --git a/src/specializations/Tetrahedron.jl b/src/specializations/Tetrahedron.jl
index 51beec43..2215819c 100644
--- a/src/specializations/Tetrahedron.jl
+++ b/src/specializations/Tetrahedron.jl
@@ -3,22 +3,25 @@
 ################################################################################
 
 function integral(
-    f::F,
-    tetrahedron::Meshes.Tetrahedron,
-    settings::GaussLegendre,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        tetrahedron::Meshes.Tetrahedron,
+        settings::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,
-    settings::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); settings.kwargs...)[1]
-    inner∫₁(w) = QuadGK.quadgk(v -> inner∫₂(v,w), FP(0), FP(1-w); settings.kwargs...)[1]
+        f::F,
+        tetrahedron::Meshes.Tetrahedron,
+        settings::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); settings.kwargs...)[1]
+    end
+    inner∫₁(w) = QuadGK.quadgk(v -> inner∫₂(v, w), FP(0), FP(1 - w); settings.kwargs...)[1]
     outer∫ = QuadGK.quadgk(w -> inner∫₁(w), FP(0), FP(1); settings.kwargs...)[1]
 
     # Apply barycentric domain correction (volume: 1/6 → actual)
@@ -26,10 +29,10 @@ function integral(
 end
 
 function integral(
-    f::F,
-    tetrahedron::Meshes.Tetrahedron,
-    settings::HAdaptiveCubature,
-    FP::Type{T} = Float64,
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        tetrahedron::Meshes.Tetrahedron,
+        settings::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 53982dc9..b66bf08a 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},
-    settings::GaussLegendre,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        triangle::Meshes.Ngon{3},
+        settings::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, settings.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},
-    settings::GaussKronrod,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        triangle::Meshes.Ngon{3},
+        settings::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)), FP(0), FP(1-u); settings.kwargs...)[1]
+    function inner∫(u)
+        QuadGK.quadgk(v -> f(triangle(u, v)), FP(0), FP(1 - u); settings.kwargs...)[1]
+    end
     outer∫ = QuadGK.quadgk(inner∫, FP(0), FP(1); settings.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},
-    settings::HAdaptiveCubature,
-    FP::Type{T} = Float64
-) where {F<:Function, T<:AbstractFloat}
+        f::F,
+        triangle::Meshes.Ngon{3},
+        settings::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], settings.kwargs...)[1]
+    intval = HCubature.hcubature(integrand, FP[0, 0], FP[1, π / 2], settings.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..14d062fe 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -8,11 +8,7 @@
 Calculate the Jacobian of a geometry at some parametric point `ts` using a simple
 central-finite-difference approximation with step size `ε`.
 """
-function jacobian(
-    geometry,
-    ts;
-    ε=1e-6
-)
+function jacobian(geometry, ts; ε = 1e-6)
     T = eltype(ts)
 
     # Get the partial derivative along the n'th axis via finite difference approximation
@@ -54,7 +50,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
 
@@ -64,10 +60,7 @@ end
 Calculate the differential element (length, area, volume, etc) of the parametric
 function for `geometry` at arguments `ts`.
 """
-function differential(
-    geometry,
-    ts
-)
+function differential(geometry, ts)
     J = jacobian(geometry, ts)
 
     # TODO generalize this with geometric algebra, e.g.: norm(foldl(∧, J))
@@ -88,10 +81,7 @@ end
 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
-)
+function derivative(bz::Meshes.BezierCurve, t)
     # Parameter t restricted to domain [0,1] by definition
     if t < 0 || t > 1
         throw(DomainError(t, "b(t) is not defined for t outside [0, 1]."))
@@ -107,12 +97,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
 
 """
@@ -121,17 +111,14 @@ end
 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
-)
+function unitdirection(bz::Meshes.BezierCurve, t)
     # Parameter t restricted to domain [0,1] by definition
     if t < 0 || t > 1
         throw(DomainError(t, "b(t) is not defined for t outside [0, 1]."))
     end
 
     # Normalize the derivative of the curve
-    u = derivative(bz,t)
+    u = derivative(bz, t)
     LinearAlgebra.normalize!(u)
     return u
 end
@@ -147,4 +134,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..d6f37c5d 100644
--- a/test/floating_point_types.jl
+++ b/test/floating_point_types.jl
@@ -3,13 +3,10 @@
 
 # Base value for atol when integrating with a particular FP type
 @testsnippet BaseAtol begin
-    baseatol = Dict(
-        Float32 => 0.01f0,
-        BigFloat => BigFloat(0.001)
-    )
+    baseatol = Dict(Float32 => 0.01f0, BigFloat => BigFloat(0.001))
 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 +18,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