diff --git a/src/differentiation.jl b/src/differentiation.jl
index 95232357..4433fd49 100644
--- a/src/differentiation.jl
+++ b/src/differentiation.jl
@@ -52,7 +52,7 @@ function jacobian(
         geometry::G,
         ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}}
 ) where {G <: Geometry, T <: AbstractFloat}
-    return jacobian(geometry, ts, _default_method(G))
+    return jacobian(geometry, ts, _default_diff_method(G))
 end
 
 function jacobian(
@@ -107,7 +107,7 @@ possible and finite difference approximations otherwise.
 function differential(
         geometry::G,
         ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}},
-        diff_method::DifferentiationMethod = _default_method(G)
+        diff_method::DifferentiationMethod = _default_diff_method(G)
 ) where {G <: Geometry, T <: AbstractFloat}
     J = Iterators.map(_KVector, jacobian(geometry, ts, diff_method))
     return LinearAlgebra.norm(foldl(∧, J))
diff --git a/src/integral.jl b/src/integral.jl
index 8777b86b..1fc77d44 100644
--- a/src/integral.jl
+++ b/src/integral.jl
@@ -42,7 +42,7 @@ function _integral(
         geometry,
         rule::GaussKronrod;
         FP::Type{T} = Float64,
-        diff_method::DM = _default_method(geometry)
+        diff_method::DM = _default_diff_method(geometry)
 ) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     # Implementation depends on number of parametric dimensions over which to integrate
     N = Meshes.paramdim(geometry)
@@ -70,24 +70,28 @@ function _integral(
         geometry,
         rule::GaussLegendre;
         FP::Type{T} = Float64,
-        diff_method::DM = _default_method(geometry)
+        diff_method::DM = _default_diff_method(geometry)
 ) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     N = Meshes.paramdim(geometry)
 
