diff --git a/scripts/Embedded/Examples/FCM_2d_thermal_with_NonDesignDom.jl b/scripts/Embedded/Examples/FCM_2d_thermal_with_NonDesignDom.jl
new file mode 100644
index 0000000..d4c7162
--- /dev/null
+++ b/scripts/Embedded/Examples/FCM_2d_thermal_with_NonDesignDom.jl
@@ -0,0 +1,114 @@
+using Gridap,GridapTopOpt, GridapSolvers
+using Gridap.Adaptivity, Gridap.Geometry
+using GridapEmbedded, GridapEmbedded.LevelSetCutters
+
+using GridapTopOpt: StateParamIntegrandWithMeasure
+
+path="./results/FCM_thermal_compliance_ALM/"
+rm(path,force=true,recursive=true)
+mkpath(path)
+n = 50
+order = 1
+γ = 0.1
+max_steps = floor(Int,order*n/5)
+vf = 0.4
+α_coeff = 4max_steps*γ
+iter_mod = 1
+
+_model = CartesianDiscreteModel((0,1,0,1),(n,n))
+base_model = UnstructuredDiscreteModel(_model)
+ref_model = refine(base_model, refinement_method = "barycentric")
+model = ref_model.model
+el_Δ = get_el_Δ(_model)
+h = maximum(el_Δ)
+h_refine = maximum(el_Δ)/2
+f_Γ_D(x) = (x[1] ≈ 0.0 && (x[2] <= 0.2 + eps() || x[2] >= 0.8 - eps()))
+f_Γ_N(x) = (x[1] ≈ 1 && 0.4 - eps() <= x[2] <= 0.6 + eps())
+f_NonDesign(x) = (0.4 <= x[1] <= 0.6 && 0.0 <= x[2] <= 0.1)
+update_labels!(1,model,f_Γ_D,"Gamma_D")
+update_labels!(2,model,f_Γ_N,"Gamma_N")
+update_labels!(3,model,f_NonDesign,"NonDesignDom")
+
+## Triangulations and measures
+Ω = Triangulation(model)
+Γ_N = BoundaryTriangulation(model,tags="Gamma_N")
+dΩ = Measure(Ω,2*order)
+dΓ_N = Measure(Γ_N,2*order)
+vol_D = sum(∫(1)dΩ)
+
+## Levet-set function space and derivative regularisation space
+reffe_scalar = ReferenceFE(lagrangian,Float64,order)
+V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_N","NonDesignDom"])
+U_reg = TrialFESpace(V_reg,0)
+V_φ = TestFESpace(model,reffe_scalar)
+
+## Levet-set function
+φh = interpolate(x->-cos(4π*x[1])*cos(4π*x[2])-0.4,V_φ)
+Ωs = EmbeddedCollection(model,φh) do cutgeo,_
+  Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_IN),V_φ)
+  Ωout = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_OUT),V_φ)
+  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
+  Γg = GhostSkeleton(cutgeo)
+  Ωact = Triangulation(cutgeo,ACTIVE)
+  (;
+    :Ωin  => Ωin,
+    :dΩin => Measure(Ωin,2*order),
+    :Ωout  => Ωout,
+    :dΩout => Measure(Ωout,2*order),
+    :Γg   => Γg,
+    :dΓg  => Measure(Γg,2*order),
+    :n_Γg => get_normal_vector(Γg),
+    :Γ    => Γ,
+    :dΓ   => Measure(Γ,2*order),
+    :Ωact => Ωact
+  )
+end
+
+## Weak form
+const ϵ = 1e-3
+a(u,v,φ) = ∫(∇(v)⋅∇(u))Ωs.dΩin + ∫(ϵ*∇(v)⋅∇(u))Ωs.dΩout
+l(v,φ) = ∫(v)dΓ_N
+
+## Optimisation functionals
+J(u,φ) = a(u,u,φ)
+Vol(u,φ) = ∫(1/vol_D)Ωs.dΩin - ∫(vf/vol_D)dΩ
+dVol(q,u,φ) = ∫(-1/vol_D*q/(norm ∘ (∇(φ))))Ωs.dΓ
+
+## Setup solver and FE operators
+V = TestFESpace(Ω,reffe_scalar;dirichlet_tags=["Gamma_D"])
+U = TrialFESpace(V,0.0)
+state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
+pcfs = PDEConstrainedFunctionals(J,[Vol],state_map;analytic_dC=(dVol,))
+
+## Evolution Method
+evo = CutFEMEvolve(V_φ,Ωs,dΩ,h;max_steps,γg=0.5)
+reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=ArtificialViscosity(3.0))
+ls_evo = UnfittedFEEvolution(evo,reinit)
+reinit!(ls_evo,φh)
+
+## Hilbertian extension-regularisation problems
+α = α_coeff*(h_refine/order)^2
+a_hilb(p,q) =∫(α*∇(p)⋅∇(q) + p*q)dΩ;
+vel_ext = VelocityExtension(a_hilb,U_reg,V_reg)
+
+## Optimiser
+converged(m) = GridapTopOpt.default_al_converged(
+  m;
+  L_tol = 0.01*h_refine,
+  C_tol = 0.01
+)
+optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;debug=true,
+  γ,verbose=true,constraint_names=[:Vol],converged)
+for (it,uh,φh,state) in optimiser
+  x_φ = get_free_dof_values(φh)
+  idx = findall(isapprox(0.0;atol=10^-10),x_φ)
+  !isempty(idx) && @warn "Boundary intersects nodes!"
+  if iszero(it % iter_mod)
+    writevtk(Ω,path*"Omega$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh,"velh"=>FEFunction(V_φ,state.vel)])
+    writevtk(Ωs.Ωin,path*"Omega_in$it",cellfields=["uh"=>uh])
+  end
+  write_history(path*"/history.txt",optimiser.history)
+end
+it = get_history(optimiser).niter; uh = get_state(pcfs)
+writevtk(Ω,path*"Omega$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+writevtk(Ωs.Ωin,path*"Omega_in$it",cellfields=["uh"=>uh])
\ No newline at end of file
diff --git a/scripts/Embedded/Examples/fsi/TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl b/scripts/Embedded/Examples/fsi/TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl
new file mode 100644
index 0000000..645c2be
--- /dev/null
+++ b/scripts/Embedded/Examples/fsi/TO-6-Brinkmann_stokes_P2-P1_Ersatz_elast_fsi.jl
@@ -0,0 +1,141 @@
+using Gridap, Gridap.Geometry, Gridap.Adaptivity
+using GridapEmbedded, GridapEmbedded.LevelSetCutters
+using GridapTopOpt
+
+path = "./results/fsi testing/"
+mkpath(path)
+
+# Cut the background model
+L = 1.0#2.0
+H = 1.0#0.5
+x0 = 0.5
+l = 0.4
+w = 0.1
+a = 0.5#0.3
+b = 0.01
+
+n = 100
+partition = (n,n)
+D = length(partition)
+_model = CartesianDiscreteModel((0,L,0,H),partition)
+base_model = UnstructuredDiscreteModel(_model)
+ref_model = refine(base_model, refinement_method = "barycentric")
+model = ref_model.model
+
+el_Δ = get_el_Δ(_model)
+h = maximum(el_Δ)
+f_Γ_D(x) = x[1] ≈ 0
+f_Γ_NoSlipTop(x) = x[2] ≈ H
+f_Γ_NoSlipBottom(x) = x[2] ≈ 0
+f_NonDesign(x) = ((x0 - w/2 <= x[1] <= x0 + w/2 && 0.0 <= x[2] <= a) &&
+  x0 - w/2 <= x[1] <= x0 + w/2 && 0.0 <= x[2] <= b)
+update_labels!(1,model,f_Γ_D,"Gamma_D")
+update_labels!(2,model,f_Γ_NoSlipTop,"Gamma_NoSlipTop")
+update_labels!(3,model,f_Γ_NoSlipBottom,"Gamma_NoSlipBottom")
+update_labels!(4,model,f_NonDesign,"NonDesign")
+
+# Cut the background model
+reffe_scalar = ReferenceFE(lagrangian,Float64,1)
+V_φ = TestFESpace(model,reffe_scalar)
+φh = interpolate(x->-max(20*abs(x[1]-0.5),3*abs(x[2]-0.2))+1,V_φ)
+geo = DiscreteGeometry(φh,model)
+cutgeo = cut(model,geo)
+cutgeo_facets = cut_facets(model,geo)
+
+# Generate the "active" model
+Ω_act = Triangulation(model)
+
+# Setup integration meshes
+Ω = Triangulation(cutgeo,PHYSICAL)
+Ωout = Triangulation(cutgeo,PHYSICAL_OUT)
+Γ = EmbeddedBoundary(cutgeo)
+Γg = GhostSkeleton(cutgeo)
+Γi = SkeletonTriangulation(cutgeo_facets)
+
+# Setup normal vectors
+n_Γ = get_normal_vector(Γ)
+n_Γg = get_normal_vector(Γg)
+n_Γi = get_normal_vector(Γi)
+
+# Setup Lebesgue measures
+order = 2
+degree = 2*order
+dΩ = Measure(Ω,degree)
+dΩout = Measure(Ωout,degree)
+dΓ = Measure(Γ,degree)
+dΓg = Measure(Γg,degree)
+dΓi = Measure(Γi,degree)
+
+# Setup FESpace
+
+uin(x) = VectorValue(x[2]*(1-x[2]),0.0)
+
+reffe_u = ReferenceFE(lagrangian,VectorValue{D,Float64},order,space=:P)
+reffe_p = ReferenceFE(lagrangian,Float64,order-1,space=:P)
+reffe_d = ReferenceFE(lagrangian,VectorValue{D,Float64},order)
+
+V = TestFESpace(Ω_act,reffe_u,conformity=:H1,dirichlet_tags=["Gamma_D","Gamma_NoSlipTop","Gamma_NoSlipBottom"])
+Q = TestFESpace(Ω_act,reffe_p,conformity=:C0)
+T = TestFESpace(Ω_act ,reffe_d,conformity=:H1,dirichlet_tags=["Gamma_NoSlipBottom"])
+
+U = TrialFESpace(V,[uin,VectorValue(0.0,0.0),VectorValue(0.0,0.0)])
+P = TrialFESpace(Q)
+R = TrialFESpace(T)
+
+X = MultiFieldFESpace([U,P,R])
+Y = MultiFieldFESpace([V,Q,T])
+
+# Weak form
+## Fluid
+# Properties
+Re = 60 # Reynolds number
+ρ = 1.0 # Density
+L = 1.0 # Characteristic length
+u0_max = maximum(abs,get_dirichlet_dof_values(U))
+μ = ρ*L*u0_max/Re # Viscosity
+# Stabilization parameters
+γ = 1000.0
+
+# Terms
+σf_n(u,p) = μ*∇(u)⋅n_Γ - p*n_Γ
+a_Ω(u,v) = μ*(∇(u) ⊙ ∇(v))
+b_Ω(v,p) = - (∇⋅v)*p
+
+a_fluid((u,p),(v,q)) =
+  ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p)) * dΩ +
+  ∫( a_Ω(u,v)+b_Ω(u,q)+b_Ω(v,p) + (γ/h)*u⋅v ) * dΩout
+
+## Structure
+# Stabilization and material parameters
+function lame_parameters(E,ν)
+  λ = (E*ν)/((1+ν)*(1-2*ν))
+  μ = E/(2*(1+ν))
+  (λ, μ)
+end
+λs, μs = lame_parameters(1.0,0.3)
+ϵ = (λs + 2μs)*1e-3
+# Terms
+σ(ε) = λs*tr(ε)*one(ε) + 2*μs*ε
+a_solid(d,s) = ∫(ε(s) ⊙ (σ ∘ ε(d)))dΩout +
+  ∫(ϵ*(ε(s) ⊙ (σ ∘ ε(d))))dΩ # Ersatz
+
+## Full problem
+a((u,p,d),(v,q,s)) = a_fluid((u,p),(v,q)) + a_solid(d,s) +
+  ∫(σf_n(u,p) ⋅ s)dΓ # plus sign because of the normal direction
+l((v,q,s)) = 0.0
+
+op = AffineFEOperator(a,l,X,Y)
+
+uh, ph, dh = solve(op)
+
+# Mass flow rate through surface (this should be close to zero)
+@show m = sum(∫(ρ*uh⋅n_Γ)dΓ)
+
+writevtk(Ω_act,path*"fsi-stokes-brinkmann-P2P1_elast-ersatz_full",
+  cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh])
+writevtk(Ω,path*"fsi-stokes-brinkmann-P2P1_elast-ersatz_fluid",
+  cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh])
+writevtk(Ωout,path*"fsi-stokes-brinkmann-P2P1_elast-ersatz_solid",
+  cellfields=["uh"=>uh,"ph"=>ph,"dh"=>dh])
+
+writevtk(Γ,path*"fsi-stokes-brinkmann-P2P1_elast-ersatz_interface",cellfields=["σ⋅n"=>(σ ∘ ε(dh))⋅n_Γ,"σf_n"=>σf_n(uh,ph)])
\ No newline at end of file
diff --git a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve_new.jl b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve_new.jl
new file mode 100644
index 0000000..d1cb215
--- /dev/null
+++ b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve_new.jl
@@ -0,0 +1,174 @@
+"""
+    mutable struct CutFEMEvolve{V,M} <: Evolver
+
+CutFEM method for level-set evolution based on method developed by
+Burman et al. (2017). DOI: `10.1016/j.cma.2017.09.005`.
+
+# Parameters
+- `ode_solver::ODESolver`: ODE solver
+- `Ωs::B`: `EmbeddedCollection` holding updatable triangulation and measures from GridapEmbedded
+- `dΩ_bg::C`: Measure for integration
+- `space::B`: Level-set FE space
+- `assembler::Assembler`: FE assembler
+- `params::D`: Tuple of stabilisation parameter `γg`, mesh size `h`, and
+  max steps `max_steps`, and background mesh skeleton parameters
+- `cache`: Cache for evolver, initially `nothing`.
+
+# Note
+- The stepsize `dt = 0.1` in `RungeKutta` is a place-holder and is updated using
+  the `γ` passed to `solve!`.
+"""
+mutable struct CutFEMEvolve{A,B} <: Evolver
+  ode_solver::ODESolver
+  Ωs::EmbeddedCollection
+  space::A
+  assembler::Assembler
+  params::B
+  cache
+  function CutFEMEvolve(V_φ::A,Ωs::EmbeddedCollection,h::Real;
+      max_steps=10,
+      γg = 0.1,
+      ode_ls = LUSolver(),
+      ode_nl = NLSolver(ode_ls, show_trace=false, method=:newton, iterations=10),
+      ode_solver = MutableRungeKutta(ode_nl, ode_ls, 0.1, :DIRK_CrankNicolson_2_2),
+      assembler=SparseMatrixAssembler(V_φ,V_φ),
+      Γg_quad_degree = 2get_order(V_φ),
+      Ωφ_quad_degree = 2get_order(V_φ)) where A
+
+    Γg = SkeletonTriangulation(get_triangulation(V_φ))
+    dΓg = Measure(Γg,Γg_quad_degree)
+    n_Γg = get_normal_vector(Γg)
+    Ωφ = get_triangulation(V_φ)
+    dΩφ = Measure(Ωφ,Ωφ_quad_degree)
+    params = (;γg,h,max_steps,dΓg,n_Γg,dΩφ)
+    new{A,typeof(params)}(ode_solver,Ωs,V_φ,assembler,params,nothing)
+  end
+end
+
+get_ode_solver(s::CutFEMEvolve) = s.ode_solver
+get_assembler(s::CutFEMEvolve) = s.assembler
+get_space(s::CutFEMEvolve) = s.space
+get_measure(s::CutFEMEvolve) = get_params(s).dΩφ
+get_params(s::CutFEMEvolve) = s.params
+get_cache(s::CutFEMEvolve) = s.cache
+
+function get_transient_operator(φh::CellField,velh::CellField,s::CutFEMEvolve)
+  V_φ, assembler, params = s.space, s.assembler, s.params
+  γg, h, dΓg, n_Γg, dΩ_bg = params.γg, params.h, params.dΓg, params.n_Γg, params.dΩφ
+  ϵ = 1e-20
+
+  β = velh*∇(φh)/(ϵ + norm ∘ ∇(φh))
+  stiffness(t,u,v) = ∫((β ⋅ ∇(u)) * v)dΩ_bg + ∫(γg*h^2*jump(∇(u) ⋅ n_Γg)*jump(∇(v) ⋅ n_Γg))dΓg
+  mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ_bg
+  forcing(t,v) = ∫(0v)dΩ_bg + ∫(0*jump(∇(v) ⋅ n_Γg))dΓg
+  # Second term is added to address the following issue:
+  #  - ODEs is allocating separately the residual and jacobian
+  #  - This is fine in serial, but in parallel there are some instances where the the following happens:
+  #     - The residual is touched by less ghost entries than the columns of the matrix
+  #     - If we assemble both jac and res together, we communicate the extra ghost ids to
+  #       the residual, so everything is consistent.
+  #     - However, if we assemble the residual and jacobian separately,
+  #       the residual is not aware of the extra ghost ids
+  # This happens when there are touched ghost entries that do not belong to the local domain.
+  # In particular, this happens when we have jumps, where some contributions come from two
+  # cells away. Boundary cells then get contributions from cells which are not in the local domain.
+  Ut_φ = TransientTrialFESpace(V_φ)
+  ode_op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ;
+    constant_forms=(false,true),assembler)
+  return ode_op
+end
+
+function update_reuse!(state,reuse_new;zero_tF=false)
+  U, (tF, stateF, state0, uF, odecache) = state
+  odeslvrcache, odeopcache = odecache
+  _, ui_pre, slopes, J, r, sysslvrcaches = odeslvrcache
+
+  odeslvrcache_new = (reuse_new, ui_pre, slopes, J, r, sysslvrcaches)
+  odecache_new = odeslvrcache_new, odeopcache
+  _tF = zero_tF ? 0.0 : tF
+  return U, (_tF, stateF, state0, uF, odecache_new)
+end
+
+function solve!(s::CutFEMEvolve,φh::CellField,velh::CellField,γ,cache::Nothing)
+  ode_solver = get_ode_solver(s)
+  params = get_params(s)
+  h, max_steps = params.h, params.max_steps
+
+  # Setup FE operator and solver
+  ode_op = get_transient_operator(φh,velh,s)
+  dt = _compute_Δt(h,γ,get_free_dof_values(velh))
+  ode_solver.dt = dt
+  ode_sol = solve(ode_solver,ode_op,0.0,dt*max_steps,interpolate(φh,s.space))
+
+  # March
+  march = Base.iterate(ode_sol)
+  data, state = march
+  state_new = update_reuse!(state,true)
+  march_new = data, state_new
+  while march_new !== nothing
+    data, state_new = march_new
+    march_new = Base.iterate(ode_sol,state_new)
+  end
+
+  # Update φh and cache
+  _, φhF = data
+  if get_triangulation(φh) === get_triangulation(φhF)
+    copy!(get_free_dof_values(φh),get_free_dof_values(φhF))
+  else
+    interpolate!(φhF,get_free_dof_values(φh),get_fe_space(φh))
+  end
+  writevtk(get_triangulation(φhF),"./results/FCM_thermal_compliance_ALM/test",cellfields=["φ"=>φhF])
+  writevtk(get_triangulation(φh),"./results/FCM_thermal_compliance_ALM/test2",cellfields=["φ"=>φh])
+  s.cache = state_new
+  update_collection!(s.Ωs,φh) # TODO: remove?
+  return φh
+end
+
+function solve!(s::CutFEMEvolve,φh::CellField,velh::CellField,γ,cache)
+  ode_solver = get_ode_solver(s)
+  params = get_params(s)
+  h, max_steps = params.h, params.max_steps
+
+  ## Update state
+  # `get_transient_operator` re-creates the entire TransientLinearFEOperator wrapper.
+  #   We do this so that the first iterate of ODESolution always recomputes the
+  #   stiffness matrix and associated the Jacboian, numerical setups, etc via
+  #   `constant_forms = (false,true)`.
+  ode_op = get_transient_operator(φh,velh,s)
+  # Between the first iterate and subsequent iterates we use the function
+  #   `update_reuse!` to update the iterator state so that we re-use
+  #   the stiffness matrix, etc. The Optional argument `zero_tF` indicates
+  #   whether we are solving a new ODE with the same functional form but
+  #   updated coefficients in the weak form. If so, we want to re-use the cache.
+  state_inter = update_reuse!(cache,false;zero_tF=true)
+  dt = _compute_Δt(h,γ,get_free_dof_values(velh))
+  ode_solver.dt = dt
+
+  ## March
+  ode_sol = solve(ode_solver,ode_op,0.0,dt*max_steps,interpolate(φh,s.space))
+  march = Base.iterate(ode_sol,state_inter) # First step includes stiffness matrix update
+  data, state = march
+  state_updated = update_reuse!(state,true) # Fix the stiffness matrix for remaining march
+  march_updated = data, state_updated
+  while march_updated !== nothing
+    data, state_updated = march_updated
+    march_updated = Base.iterate(ode_sol,state_updated)
+  end
+
+  ## Update φh and cache
+  _, φhF = data
+  if get_triangulation(φh) === get_triangulation(φhF)
+    copy!(get_free_dof_values(φh),get_free_dof_values(φhF))
+  else
+    interpolate!(φhF,get_free_dof_values(φh),get_fe_space(φh))
+  end
+  s.cache = state_updated
+  update_collection!(s.Ωs,φh) # TODO: remove?
+  return φh
+end
+
+function _compute_Δt(h,γ,vel)
+  T = eltype(γ)
+  v_norm = maximum(abs,vel)
+  return γ * h / (eps(T)^2 + v_norm)
+end
\ No newline at end of file