diff --git a/src/AlgoimUtils/AlgoimUtils.jl b/src/AlgoimUtils/AlgoimUtils.jl
index 95339d9..8b994c7 100644
--- a/src/AlgoimUtils/AlgoimUtils.jl
+++ b/src/AlgoimUtils/AlgoimUtils.jl
@@ -516,6 +516,22 @@ function compute_normal_displacement(
   disps
 end
 
+function compute_normal_displacement(
+    cps::AbstractArray,
+    phi::AlgoimCallLevelSetFunction,
+    fun::DistributedCellField,
+    dt::Float64,
+    Ω::DistributedTriangulation)
+  # TO-DO: Optimize this function such that it reuses caches
+  map(cps) do cp
+    fun_xs = evaluate(fun,cp)
+    nΓ_xs = evaluate(normal(phi,Ω),cp)
+    lazy_map(fun_xs,nΓ_xs) do fun_x,nΓ_x
+      dt * ( fun_x ⋅ nΓ_x )
+    end
+  end
+end
+
 function compute_normal_displacement!(
     cache,
     cps::AbstractVector{<:Point},
diff --git a/test/AlgoimUtilsTests/ClosestPointTests.jl b/test/AlgoimUtilsTests/ClosestPointTests.jl
index c7042c1..8e1027e 100644
--- a/test/AlgoimUtilsTests/ClosestPointTests.jl
+++ b/test/AlgoimUtilsTests/ClosestPointTests.jl
@@ -16,9 +16,9 @@ const CUT = 0
 function run_case(distribute,parts,degree)
 
   ranks = distribute(LinearIndices((prod(parts),)))
-  cells = Int32[32,32,32]
+  cells = Int32[4,4]
   
-  domain = (-1.1,1.1,-1.1,1.1,-1.1,1.1)
+  domain = (-1.1,1.1,-1.1,1.1)
   bgmodel = CartesianDiscreteModel(ranks,parts,domain,cells)
   Ω = Triangulation(bgmodel)
   
@@ -26,35 +26,25 @@ function run_case(distribute,parts,degree)
   reffeᵠ = ReferenceFE(lagrangian,Float64,order)
   V = TestFESpace(Ω,reffeᵠ)
   
-  # f(x) = x[1]*x[1]+x[2]*x[2] - 0.25
-  f(x) = x[1]*x[1]+x[2]*x[2]+x[3]*x[3] - 0.25
+  f(x) = x[1]*x[1]+x[2]*x[2] - 0.25
+  # f(x) = x[1]*x[1]+x[2]*x[2]+x[3]*x[3] - 0.25
   ∇f(x) = ∇(f)(x)
   phi = AlgoimCallLevelSetFunction(f,∇f)
 
   fₕ = compute_distance_fe_function(bgmodel,V,phi,order,cppdegree=degree)
   writevtk(Ω,"distance",cellfields=["f"=>fₕ])
 
-  # coords = map(local_views(bgmodel)) do t
-  #   _p = reduce(vcat,collect(get_node_coordinates(t)))
-  #   [ Point(_p[i]...) for i in 1:length(_p) ]
-  # end
-  # _cpps = map(coords) do c
-  #   [ p-f(p)*∇f(p) for p in c ]
-  # end
-  # cpps = compute_closest_point_projections(bgmodel,phi,cppdegree=degree)
-  # map(coords,cpps,_cpps,ranks) do c,p,_p,i
-  #   writevtk(c,"cpps_$i",nodaldata=["f"=>f.(p),"cp1"=>p,"cp2"=>_p])
-  # end
-
-  # partials = map(cpps) do c
-  #   maximum(f.(c))
-  # end
-  # maxf = reduce(max,partials,init=typemin(Float64))
-  # partials = map(cpps) do c
-  #   minimum(f.(c))
-  # end
-  # minf = reduce(min,partials,init=typemax(Float64))
-  # @show maxf,minf
+  φₕ = AlgoimCallLevelSetFunction(fₕ,∇(fₕ))
+
+  # cp₋ = compute_closest_point_projections(V,φₕ,order)
+  cpps = compute_closest_point_projections(V,phi,order)
+  nΓ = normal(phi,Ω)
+  dt = 0.1
+
+  dₕ = compute_normal_displacement(cpps,phi,nΓ,dt,Ω)
+  map(dₕ) do d
+    @show d
+  end
 
 end