-    # Get Gauss-Legendre nodes and weights for a region [-1,1]^N
-    xs, ws = _gausslegendre(FP, rule.n)
-    weights = Iterators.product(ntuple(Returns(ws), N)...)
-    nodes = Iterators.product(ntuple(Returns(xs), N)...)
+    # Get Gauss-Legendre nodes and weights of type FP for a region [-1,1]ᴺ
+    xs, ws = FastGaussQuadrature.gausslegendre(rule.n)
+    xsFP = Iterators.map(FP, xs)
+    wsFP = Iterators.map(FP, ws)
+    weight_grid = Iterators.product(ntuple(Returns(wsFP), N)...)
+    node_grid = Iterators.product(ntuple(Returns(xsFP), N)...)
 
     # Domain transformation: x [-1,1] ↦ t [0,1]
     t(x) = (1 // 2) * x + (1 // 2)
 
     function integrand((weights, nodes))
-        ts = t.(nodes)
+        # Transforms nodes/xs, store in an NTuple 
+        ts = ntuple(i -> t(nodes[i]), length(nodes))
+        # Integrand function
         prod(weights) * f(geometry(ts...)) * differential(geometry, ts, diff_method)
     end
 
-    return (1 // (2^N)) .* sum(integrand, zip(weights, nodes))
+    return (1 // (2^N)) .* sum(integrand, zip(weight_grid, node_grid))
 end
 
 # HAdaptiveCubature
@@ -96,7 +100,7 @@ function _integral(
         geometry,
         rule::HAdaptiveCubature;
         FP::Type{T} = Float64,
-        diff_method::DM = _default_method(geometry)
+        diff_method::DM = _default_diff_method(geometry)
 ) where {DM <: DifferentiationMethod, T <: AbstractFloat}
     N = Meshes.paramdim(geometry)
 
diff --git a/src/utils.jl b/src/utils.jl
index 1e5292f8..6cf5ee14 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -2,50 +2,62 @@
 #                           Misc. Internal Tools
 ################################################################################
 
-# Calculate Gauss-Legendre nodes/weights and convert to type T
-function _gausslegendre(T, n)
-    xs, ws = FastGaussQuadrature.gausslegendre(n)
-    return T.(xs), T.(ws)
-end
-
 # Common error message structure
 function _error_unsupported_combination(geometry, rule)
     msg = "Integrating a $geometry using a $rule rule not supported."
     throw(ArgumentError(msg))
 end
 
-# Return an NTuple{N, T} of zeros; same interface as Base.zeros() but faster
-_zeros(T::DataType, N::Int64) = ntuple(_ -> zero(T), N)
-_zeros(N::Int) = _zeros(Float64, N)
-
-# Return an NTuple{N, T} of ones; same interface as Base.ones() but faster
-_ones(T::DataType, N::Int64) = ntuple(_ -> one(T), N)
-_ones(N::Int) = _ones(Float64, N)
-
 ################################################################################
 #                           DifferentiationMethod
 ################################################################################
 
 # Return the default DifferentiationMethod instance for a particular geometry type
-function _default_method(
+function _default_diff_method(
         g::Type{G}
 ) where {G <: Geometry}
     return FiniteDifference()
 end
 
 # Return the default DifferentiationMethod instance for a particular geometry instance
-_default_method(g::G) where {G <: Geometry} = _default_method(G)
+_default_diff_method(g::G) where {G <: Geometry} = _default_diff_method(G)
 
 ################################################################################
-#                        CliffordNumbers and Units
+#                           Numerical Tools
 ################################################################################
 
-# Meshes.Vec -> Unitful.Quantity{CliffordNumber.KVector}
+"""
+    _KVector(v::Meshes.Vec) -> Unitful.Quantity{CliffordNumbers.KVector}
+
+Convert a `Vec` to a Unitful KVector.
+"""
 function _KVector(v::Meshes.Vec{Dim, T}) where {Dim, T}
     ucoords = Iterators.map(Unitful.ustrip, v.coords)
     return CliffordNumbers.KVector{1, VGA(Dim)}(ucoords...) * _units(v)
 end
 
-# Extract the length units used by the CRS of a Geometry
+"""
+    _units(geometry)
+
+Return the Unitful.jl units associated with a particular `geometry`.
+"""
 _units(::Geometry{M, CRS}) where {M, CRS} = Unitful.unit(CoordRefSystems.lentype(CRS))
 _units(::Meshes.Vec{Dim, T}) where {Dim, T} = Unitful.unit(T)
+
+"""
+    _zeros(T = Float64, N)
+
+Return an `NTuple{N, T}` filled with zeros. This method avoids allocating an array, which
+happens when using `Base.zeros`.
+"""
+_zeros(T::DataType, N::Int64) = ntuple(_ -> zero(T), N)
+_zeros(N::Int) = _zeros(Float64, N)
+
+"""
+    _ones(T = Float64, N)
+
+Return an `NTuple{N, T}` filled with ones. This method avoids allocating an array, which
+happens when using `Base.ones`.
+"""
+_ones(T::DataType, N::Int64) = ntuple(_ -> one(T), N)
+_ones(N::Int) = _ones(Float64, N)
diff --git a/test/combinations.jl b/test/combinations.jl
index ed32a0a7..2ffb1b9a 100644
--- a/test/combinations.jl
+++ b/test/combinations.jl
@@ -66,9 +66,7 @@ end
 end
 
 @testitem "Meshes.BezierCurve" setup=[Setup] begin
-    curve = BezierCurve(
-        [Point(t * u"m", sin(t) * u"m", 0.0u"m") for t in range(-pi, pi, length = 361)]
-    )
+    curve = BezierCurve([Point(t, sin(t), 0) for t in range(-pi, pi, length = 361)])
 
     function f(p::P) where {P <: Meshes.Point}
         ux = ustrip(p.coords.x)
diff --git a/test/utils.jl b/test/utils.jl
index 545b2540..50cb4ea1 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -20,15 +20,12 @@
 end
 
 @testitem "DifferentiationMethod" setup=[Setup] begin
-    using MeshIntegrals: _default_method
+    using MeshIntegrals: _default_diff_method
 
-    # Test geometries
+    # _default_diff_method
     sphere = Sphere(Point(0, 0, 0), 1.0)
-    triangle = Triangle(Point(0, 0, 0), Point(0, 1, 0), Point(1, 0, 0))
-
-    # _default_method
-    @test _default_method(Meshes.Sphere) isa FiniteDifference
-    @test _default_method(sphere) isa FiniteDifference
+    @test _default_diff_method(Meshes.Sphere) isa FiniteDifference
+    @test _default_diff_method(sphere) isa FiniteDifference
 
     # FiniteDifference
     @test FiniteDifference().ε ≈ 1e-6