diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
deleted file mode 100644
index e7557e9..0000000
--- a/.JuliaFormatter.toml
+++ /dev/null
@@ -1,31 +0,0 @@
-#
-# alignement options
-align_assignment = true
-align_conditional = true
-align_matrix = false
-align_pair_arrow = true
-align_struct_field = true
-always_for_in = true
-always_use_return = true
-annotate_untyped_fields_with_any = false
-conditional_to_if = true
-for_in_replacement = "in"
-format_docstrings = true
-ignore = ["docs/*"]
-import_to_using = false
-indent = 4
-indent_submodule = true
-join_lines_based_on_source = false
-long_to_short_function_def = false
-margin = 92
-normalize_line_endings = "auto"
-pipe_to_function_call = false
-remove_extra_newlines = true
-separate_kwargs_with_semicolon = true
-short_to_long_function_def = true
-surround_whereop_typeparameters = true
-trailing_comma = true
-trailing_zero = true
-whitespace_in_kwargs = true
-whitespace_ops_in_indices = false
-whitespace_typedefs = false
diff --git a/.github/workflows/Lint.yml b/.github/workflows/Lint.yml
index c1c1ef9..461ab7d 100644
--- a/.github/workflows/Lint.yml
+++ b/.github/workflows/Lint.yml
@@ -26,8 +26,8 @@ jobs:
           version: "1"
       - name: Use Julia cache
         uses: julia-actions/cache@v2
-      - name: Install JuliaFormatter.jl
-        run: julia -e 'using Pkg; pkg"add JuliaFormatter"'
+      - name: Install Runic.jl
+        run: julia -e 'using Pkg; pkg"add Runic"'
       - name: Hack for setup-python cache # https://github.com/actions/setup-python/issues/807
         run: touch requirements.txt
       - name: Setup Python
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 69e72e1..ce990fd 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -14,9 +14,6 @@ repos:
       - id: check-toml
       - id: check-yaml
       - id: end-of-file-fixer
-      - id: file-contents-sorter
-        files: .JuliaFormatter.toml
-        args: [--unique]
       - id: mixed-line-ending
         args: [--fix=lf]
       # - id: no-commit-to-branch
@@ -43,7 +40,7 @@ repos:
     rev: v1.37.1
     hooks:
       - id: yamllint
-  - repo: https://github.com/domluna/JuliaFormatter.jl
-    rev: v2.1.6
+  - repo: https://github.com/fredrikekre/runic-pre-commit
+    rev: v2.0.1
     hooks:
