diff --git a/src/AlgoimUtils/AlgoimUtils.jl b/src/AlgoimUtils/AlgoimUtils.jl
index c5831e6..3e8f64a 100644
--- a/src/AlgoimUtils/AlgoimUtils.jl
+++ b/src/AlgoimUtils/AlgoimUtils.jl
@@ -44,6 +44,7 @@ export init_bboxes
 export fill_cpp_data
 export fill_cpp_data_raw
 export compute_closest_point_projections
+export compute_displacement
 export compute_normal_displacement
 export compute_distance_fe_function
 export delaunaytrian
@@ -290,7 +291,6 @@ function compute_closest_point_projections(model::AdaptedDiscreteModel,
                                            limitstol::Float64=1.0e-8)
   reffe = ReferenceFE(lagrangian,Float64,order)
   rfespace = TestFESpace(model,reffe)
-  _φ = φ.φ
   _rφ = interpolate_everywhere(φ.φ,rfespace)
   rφ = AlgoimCallLevelSetFunction(_rφ,∇(_rφ))
   cdesc = get_cartesian_descriptor(get_model(model))
@@ -395,9 +395,54 @@ function compute_normal_displacement(
   disps
 end
 
-_dist(x,y) = begin
-  sign = norm(x) > norm(y) ? -1 : 1
-  sign * √( (x[1]-y[1])^2 + (x[2]-y[2])^2 )
+
+function compute_displacement(
+    cps::AbstractVector{<:Point},
+    phi::AlgoimCallLevelSetFunction,
+    fun,
+    dt::Float64,
+    Ω::Triangulation)
+  # Note that cps must be (scalar) DoF-numbered, not lexicographic-numbered
+  searchmethod = KDTreeSearch()
+  cache1 = _point_to_cell_cache(searchmethod,Ω)
+  x_to_cell(x) = _point_to_cell!(cache1, x)
+  point_to_cell = lazy_map(x_to_cell, cps)
+  cell_to_points, _ = make_inverse_table(point_to_cell, num_cells(Ω))
+  cell_to_xs = lazy_map(Broadcasting(Reindex(cps)), cell_to_points)
+  cell_point_xs = CellPoint(cell_to_xs, Ω, PhysicalDomain())
+  fun_xs = evaluate(fun,cell_point_xs)
+  ∇φ_xs = evaluate(phi.∇φ,cell_point_xs)
+  cell_point_disp = lazy_map(Broadcasting(⋅),fun_xs,∇φ_xs)
+  cache_vals = array_cache(cell_point_disp)
+  cache_ctop = array_cache(cell_to_points)
+  disps = zeros(Float64,length(cps))
+  for cell in 1:length(cell_to_points)
+    pts = getindex!(cache_ctop,cell_to_points,cell)
+    vals = getindex!(cache_vals,cell_point_disp,cell)
+    for (i,pt) in enumerate(pts)
+      val = vals[i]
+      disps[pt] = dt * val
+    end
+  end
+  disps
+end
+
+signed_distance(φ::Function,x,y) = sign(φ(y))*norm(x-y)
+
+signed_distance(φ::T,x,y) where {T<:Number} = sign(φ)*norm(x-y)
+
+function _compute_signed_distance(
+    φ::AlgoimCallLevelSetFunction{<:Function,<:Function},
+    cps::Vector{<:Point},cos::Matrix{<:Point})
+  _dist(x,y) = signed_distance(φ.φ,x,y)
+  lazy_map(_dist,cps,cos)
+end
+
+function _compute_signed_distance(
+    φ::AlgoimCallLevelSetFunction{<:CellField,<:CellField},
+    cps::Vector{<:Point},cos::Matrix{<:Point})
+  φs = get_free_dof_values(φ.φ)
+  lazy_map(signed_distance,φs,cps,cos)
 end
 
 function compute_distance_fe_function(
@@ -411,7 +456,7 @@ function compute_distance_fe_function(
   rmodel = refine(bgmodel,order)
   cos = get_node_coordinates(rmodel)
   cos = node_to_dof_order(cos,fespace,rmodel,order)
-  dists = lazy_map(_dist,cps,cos)
+  dists = _compute_signed_distance(φ,cps,cos) 
   FEFunction(fespace,dists)
 end
 
diff --git a/src/Exports.jl b/src/Exports.jl
index 0b22391..37a6fcc 100644
--- a/src/Exports.jl
+++ b/src/Exports.jl
@@ -58,6 +58,7 @@ end
 @publish AlgoimUtils fill_cpp_data
 @publish AlgoimUtils fill_cpp_data_raw
 @publish AlgoimUtils compute_closest_point_projections
+@publish AlgoimUtils compute_displacement
 @publish AlgoimUtils compute_normal_displacement
 @publish AlgoimUtils compute_distance_fe_function
 @publish AlgoimUtils delaunaytrian
diff --git a/test/AlgoimUtilsTests/DualQuadraturesTests.jl b/test/AlgoimUtilsTests/DualQuadraturesTests.jl
index 6d5b61c..df0b269 100644
--- a/test/AlgoimUtilsTests/DualQuadraturesTests.jl
+++ b/test/AlgoimUtilsTests/DualQuadraturesTests.jl
@@ -1,4 +1,4 @@
-module AlgoimInterfaceTests
+module DualQuadraturesTests
 
 using Test
 using Gridap
diff --git a/test/AlgoimUtilsTests/TwoPhaseStokesAlgoimTests.jl b/test/AlgoimUtilsTests/TwoPhaseStokesAlgoimTests.jl
index b0f4da8..2fff55f 100644
--- a/test/AlgoimUtilsTests/TwoPhaseStokesAlgoimTests.jl
+++ b/test/AlgoimUtilsTests/TwoPhaseStokesAlgoimTests.jl
@@ -8,7 +8,7 @@ using GridapPardiso
 using SparseMatricesCSR
 using Test
 
-function main(n::Int,w::Float64,νˡ::Float64,νˢ::Float64,γ::Float64)
+function steady_state(n::Int,w::Float64,νˡ::Float64,νˢ::Float64,γ::Float64)
 
   pmin = Point(-1.5,-1.0)
   pmax = Point( 4.5, 1.0)
@@ -147,14 +147,219 @@ function main(n::Int,w::Float64,νˡ::Float64,νˢ::Float64,γ::Float64)
 
 end
 
-n = 80
-w = 1.0e-1
+function dynamic(n::Int,w::Float64,νˡ::Float64,νˢ::Float64,γ::Float64,Δt₀::Float64)
+
+  pmin = Point(-1.5,-2.5)
+  pmax = Point( 3.5, 2.5)
+  partition = (n,n)
+
+  model = CartesianDiscreteModel(pmin,pmax,partition)
+  Ω = Triangulation(model)
+  dp = pmax - pmin
+  h = dp[2]/n
+
+  order = 2
+  degree = order == 1 ? 3 : 2*order
+
+  R = 0.12
+  φ(x) = 1.0 - ( ( x[1]*x[1] + x[2]*x[2] ) / R^2 )
+  reffeᵠ = ReferenceFE(lagrangian,Float64,order+1)
+  Vbg = TestFESpace(Ω,reffeᵠ)
+
+  # Buffer of active model and integration objects
+  buffer = Ref{Any}((Ωˡ=nothing,Ωˢ=nothing,
+                     dΩˡ=nothing,dΩˢ=nothing,
+                     dΓ=nothing,n_Γ=nothing,
+                     aggsˡ=nothing,aggsˢ=nothing,
+                     φ₋=nothing,cp₋=nothing,t=nothing))
+
+  function update_buffer!(t,dt,v₋₂)
+
+    if buffer[].t == t
+      return true
+    else
+
+      if buffer[].φ₋ === nothing
+        __φ₋ = AlgoimCallLevelSetFunction(φ,∇(φ))
+        _φ₋ = compute_distance_fe_function(model,Vbg,__φ₋,order+1,cppdegree=3)
+      else
+        cp₋₂ = buffer[].cp₋
+        φ₋₂ = buffer[].φ₋
+        ϕ₋ = compute_displacement(cp₋₂,φ₋₂,v₋₂,dt,Ω)
+        ϕ₋ = get_free_dof_values(φ₋₂.φ) - ϕ₋
+        _φ₋ = FEFunction(Vbg,ϕ₋)
+      end
+
+      φ₋ = AlgoimCallLevelSetFunction(_φ₋,∇(_φ₋))
+      cp₋ = compute_closest_point_projections(Vbg,φ₋,order+1,cppdegree=3)
+
+      lquad = Quadrature(algoim,φ₋,degree,phase=IN)
+      Ωˡ,dΩˡ,cell_to_is_liquid = TriangulationAndMeasure(Ω,lquad)
+      
+      squad = Quadrature(algoim,φ₋,degree,phase=OUT)
+      Ωˢ,dΩˢ,cell_to_is_solid = TriangulationAndMeasure(Ω,squad)
+    
+      iquad = Quadrature(algoim,φ₋,degree)
+      _,dΓ,cell_to_is_cut = TriangulationAndMeasure(Ω,iquad)
+      n_Γ = normal(φ₋,Ω) # Exterior to liquid
+    
+      aggsˡ = aggregate(Ω,cell_to_is_liquid,cell_to_is_cut,IN)
+      aggsˢ = aggregate(Ω,cell_to_is_solid,cell_to_is_cut,OUT)
+
+      buffer[] = (Ωˡ=Ωˡ,Ωˢ=Ωˢ,dΩˡ=dΩˡ,dΩˢ=dΩˢ,
+                  dΓ=dΓ,n_Γ=n_Γ,aggsˡ=aggsˡ,aggsˢ=aggsˢ,
+                  φ₋=φ₋,cp₋=cp₋,t=t)
+      return true
+
+    end
+
+  end
+
+  reffeᵘ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
+  reffeˢ = ReferenceFE(lagrangian,VectorValue{2,Float64},order,space=:S)
+  reffeᵖ = ReferenceFE(lagrangian,Float64,order-1,space=:P)
+
+  labels = get_face_labeling(model)
+  add_tag_from_tags!(labels,"inlet",[7])
+
+  function update_all!(t::Real,dt::Real,disp)
+
+    update_buffer!(t,dt,disp)
+
+    Ωˡ  = buffer[].Ωˡ
+    Ωˢ  = buffer[].Ωˢ
+ 
+    dΩˡ = buffer[].dΩˡ
+    dΩˢ = buffer[].dΩˢ
+ 
+    dΓ  = buffer[].dΓ
+    n_Γ = buffer[].n_Γ
+
+    aggsˡ = buffer[].aggsˡ
+    aggsˢ = buffer[].aggsˢ
+
+    Vˡstd = TestFESpace(Ωˡ,reffeᵘ,dirichlet_tags=["inlet"])
+    Vˡser = TestFESpace(Ωˡ,reffeˢ,conformity=:L2)
+    Qˡstd = TestFESpace(Ωˡ,reffeᵖ)
+  
+    Vˢstd = TestFESpace(Ωˢ,reffeᵘ)
+    Vˢser = TestFESpace(Ωˢ,reffeˢ,conformity=:L2)
+    Qˢstd = TestFESpace(Ωˢ,reffeᵖ)
+  
+    Vˡ = AgFEMSpace(Vˡstd,aggsˡ,Vˡser)
+    Qˡ = AgFEMSpace(Qˡstd,aggsˡ)
+    Vˢ = AgFEMSpace(Vˢstd,aggsˢ,Vˢser)
+    Qˢ = AgFEMSpace(Qˢstd,aggsˢ)  
+    K = ConstantFESpace(model)
+  
+    uᵢ(x) = VectorValue(w*(1.0-x[2]*x[2]),0.0)
+    Uˡ = TrialFESpace(Vˡ,[uᵢ])
+    Pˡ = TrialFESpace(Qˡ)
+    Uˢ = TrialFESpace(Vˢ)
+    Pˢ = TrialFESpace(Qˢ)
+    L = TrialFESpace(K)
+  
+    Y = MultiFieldFESpace([Vˡ,Qˡ,K,Vˢ,Qˢ,K])
+    X = MultiFieldFESpace([Uˡ,Pˡ,L,Uˢ,Pˢ,L])
+
+    X,Y,dΩˡ,dΩˢ,Ωˡ,Ωˢ,dΓ,n_Γ
+
+  end
+
+  t₀ = 0.0
+  Δt = Δt₀
+
+  X,Y,dΩˡ,dΩˢ,Ωˡ,Ωˢ,dΓ,n_Γ = update_all!(t₀,Δt,VectorValue(0.0,0.0))
+
+  wˡ = CellField(νˢ/(νˡ+νˢ),Ω)
+  wˢ = CellField(νˡ/(νˡ+νˢ),Ω)
+  νᵞ = 2*νˡ*νˢ/(νˡ+νˢ)
+
+  jumpᵘ(uˡ,uˢ) = uˡ-uˢ
+  σ(ε,ν) = 2*ν*ε
+  τ(ε) = one(ε)
+  meanᵗ(uˡ,pˡ,νˡ,uˢ,pˢ,νˢ) = 
+    wˡ*(σ(ε(uˡ),νˡ))-pˡ*(τ∘(ε(uˡ))) + 
+    wˢ*(σ(ε(uˢ),νˢ))-pˢ*(τ∘(ε(uˢ)))
+
+  a((uˡ,pˡ,lˡ,uˢ,pˢ,lˢ),(vˡ,qˡ,kˡ,vˢ,qˢ,kˢ)) =
+    ∫( 2*νˡ*(ε(uˡ)⊙ε(vˡ)) - qˡ*(∇⋅uˡ) - (∇⋅vˡ)*pˡ )dΩˡ +
+    ∫( 2*νˢ*(ε(uˢ)⊙ε(vˢ)) - qˢ*(∇⋅uˢ) - (∇⋅vˢ)*pˢ )dΩˢ +
+    ∫( pˡ*kˡ )dΩˡ + ∫( qˡ*lˡ )dΩˡ +
+    ∫( pˢ*kˢ )dΩˢ + ∫( qˢ*lˢ )dΩˢ +
+    ∫( (γ*νᵞ/h)*(jumpᵘ(uˡ,uˢ)⋅jumpᵘ(vˡ,vˢ)) - 
+       (jumpᵘ(vˡ,vˢ)⋅(n_Γ⋅meanᵗ(uˡ,pˡ,νˡ,uˢ,pˢ,νˢ))) - 
+       ((n_Γ⋅meanᵗ(vˡ,qˡ,νˡ,vˢ,qˢ,νˢ))⋅jumpᵘ(uˡ,uˢ)) )dΓ
+
+  l((vˡ,qˡ,kˡ,vˢ,qˢ,kˢ)) = 0.0
+
+  assem = SparseMatrixAssembler(SymSparseMatrixCSR{1,Float64,Int},Vector{Float64},X,Y)
+  op = AffineFEOperator(a,l,X,Y,assem)
+
+  A = get_matrix(op)
+  b = get_vector(op)
+  x = similar(b)
+  msglvl = 0
+  # Customizing solver for real symmetric indefinite matrices
+  # See https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface/pardiso-iparm-parameter.html
+  iparm = new_iparm()
+  iparm[0+1] = 1 # Use default values (0 = enabled)
+  iparm[1+1] = 3 # Fill-in reducing ordering for the input matrix
+  iparm[9+1] = 8 # Pivoting perturbation
+  iparm[10+1] = 1 # Scaling vectors
+  iparm[12+1] = 1 # Improved accuracy using (non-) symmetric weighted matching
+  # Note that in the analysis phase (phase=11) you must provide the numerical values 
+  # of the matrix A in array a in case of scaling and symmetric weighted matching.
+  iparm[17+1] = -1 # Report the number of non-zero elements in the factors. 
+  iparm[18+1] = -1 # Report number of floating point operations (in 10^6 floating point operations) that are necessary to factor the matrix A.
+  iparm[20+1] = 1 # Pivoting for symmetric indefinite matrices. 
+  ps = PardisoSolver(GridapPardiso.MTYPE_REAL_SYMMETRIC_INDEFINITE, iparm, msglvl)
+  ss = symbolic_setup(ps, A)
+  ns = numerical_setup(ss, A)
+  solve!(x, ns, b)
+  xh = FEFunction(X,x)
+  uhl, phl, _, uhs, phs, _ = xh
+
+  _A = get_matrix(op)
+  _b = get_vector(op)
+  _x = get_free_dof_values(xh)
+  _r = _A*_x - _b
+  nr = norm(_r)
+  nb = norm(_b)
+  nx = norm(_x)
+  # @show nr, nr/nb, nr/nx
+  tol_warn  = 1.0e-10
+  if nr > tol_warn && nr/nb > tol_warn && nr/nx > tol_warn
+    @warn "Solver not accurate"
+  end
+
+  # colors = color_aggregates(aggsˡ,model)
+  # writevtk(Ω,"res_bg_l",celldata=["aggregate"=>aggsˡ,"color"=>colors])
+  # writevtk(Ω,"res_bg_s",celldata=["aggregate"=>aggsˢ,"color"=>colors])
+  writevtk(Ωˡ,"res_l",cellfields=["uhl"=>uhl,"phl"=>phl])
+  writevtk(dΩˢ,"res_s",cellfields=["uhs"=>uhs,"phs"=>phs])
+  writevtk(dΓ,"res_gam",cellfields=["uhl"=>uhl,"uhs"=>uhs],qhulltype=convexhull)
+  # nh = interpolate_everywhere(n_Γ,Vˡstd)  
+  # σn = νˡ*(ε(uhl)⋅nh) # - phl*nh
+  # writevtk(dΓ,"res_gam",cellfields=["uhl"=>uhl,"phl"=>phl,"sn"=>σn],qhulltype=convexhull)
+
+  X,Y,dΩˡ,dΩˢ,Ωˡ,Ωˢ,dΓ,n_Γ = update_all!(t₀+Δt,Δt,uhs)
+  
+  writevtk(Ωˡ,"res_l_2")
+  writevtk(dΩˢ,"res_s_2")
+  writevtk(dΓ,"res_gam_2",qhulltype=convexhull)
+
+end
+
+n  = 80
+w  = 1.0e-1
 νˡ = 1.0e-1
 νˢ = 1.0e+2
-γ = 1.0e-1 # Scale with νˢ
+γ  = 1.0e-1 # Scale with νˢ
+Δt = 1.0
 
 @info "Values: n = $n, w = $w, νˡ = $νˡ, νˢ = $νˢ, γ = $γ"
 
-main(n,w,νˡ,νˢ,γ)
+dynamic(n,w,νˡ,νˢ,γ,Δt)
 
 end # module
\ No newline at end of file