-      - id: julia-formatter
+      - id: runic
diff --git a/README.md b/README.md
index 20c8d7a..017ce98 100644
--- a/README.md
+++ b/README.md
@@ -5,6 +5,7 @@
 [![Build Status](https://github.com/maltezfaria/LevelSetMethods.jl/workflows/Test/badge.svg)](https://github.com/maltezfaria/LevelSetMethods.jl/actions)
 [![Test workflow status](https://github.com/maltezfaria/LevelSetMethods.jl/actions/workflows/Test.yml/badge.svg?branch=main)](https://github.com/maltezfaria/LevelSetMethods.jl/actions/workflows/Test.yml?query=branch%3Amain)
 [![Lint workflow Status](https://github.com/maltezfaria/LevelSetMethods.jl/actions/workflows/Lint.yml/badge.svg?branch=main)](https://github.com/maltezfaria/LevelSetMethods.jl/actions/workflows/Lint.yml?query=branch%3Amain)
+[![Runic Badge](https://img.shields.io/badge/Runic-formatting-brightgreen.svg)](https://github.com/fredrikekre/Runic.jl)
 [![Docs workflow Status](https://github.com/maltezfaria/LevelSetMethods.jl/actions/workflows/Docs.yml/badge.svg?branch=main)](https://github.com/maltezfaria/LevelSetMethods.jl/actions/workflows/Docs.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/maltezfaria/LevelSetMethods.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/maltezfaria/LevelSetMethods.jl)
 [![BestieTemplate](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/JuliaBesties/BestieTemplate.jl/main/docs/src/assets/badge.json)](https://github.com/JuliaBesties/BestieTemplate.jl)
diff --git a/docs/make.jl b/docs/make.jl
index eb0ebf5..f30788a 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -19,7 +19,7 @@ DocMeta.setdocmeta!(
 const page_rename = Dict("developer.md" => "Developer docs") # Without the numbers
 const numbered_pages = [
     file for file in readdir(joinpath(@__DIR__, "src")) if
-    file != "index.md" && splitext(file)[2] == ".md" && occursin(r"^\d", file)
+        file != "index.md" && splitext(file)[2] == ".md" && occursin(r"^\d", file)
 ]
 
 modules = [LevelSetMethods]
diff --git a/ext/InterpolationsExt.jl b/ext/InterpolationsExt.jl
index a3f5b86..53c863e 100644
--- a/ext/InterpolationsExt.jl
+++ b/ext/InterpolationsExt.jl
@@ -4,7 +4,7 @@ import Interpolations as Itp
 import LevelSetMethods as LSM
 
 function __init__()
-    @info "Loading Interpolations extension for LevelSetMethods.jl"
+    return @info "Loading Interpolations extension for LevelSetMethods.jl"
 end
 
 Itp.interpolate(ϕ::LSM.LevelSet) = Itp.interpolate(ϕ, Itp.BSpline(Itp.Cubic()))
diff --git a/ext/MMGSurfaceExt.jl b/ext/MMGSurfaceExt.jl
index c2b3606..d69fde9 100644
--- a/ext/MMGSurfaceExt.jl
+++ b/ext/MMGSurfaceExt.jl
@@ -5,7 +5,7 @@ import MarchingCubes
 import LevelSetMethods as LSM
 
 function __init__()
-    @info "Loading MMGSurface extension for LevelSetMethods.jl"
+    return @info "Loading MMGSurface extension for LevelSetMethods.jl"
 end
 
 """
@@ -37,13 +37,13 @@ boundary and the reconstructed ideal boundary
     Only works for 3 dimensional level-set.
 """
 function LSM.export_surface_mesh(
-    ϕ::LSM.LevelSet,
-    output::String;
-    hgrad = nothing,
-    hmin = nothing,
-    hmax = nothing,
-    hausd = nothing,
-)
+        ϕ::LSM.LevelSet,
+        output::String;
+        hgrad = nothing,
+        hmin = nothing,
+        hmax = nothing,
+        hausd = nothing,
+    )
     N = LSM.dimension(ϕ)
     if N != 3
         throw(ArgumentError("export_mesh of $N dimensional level-set not supported."))
@@ -85,7 +85,7 @@ function LSM.export_surface_mesh(
 end
 
 function _write_3D_triangular_mesh(path, vertices, triangles)
-    open(path, "w") do file
+    return open(path, "w") do file
         write(file, "MeshVersionFormatted 1\n")
         write(file, "Dimension 3\n")
 
diff --git a/ext/MMGVolumeExt.jl b/ext/MMGVolumeExt.jl
index 0f2ab86..63d8eb5 100644
--- a/ext/MMGVolumeExt.jl
+++ b/ext/MMGVolumeExt.jl
@@ -5,7 +5,7 @@ import LevelSetMethods as LSM
 # using DelimitedFiles # using or import ? # faster with writedlm ?
 
 function __init__()
-    @info "Loading MMGVolume extension for LevelSetMethods.jl"
+    return @info "Loading MMGVolume extension for LevelSetMethods.jl"
 end
 
 """
@@ -39,13 +39,13 @@ For more information, see the official [MMG documentation](http://www.mmgtools.o
     Only works for 2 and 3 dimensional level-set.
 """
 function LSM.export_volume_mesh(
-    ϕ::LSM.LevelSet,
-    output::String;
-    hgrad = nothing,
-    hmin = nothing,
-    hmax = nothing,
-    hausd = nothing,
-)
+        ϕ::LSM.LevelSet,
+        output::String;
+        hgrad = nothing,
+        hmin = nothing,
+        hmax = nothing,
+        hausd = nothing,
+    )
     N = LSM.dimension(ϕ)
     if N != 2 && N != 3
         throw(ArgumentError("export_mesh of $N dimensional level-set not supported."))
@@ -78,7 +78,7 @@ function LSM.export_volume_mesh(
                 num_triangles = 2(nx - 1) * (ny - 1)
                 write(file, "\nTriangles\n")
                 write(file, "$num_triangles\n")
-                for x_id in 1:(nx-1), y_id in 1:(ny-1)
+                for x_id in 1:(nx - 1), y_id in 1:(ny - 1)
                     c00 = (y_id - 1) * nx + x_id
                     c10 = c00 + 1
                     c01 = c00 + nx
@@ -96,7 +96,7 @@ function LSM.export_volume_mesh(
                 num_tetrahedrons = 6(nx - 1) * (ny - 1) * (nz - 1)
                 write(file, "\nTetrahedra\n")
                 write(file, "$num_tetrahedrons\n")
-                for x_id in 1:(nx-1), y_id in 1:(ny-1), z_id in 1:(nz-1)
+                for x_id in 1:(nx - 1), y_id in 1:(ny - 1), z_id in 1:(nz - 1)
                     c000 = (z_id - 1) * nx * ny + (y_id - 1) * nx + x_id
                     c100 = c000 + 1
                     c010 = c000 + nx
diff --git a/ext/MakieExt.jl b/ext/MakieExt.jl
index 4940a8c..2b5d188 100644
--- a/ext/MakieExt.jl
+++ b/ext/MakieExt.jl
@@ -4,7 +4,7 @@ using Makie
 import LevelSetMethods as LSM
 
 function __init__()
-    @info "Loading Makie extension for LevelSetMethods.jl"
+    return @info "Loading Makie extension for LevelSetMethods.jl"
 end
 
 # NOTE: Makie recipes currently can't modify the Axis (https://discourse.julialang.org/t/makie-plot-recipe-collections/86434)
@@ -37,7 +37,10 @@ end
 #     end
 # end
 
-function Makie.convert_arguments(::Union{Type{<:Contour},Type{<:Contourf}}, ϕ::LSM.LevelSet)
+function Makie.convert_arguments(
+        ::Union{Type{<:Contour}, Type{<:Contourf}, Type{<:Heatmap}},
+        ϕ::LSM.LevelSet,
+    )
     LSM.dimension(ϕ) == 2 ||
         throw(ArgumentError("Contour plot only supported for 2D level-sets."))
     return _contour_plot(ϕ)
diff --git a/src/derivatives.jl b/src/derivatives.jl
index 3125db0..dcf1718 100644
--- a/src/derivatives.jl
+++ b/src/derivatives.jl
@@ -11,7 +11,7 @@ Centered finite difference scheme for first order derivative at grid point `I`
 along dimension `dim`.
 """
 function D⁰(ϕ::CartesianMeshField, I, dim)
-    h  = meshsize(ϕ, dim)
+    h = meshsize(ϕ, dim)
     Im = _decrement_index(I, dim)
     Ip = _increment_index(I, dim)
     return (ϕ[Ip] - ϕ[Im]) / (2h)
@@ -24,14 +24,14 @@ Forward finite difference scheme for first order derivative at grid point `I`
 along dimension `dim`.
 """
 @inline function D⁺(ϕ::CartesianMeshField, I, dim)
-    h  = meshsize(ϕ, dim)
+    h = meshsize(ϕ, dim)
     Ip = _increment_index(I, dim)
     return (ϕ[Ip] - ϕ[I]) / h
 end
 
 function D⁺⁺(ϕ::CartesianMeshField, I, dim)
-    h   = meshsize(ϕ, dim)
-    Ip  = _increment_index(I, dim)
+    h = meshsize(ϕ, dim)
+    Ip = _increment_index(I, dim)
     Ipp = _increment_index(I, dim, 2)
     return (-1.5 * ϕ[I] + 2 * ϕ[Ip] - 1 / 2 * ϕ[Ipp]) / h
 end
@@ -43,7 +43,7 @@ Backward finite difference scheme for first order derivative at grid point `I`
 along dimension `dim`.
 """
 function D⁻(ϕ::CartesianMeshField, I, dim)
-    h  = meshsize(ϕ, dim)
+    h = meshsize(ϕ, dim)
     Im = _decrement_index(I, dim)
     return (ϕ[I] - ϕ[Im]) / h
 end
@@ -57,9 +57,9 @@ end
 
 function weno5⁻(ϕ::CartesianMeshField, I, dim)
     # see section 3.4 of Osher-Fedwik
-    Im  = _decrement_index(I, dim)
+    Im = _decrement_index(I, dim)
     Imm = _decrement_index(Im, dim)
-    Ip  = _increment_index(I, dim)
+    Ip = _increment_index(I, dim)
     Ipp = _increment_index(Ip, dim)
     # finite differences
     v1 = D⁻(ϕ, Imm, dim)
@@ -76,7 +76,7 @@ function weno5⁻(ϕ::CartesianMeshField, I, dim)
     S2 = (13 / 12) * (v2 - 2 * v3 + v4)^2 + (1 / 4) * (v2 - v4)^2
     S3 = (13 / 12) * (v3 - 2 * v4 + v5)^2 + (1 / 4) * (3 * v3 - 4 * v4 + v5)^2
     # fudge factor
-    ϵ = 1e-6 * max(v1^2, v2^2, v3^2, v4^2, v5^2) + 1e-99
+    ϵ = 1.0e-6 * max(v1^2, v2^2, v3^2, v4^2, v5^2) + 1.0e-99
     # weights
     α1 = 0.1 / (S1 + ϵ)^2
     α2 = 0.6 / (S2 + ϵ)^2
@@ -90,9 +90,9 @@ end
 
 function weno5⁺(ϕ::CartesianMeshField, I, dim)
     # see section 3.4 of Osher-Fedwik
-    Im  = _decrement_index(I, dim)
+    Im = _decrement_index(I, dim)
     Imm = _decrement_index(Im, dim)
-    Ip  = _increment_index(I, dim)
+    Ip = _increment_index(I, dim)
     Ipp = _increment_index(Ip, dim)
     # finite differences
     v1 = D⁺(ϕ, Ipp, dim)
@@ -109,7 +109,7 @@ function weno5⁺(ϕ::CartesianMeshField, I, dim)
     S2 = (13 / 12) * (v2 - 2 * v3 + v4)^2 + (1 / 4) * (v2 - v4)^2
     S3 = (13 / 12) * (v3 - 2 * v4 + v5)^2 + (1 / 4) * (3 * v3 - 4 * v4 + v5)^2
     # fudge factor
-    ϵ = 1e-6 * max(v1^2, v2^2, v3^2, v4^2, v5^2) + 1e-99
+    ϵ = 1.0e-6 * max(v1^2, v2^2, v3^2, v4^2, v5^2) + 1.0e-99
     # weights
     α1 = 0.1 / (S1 + ϵ)^2
     α2 = 0.6 / (S2 + ϵ)^2
@@ -128,7 +128,7 @@ Centered finite difference scheme for second order derivative at grid point `I`
 along dimension `dim`. E.g. if `dim=1`, this approximates `∂ₓₓ`.
 """
 function D2⁰(ϕ::CartesianMeshField, I, dim)
-    h  = meshsize(ϕ, dim)
+    h = meshsize(ϕ, dim)
     Im = _decrement_index(I, dim)
     Ip = _increment_index(I, dim)
     return (ϕ[Ip] - 2ϕ[I] + ϕ[Im]) / h^2
@@ -143,7 +143,7 @@ along the dimensions `dims`.
 If `dims[1] == dims[2]`, it is more efficient to call `D2⁰(ϕ,I,dims[1])`.
 """
 function D2(ϕ, I, dims)
-    h  = meshsize(ϕ)
+    h = meshsize(ϕ)
     Ip = _increment_index(I, dims[1])
     Im = _decrement_index(I, dims[1])
     return (D⁰(ϕ, Ip, dims[2]) - D⁰(ϕ, Im, dims[2])) / (2 * h[dims[1]])
@@ -156,8 +156,8 @@ Upward finite difference scheme for second order derivative at grid point `I`
 along dimension `dim`. E.g. if `dim=1`, this approximates `∂ₓₓ`.
 """
 function D2⁺⁺(ϕ::CartesianMeshField, I, dim)
-    h   = meshsize(ϕ, dim)
-    Ip  = _increment_index(I, dim, 1)
+    h = meshsize(ϕ, dim)
+    Ip = _increment_index(I, dim, 1)
     Ipp = _increment_index(I, dim, 2)
     return (ϕ[I] - 2ϕ[Ip] + ϕ[Ipp]) / h^2
 end
@@ -169,8 +169,8 @@ Backward finite difference scheme for second order derivative at grid point `I`
 along dimension `dim`. E.g. if `dim=1`, this approximates `∂ₓₓ`.
 """
 function D2⁻⁻(ϕ::CartesianMeshField, I, dim)
-    h   = meshsize(ϕ, dim)
-    Im  = _decrement_index(I, dim, 1)
+    h = meshsize(ϕ, dim)
+    Im = _decrement_index(I, dim, 1)
     Imm = _decrement_index(I, dim, 2)
     return (ϕ[Imm] - 2ϕ[Im] + ϕ[I]) / h^2
 end
diff --git a/src/levelset.jl b/src/levelset.jl
index 16b55b7..4c71771 100644
--- a/src/levelset.jl
+++ b/src/levelset.jl
@@ -3,7 +3,7 @@
 
 Alias for [`MeshField`](@ref) with `vals` as an `AbstractArray` of `Real`s.
 """
-const LevelSet{V<:AbstractArray{<:Real},M,B} = MeshField{V,M,B}
+const LevelSet{V <: AbstractArray{<:Real}, M, B} = MeshField{V, M, B}
 
 function LevelSet(f::Function, m)
     vals = map(f, m)
@@ -113,24 +113,24 @@ https://en.wikipedia.org/wiki/Mean_curvature#Implicit_form_of_mean_curvature
 function curvature(ϕ::LevelSet, I)
     N = dimension(ϕ)
     if N == 2
-        ϕx  = D⁰(ϕ, I, 1)
-        ϕy  = D⁰(ϕ, I, 2)
+        ϕx = D⁰(ϕ, I, 1)
+        ϕy = D⁰(ϕ, I, 2)
         ϕxx = D2⁰(ϕ, I, 1)
         ϕyy = D2⁰(ϕ, I, 2)
         ϕxy = D2(ϕ, I, (2, 1))
-        κ   = (ϕxx * (ϕy)^2 - 2 * ϕy * ϕx * ϕxy + ϕyy * ϕx^2) / (ϕx^2 + ϕy^2)^(3 / 2)
+        κ = (ϕxx * (ϕy)^2 - 2 * ϕy * ϕx * ϕxy + ϕyy * ϕx^2) / (ϕx^2 + ϕy^2)^(3 / 2)
         return κ
     elseif N == 3
-        ϕx  = D⁰(ϕ, I, 1)
-        ϕy  = D⁰(ϕ, I, 2)
-        ϕz  = D⁰(ϕ, I, 3)
+        ϕx = D⁰(ϕ, I, 1)
+        ϕy = D⁰(ϕ, I, 2)
+        ϕz = D⁰(ϕ, I, 3)
         ϕxx = D2⁰(ϕ, I, 1)
         ϕyy = D2⁰(ϕ, I, 2)
         ϕzz = D2⁰(ϕ, I, 3)
         ϕxy = D2(ϕ, I, (2, 1))
         ϕxz = D2(ϕ, I, (3, 1))
         ϕyz = D2(ϕ, I, (3, 2))
-        κ   = (ϕxx * ϕy^2 + ϕyy * ϕx^2 + ϕxx * ϕz^2 + ϕzz * ϕx^2 + ϕyy * ϕz^2 + ϕzz  * ϕy^2 - 2 * ϕx * ϕz * ϕxz - 2 * ϕy * ϕz * ϕyz  - 2 * ϕy * ϕx * ϕxy) / (ϕx^2 + ϕy^2 + ϕz^2)^(3 / 2)
+        κ = (ϕxx * ϕy^2 + ϕyy * ϕx^2 + ϕxx * ϕz^2 + ϕzz * ϕx^2 + ϕyy * ϕz^2 + ϕzz * ϕy^2 - 2 * ϕx * ϕz * ϕxz - 2 * ϕy * ϕz * ϕyz - 2 * ϕy * ϕx * ϕxy) / (ϕx^2 + ϕy^2 + ϕz^2)^(3 / 2)
         return κ
     else
         # generic method
diff --git a/src/levelsetequation.jl b/src/levelsetequation.jl
index ae676d2..ff9dee4 100644
--- a/src/levelsetequation.jl
+++ b/src/levelsetequation.jl
@@ -138,8 +138,8 @@ function integrate!(ls::LevelSetEquation, tf, Δt = Inf)
            the level-set equation cannot be solved back in time"
     @assert tf >= tc msg
     # append boundary conditions for integration
-    ϕ         = current_state(ls)
-    buf        = buffers(ls)
+    ϕ = current_state(ls)
+    buf = buffers(ls)
     integrator = time_integrator(ls)
     # dynamic dispatch. Should not be a problem provided enough computation is
     # done inside of the function below
@@ -229,7 +229,7 @@ function _integrate!(ϕ::LevelSet, buffers, integrator::RK3, terms, tc, tf, Δt)
 end
 
 function _compute_terms(terms, ϕ, I, t)
-    sum(terms) do term
+    return sum(terms) do term
         return _compute_term(term, ϕ, I, t)
     end
 end
diff --git a/src/levelsetterms.jl b/src/levelsetterms.jl
index 341b5c4..ff943dc 100644
--- a/src/levelsetterms.jl
+++ b/src/levelsetterms.jl
@@ -6,7 +6,7 @@ A typical term in a level-set evolution equation.
 abstract type LevelSetTerm end
 
 function compute_cfl(terms, ϕ, t)
-    minimum(terms) do term
+    return minimum(terms) do term
         return _compute_cfl(term, ϕ, t)
     end
 end
@@ -21,7 +21,7 @@ function _compute_cfl(term::LevelSetTerm, ϕ, t)
     return dt
 end
 
-struct AdvectionTerm{V,S<:SpatialScheme} <: LevelSetTerm
+struct AdvectionTerm{V, S <: SpatialScheme} <: LevelSetTerm
     velocity::V
     scheme::S
 end
@@ -50,7 +50,7 @@ Base.show(io::IO, t::AdvectionTerm) = print(io, "𝐮 ⋅ ∇ ϕ")
     end
     # for dimension dim, compute the upwind derivative and multiply by the
     # velocity
-    sum(1:N) do dim
+    return sum(1:N) do dim
         v = 𝐮[dim]
         if v > 0
             if sch === Upwind()
diff --git a/src/meshes.jl b/src/meshes.jl
index 152fc7f..4a11b41 100644
--- a/src/meshes.jl
+++ b/src/meshes.jl
@@ -3,12 +3,12 @@
 
 An abstract mesh structure in dimension `N` with primite data of type `T`.
 """
-abstract type AbstractMesh{N,T} end
+abstract type AbstractMesh{N, T} end
 
-struct CartesianGrid{N,T} <: AbstractMesh{N,T}
-    lc::SVector{N,T}
-    hc::SVector{N,T}
-    n::NTuple{N,Int}
+struct CartesianGrid{N, T} <: AbstractMesh{N, T}
+    lc::SVector{N, T}
+    hc::SVector{N, T}
+    n::NTuple{N, Int}
 end
 
 """
@@ -34,10 +34,10 @@ CartesianGrid{2, Int64}([0, 0], [1, 1], (10, 4))
 function CartesianGrid(lc, hc, n)
     length(lc) == length(hc) == length(n) ||
         throw(ArgumentError("all arguments must have the same length"))
-    N   = length(lc)
-    lc_ = SVector{N,eltype(lc)}(lc...)
-    hc_ = SVector{N,eltype(hc)}(hc...)
-    n   = ntuple(i -> Int(n[i]), N)
+    N = length(lc)
+    lc_ = SVector{N, eltype(lc)}(lc...)
+    hc_ = SVector{N, eltype(hc)}(hc...)
+    n = ntuple(i -> Int(n[i]), N)
     return CartesianGrid(promote(lc_, hc_)..., n)
 end
 
@@ -50,7 +50,7 @@ xgrid(g::CartesianGrid) = grid1d(g, 1)
 ygrid(g::CartesianGrid) = grid1d(g, 2)
 zgrid(g::CartesianGrid) = grid1d(g, 3)
 
-meshsize(g::CartesianGrid)      = (g.hc .- g.lc) ./ (g.n .- 1)
+meshsize(g::CartesianGrid) = (g.hc .- g.lc) ./ (g.n .- 1)
 meshsize(g::CartesianGrid, dim) = (g.hc[dim] - g.lc[dim]) / (g.n[dim] - 1)
 
 Base.size(g::CartesianGrid) = g.n
@@ -63,10 +63,12 @@ end
 
 Base.getindex(g::CartesianGrid, I::Int...) = g[CartesianIndex(I...)]
 
+Base.eltype(g::CartesianGrid) = typeof(g.lc)
+
 function _getindex(g::CartesianGrid, I::CartesianIndex)
     N = dimension(g)
     @assert N == length(I)
-    ntuple(N) do dim
+    return ntuple(N) do dim
         return g.lc[dim] + (I[dim] - 1) / (g.n[dim] - 1) * (g.hc[dim] - g.lc[dim])
     end |> SVector
 end
@@ -77,10 +79,10 @@ Base.eachindex(g::CartesianGrid) = CartesianIndices(g)
 
 # NOTE: remove?
 function interior_indices(g::CartesianGrid, P)
-    N  = dimension(g)
+    N = dimension(g)
     sz = size(g)
-    I  = ntuple(N) do dim
-        return (P+1):(sz[dim]-P)
+    I = ntuple(N) do dim
+        return (P + 1):(sz[dim] - P)
     end
     return CartesianIndices(I)
 end
diff --git a/src/meshfield.jl b/src/meshfield.jl
index 26c1a7a..53d8ac6 100644
--- a/src/meshfield.jl
+++ b/src/meshfield.jl
@@ -6,7 +6,7 @@ A field described by its discrete values on a mesh.
 `Base.getindex` of an `MeshField` is overloaded to handle indices that lie outside the
 `CartesianIndices` of its `MeshField` by using `bcs`.
 """
-struct MeshField{V,M,B}
+struct MeshField{V, M, B}
     vals::V
     mesh::M
     bcs::B
@@ -38,7 +38,7 @@ dimension(f::MeshField) = dimension(mesh(f))
 
 # overload base methods for convenience
 function Base.getindex(ϕ::MeshField, I::CartesianIndex)
-    if has_boundary_conditions(ϕ)
+    return if has_boundary_conditions(ϕ)
         _getindex(ϕ, I)
     else
         getindex(values(ϕ), I)
@@ -81,7 +81,7 @@ end
 
 function _wrap_index_periodic(I::CartesianIndex{N}, ax, dim) where {N}
     i = I[dim]
-    ntuple(N) do d
+    return ntuple(N) do d
         if d == dim
             if i < first(ax)
                 return (last(ax) - (first(ax) - i))
@@ -95,7 +95,7 @@ end
 
 function _wrap_index_neumann(I::CartesianIndex{N}, ax, dim) where {N}
     i = I[dim]
-    ntuple(N) do d
+    return ntuple(N) do d
         if d == dim
             if i < first(ax)
                 return (first(ax) + (first(ax) - i))
@@ -139,7 +139,7 @@ Base.eachindex(ϕ::MeshField) = eachindex(mesh(ϕ))
 
 [`MeshField`](@ref) over a [`CartesianGrid`](@ref).
 """
-const CartesianMeshField{V,M<:CartesianGrid,B} = MeshField{V,M,B}
+const CartesianMeshField{V, M <: CartesianGrid, B} = MeshField{V, M, B}
 
 # Boundary conditions
 
@@ -172,11 +172,11 @@ function _getindex(ϕ::CartesianMeshField, I::CartesianIndex{N}, ::NeumannBC, d)
 end
 
 function _getindex(
-    ϕ::CartesianMeshField,
-    I::CartesianIndex{N},
-    ::NeumannGradientBC,
-    d,
-) where {N}
+        ϕ::CartesianMeshField,
+        I::CartesianIndex{N},
+        ::NeumannGradientBC,
+        d,
+    ) where {N}
     ax = axes(ϕ)[abs(d)]
     # compute mirror index to I[d]
     # TODO
@@ -193,11 +193,11 @@ end
 
 # TODO: test this
 function _getindex(
-    ϕ::CartesianMeshField,
-    I::CartesianIndex{N},
-    bc::DirichletBC,
-    d,
-) where {N}
+        ϕ::CartesianMeshField,
+        I::CartesianIndex{N},
+        bc::DirichletBC,
+        d,
+    ) where {N}
     # Compute the closest index to I that is within the domain and return value of bc there
     Iproj = clamp.(Tuple(I), axes(ϕ)) |> CartesianIndex
     x = mesh(ϕ)(Iproj)
diff --git a/test/test-boundaryconditions.jl b/test/test-boundaryconditions.jl
index 141a985..ef7b904 100644
--- a/test/test-boundaryconditions.jl
+++ b/test/test-boundaryconditions.jl
@@ -17,16 +17,16 @@ end
 @testset "Normalize BC" begin
     bcs = LSM.PeriodicBC()
     @test LSM._normalize_bc(bcs, 2) ==
-          ((LSM.PeriodicBC(), LSM.PeriodicBC()), (LSM.PeriodicBC(), LSM.PeriodicBC()))
+        ((LSM.PeriodicBC(), LSM.PeriodicBC()), (LSM.PeriodicBC(), LSM.PeriodicBC()))
     bcs = (LSM.PeriodicBC(), LSM.PeriodicBC())
     @test LSM._normalize_bc(bcs, 2) ==
-          ((LSM.PeriodicBC(), LSM.PeriodicBC()), (LSM.PeriodicBC(), LSM.PeriodicBC()))
+        ((LSM.PeriodicBC(), LSM.PeriodicBC()), (LSM.PeriodicBC(), LSM.PeriodicBC()))
     bcs = (LSM.PeriodicBC(), LSM.NeumannBC())
     @test LSM._normalize_bc(bcs, 2) ==
-          ((LSM.PeriodicBC(), LSM.PeriodicBC()), (LSM.NeumannBC(), LSM.NeumannBC()))
+        ((LSM.PeriodicBC(), LSM.PeriodicBC()), (LSM.NeumannBC(), LSM.NeumannBC()))
     bcs = [LSM.PeriodicBC(), (LSM.DirichletBC(), LSM.NeumannBC())]
     @test LSM._normalize_bc(bcs, 2) ==
-          ((LSM.PeriodicBC(), LSM.PeriodicBC()), (LSM.DirichletBC(), LSM.NeumannBC()))
+        ((LSM.PeriodicBC(), LSM.PeriodicBC()), (LSM.DirichletBC(), LSM.NeumannBC()))
     bcs = [(LSM.PeriodicBC(), LSM.DirichletBC()), (LSM.DirichletBC(), LSM.NeumannBC())]
     @test_throws ArgumentError LSM._normalize_bc(bcs, 2)
 end
diff --git a/test/test-meshes.jl b/test/test-meshes.jl
index 27743dc..b0e65aa 100644
--- a/test/test-meshes.jl
+++ b/test/test-meshes.jl
@@ -3,8 +3,8 @@ using LevelSetMethods
 
 @testset "Basic ops" begin
     nx, ny = 100, 50
-    a, b   = (-1, 0), (1, 3)
-    grid   = CartesianGrid(a, b, (nx, ny))
+    a, b = (-1, 0), (1, 3)
+    grid = CartesianGrid(a, b, (nx, ny))
     @test size(grid) === (nx, ny)
     @test length(CartesianIndices(grid)) == nx * ny
 end