diff --git a/scripts/Embedded/Bugs/FIXED_change_domain_bug.jl b/scripts/Embedded/Bugs/FIXED_change_domain_bug.jl
index 4f1be679..d770e5f7 100644
--- a/scripts/Embedded/Bugs/FIXED_change_domain_bug.jl
+++ b/scripts/Embedded/Bugs/FIXED_change_domain_bug.jl
@@ -25,7 +25,7 @@ dΓ = Measure(Γ,2order)
 
 
 _Γ = EmbeddedBoundary(cutgeo)
-Γ = DifferentiableTriangulation(_Γ)
+Γ = DifferentiableTriangulation(_Γ,V_φ)
 dΓ = Measure(Γ,2order)
 ∫(vh)dΓ
 ∫(2vh)dΓ
diff --git a/scripts/Embedded/Bugs/basis_isbitstype_bug.jl b/scripts/Embedded/Bugs/basis_isbitstype_bug.jl
index c19d000e..f298ca6f 100644
--- a/scripts/Embedded/Bugs/basis_isbitstype_bug.jl
+++ b/scripts/Embedded/Bugs/basis_isbitstype_bug.jl
@@ -41,7 +41,7 @@ cutgeo = cut(model,geo)
 # A.1) Volume integral
 
 Ω = Triangulation(cutgeo,PHYSICAL_IN)
-Ω_AD = DifferentiableTriangulation(Ω)
+Ω_AD = DifferentiableTriangulation(Ω,V_φ)
 dΩ = Measure(Ω_AD,2*order)
 
 Γ = EmbeddedBoundary(cutgeo)
diff --git a/scripts/Embedded/Bugs/interface_change_domain_bug.jl b/scripts/Embedded/Bugs/interface_change_domain_bug.jl
new file mode 100644
index 00000000..548bb560
--- /dev/null
+++ b/scripts/Embedded/Bugs/interface_change_domain_bug.jl
@@ -0,0 +1,45 @@
+using Test
+
+using GridapTopOpt
+using Gridap
+
+using GridapDistributed, PartitionedArrays
+
+using GridapEmbedded
+using GridapEmbedded.LevelSetCutters
+using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData, Gridap.Adaptivity
+
+using GridapTopOpt: get_subfacet_normal_vector, get_ghost_normal_vector
+using GridapTopOpt: get_conormal_vector, get_tangent_vector
+
+using GridapDistributed: DistributedTriangulation, DistributedDomainContribution
+
+order = 1
+n = 16
+
+parts = (2,2)
+ranks = DebugArray(LinearIndices((prod(parts),)))
+
+# _model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),(n,n))
+_model = CartesianDiscreteModel((0,1,0,1),(n,n))
+
+base_model = UnstructuredDiscreteModel(_model)
+ref_model = refine(base_model, refinement_method = "barycentric")
+model = Gridap.Adaptivity.get_model(ref_model)
+
+reffe = ReferenceFE(lagrangian,Float64,order)
+V_φ = TestFESpace(model,reffe)
+
+φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.5223,V_φ) # Circle
+fh = interpolate(x->cos(x[1]*x[2]),V_φ)
+
+geo = DiscreteGeometry(φh,model)
+cutgeo = cut(model,geo)
+
+Γ = EmbeddedBoundary(cutgeo)
+Γ_AD = DifferentiableTriangulation(Γ,V_φ)
+dΓ_AD = Measure(Γ_AD,2*order)
+
+J_int2(φ) = ∫(g(fh))dΓ_AD
+dJ_int_AD2 = gradient(J_int2,φh)
+dJ_int_AD_vec2 = assemble_vector(dJ_int_AD2,V_φ)
\ No newline at end of file
diff --git a/scripts/Embedded/Examples/2d_compliance_L.jl b/scripts/Embedded/Examples/2d_compliance_L.jl
index 920873a3..cc809631 100644
--- a/scripts/Embedded/Examples/2d_compliance_L.jl
+++ b/scripts/Embedded/Examples/2d_compliance_L.jl
@@ -56,11 +56,11 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-cos(8π*x[1])*cos(8π*x[2])-0.2,V_φ)
 Ωs = EmbeddedCollection(model,φh) do cutgeo,_
   Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL))
-  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
+  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
   n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Γg   => Γg,
@@ -70,7 +70,7 @@ V_φ = TestFESpace(model,reffe_scalar)
     :dΓ   => Measure(Γ,2*order),
     :Ωact => Ωact
   )
-end  
+end
 
 ## Weak form
 function lame_parameters(E,ν)
@@ -87,7 +87,7 @@ E = 1.0
 γg = 0.1
 
 g = VectorValue(0,-1)
-a(u,v,φ) = ∫(ε(v) ⊙ (σ ∘ ε(u)))Ωs.dΩin + 
+a(u,v,φ) = ∫(ε(v) ⊙ (σ ∘ ε(u)))Ωs.dΩin +
   # ∫((γg*h)*jump(nΛ_D⋅∇(v)) ⋅ jump(nΛ_D⋅∇(u)))dΛ_D + # <- this currently breaks due to `_compute_cell_ids` function
   ∫((γg*h^3)*jump(Ωs.n_Γg⋅∇(v)) ⋅ jump(Ωs.n_Γg⋅∇(u)))Ωs.dΓg
 l(v,φ) = ∫(v⋅g)dΓ_N
@@ -106,12 +106,12 @@ state_collection = EmbeddedCollection(model,φh) do _,_
   V = TestFESpace(Ωs.Ωact,reffe;dirichlet_tags=["Gamma_D"])
   U = TrialFESpace(V,VectorValue(0.0,0.0))
   state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh;ls,adjoint_ls=ls)
-  (; 
+  (;
     :state_map => state_map,
     :J => StateParamIntegrandWithMeasure(J,state_map),
     :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,])
   )
-end  
+end
 pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,))
 
 ## Evolution Method
diff --git a/scripts/Embedded/Examples/2d_thermal.jl b/scripts/Embedded/Examples/2d_thermal.jl
index 39622d1f..6ddee45f 100644
--- a/scripts/Embedded/Examples/2d_thermal.jl
+++ b/scripts/Embedded/Examples/2d_thermal.jl
@@ -44,11 +44,10 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-cos(8π*x[1])*cos(8π*x[2])-0.2,V_φ)
 Ωs = EmbeddedCollection(model,φh) do cutgeo,_
   Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL))
-  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
+  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Γg   => Γg,
@@ -58,7 +57,7 @@ V_φ = TestFESpace(model,reffe_scalar)
     :dΓ   => Measure(Γ,2*order),
     :Ωact => Ωact
   )
-end  
+end
 
 ## Weak form
 const γg = 0.1
@@ -78,12 +77,12 @@ state_collection = EmbeddedCollection(model,φh) do _,_
   V = TestFESpace(Ωs.Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"])
   U = TrialFESpace(V,0.0)
   state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh;ls,adjoint_ls=ls)
-  (; 
+  (;
     :state_map => state_map,
     :J => StateParamIntegrandWithMeasure(J,state_map),
     :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,])
   )
-end  
+end
 pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,))
 
 ## Evolution Method
diff --git a/scripts/Embedded/Examples/2d_thermal_AgFEM.jl b/scripts/Embedded/Examples/2d_thermal_AgFEM.jl
index 40a0a534..b3346c8a 100644
--- a/scripts/Embedded/Examples/2d_thermal_AgFEM.jl
+++ b/scripts/Embedded/Examples/2d_thermal_AgFEM.jl
@@ -44,11 +44,10 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-cos(8π*x[1])*cos(8π*x[2])-0.2,V_φ)
 Ωs = EmbeddedCollection(model,φh) do cutgeo,_
   Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL))
-  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
+  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Γg   => Γg,
@@ -58,7 +57,7 @@ V_φ = TestFESpace(model,reffe_scalar)
     :dΓ   => Measure(Γ,2*order),
     :Ωact => Ωact
   )
-end  
+end
 
 ## Weak form
 const γg = 0.1
@@ -81,12 +80,12 @@ state_collection = EmbeddedCollection(model,φh) do cutgeo,_
   V = AgFEMSpace(Vstd,aggregates)
   U = TrialFESpace(V,0.0)
   state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh;ls,adjoint_ls=ls)
-  (; 
+  (;
     :state_map => state_map,
     :J => StateParamIntegrandWithMeasure(J,state_map),
     :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,])
   )
-end  
+end
 pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,))
 
 ## Evolution Method
diff --git a/scripts/Embedded/Examples/2d_thermal_L.jl b/scripts/Embedded/Examples/2d_thermal_L.jl
index ba41ee31..87b388d7 100644
--- a/scripts/Embedded/Examples/2d_thermal_L.jl
+++ b/scripts/Embedded/Examples/2d_thermal_L.jl
@@ -52,11 +52,10 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-cos(8π*x[1])*cos(8π*x[2])-0.2,V_φ)
 Ωs = EmbeddedCollection(model,φh) do cutgeo,_
   Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL))
-  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
+  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Γg   => Γg,
@@ -66,7 +65,7 @@ V_φ = TestFESpace(model,reffe_scalar)
     :dΓ   => Measure(Γ,2*order),
     :Ωact => Ωact
   )
-end  
+end
 
 ## Weak form
 const γg = 0.1
@@ -86,12 +85,12 @@ state_collection = EmbeddedCollection(model,φh) do _,_
   V = TestFESpace(Ωs.Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"])
   U = TrialFESpace(V,0.0)
   state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh;ls,adjoint_ls=ls)
-  (; 
+  (;
     :state_map => state_map,
     :J => StateParamIntegrandWithMeasure(J,state_map),
     :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,])
   )
-end  
+end
 pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,))
 
 ## Evolution Method
diff --git a/scripts/Embedded/Examples/2d_thermal_gmsh_example.jl b/scripts/Embedded/Examples/2d_thermal_gmsh_example.jl
index f50d0034..ad9aabee 100644
--- a/scripts/Embedded/Examples/2d_thermal_gmsh_example.jl
+++ b/scripts/Embedded/Examples/2d_thermal_gmsh_example.jl
@@ -52,11 +52,10 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-cos(8π*x[1])*cos(8π*x[2])-0.3,V_φ)
 Ωs = EmbeddedCollection(model,φh) do cutgeo,_
   Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL))
-  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
+  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Γg   => Γg,
@@ -66,7 +65,7 @@ V_φ = TestFESpace(model,reffe_scalar)
     :dΓ   => Measure(Γ,2*order),
     :Ωact => Ωact
   )
-end  
+end
 
 ## Weak form
 const γg = 0.1
@@ -86,12 +85,12 @@ state_collection = EmbeddedCollection(model,φh) do _,_
   V = TestFESpace(Ωs.Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"])
   U = TrialFESpace(V,0.0)
   state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh;ls,adjoint_ls=ls)
-  (; 
+  (;
     :state_map => state_map,
     :J => StateParamIntegrandWithMeasure(J,state_map),
     :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,])
   )
-end  
+end
 pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,))
 
 ## Evolution Method
diff --git a/scripts/Embedded/Examples/2d_thermal_weird_background.jl b/scripts/Embedded/Examples/2d_thermal_weird_background.jl
index 2877aa93..a34dba98 100644
--- a/scripts/Embedded/Examples/2d_thermal_weird_background.jl
+++ b/scripts/Embedded/Examples/2d_thermal_weird_background.jl
@@ -52,11 +52,10 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-cos(8π*x[1])*cos(8π*x[2])-0.3,V_φ)
 Ωs = EmbeddedCollection(model,φh) do cutgeo,_
   Ωin = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL))
-  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo))
+  Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Γg   => Γg,
@@ -66,7 +65,7 @@ V_φ = TestFESpace(model,reffe_scalar)
     :dΓ   => Measure(Γ,2*order),
     :Ωact => Ωact
   )
-end  
+end
 
 ## Weak form
 const γg = 0.1
@@ -86,12 +85,12 @@ state_collection = EmbeddedCollection(model,φh) do _,_
   V = TestFESpace(Ωs.Ωact,reffe_scalar;dirichlet_tags=["Gamma_D"])
   U = TrialFESpace(V,0.0)
   state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh;ls,adjoint_ls=ls)
-  (; 
+  (;
     :state_map => state_map,
     :J => StateParamIntegrandWithMeasure(J,state_map),
     :C => map(Ci -> StateParamIntegrandWithMeasure(Ci,state_map),[Vol,])
   )
-end  
+end
 pcfs = EmbeddedPDEConstrainedFunctionals(state_collection;analytic_dC=(dVol,))
 
 ## Evolution Method
diff --git a/scripts/Embedded/Examples/FCM_2d_compliance_L.jl b/scripts/Embedded/Examples/FCM_2d_compliance_L.jl
index 5b11ea6e..a5e4f105 100644
--- a/scripts/Embedded/Examples/FCM_2d_compliance_L.jl
+++ b/scripts/Embedded/Examples/FCM_2d_compliance_L.jl
@@ -55,9 +55,8 @@ V_φ = TestFESpace(model,reffe_scalar)
   Ωout = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_OUT),V_φ)
   Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Ωout  => Ωout,
@@ -69,7 +68,7 @@ V_φ = TestFESpace(model,reffe_scalar)
     :dΓ   => Measure(Γ,2*order),
     :Ωact => Ωact
   )
-end  
+end
 
 ## Weak form
 function lame_parameters(E,ν)
@@ -84,7 +83,7 @@ E = 1.0
 g = VectorValue(0,-0.1)
 
 const ϵ = (λ + μ)*1e-3
-a(u,v,φ) = ∫(ε(v) ⊙ (σ ∘ ε(u)))Ωs.dΩin + 
+a(u,v,φ) = ∫(ε(v) ⊙ (σ ∘ ε(u)))Ωs.dΩin +
   # ∫(ϵ*(ε(v) ⊙ (σ ∘ ε(u))))Ωs.dΩout # Ersatz
   ∫(ϵ * ∇(v) ⊙ ∇(u))Ωs.dΩout # Pure diffusion
 l(v,φ) = ∫(v⋅g)dΓ_N
diff --git a/scripts/Embedded/Examples/FCM_2d_thermal.jl b/scripts/Embedded/Examples/FCM_2d_thermal.jl
index 1c9c9693..8dfd54a5 100644
--- a/scripts/Embedded/Examples/FCM_2d_thermal.jl
+++ b/scripts/Embedded/Examples/FCM_2d_thermal.jl
@@ -47,9 +47,8 @@ V_φ = TestFESpace(model,reffe_scalar)
   Ωout = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_OUT),V_φ)
   Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
   Γg = GhostSkeleton(cutgeo)
-  n_Γg = get_normal_vector(Γg)
   Ωact = Triangulation(cutgeo,ACTIVE)
-  (; 
+  (;
     :Ωin  => Ωin,
     :dΩin => Measure(Ωin,2*order),
     :Ωout  => Ωout,
diff --git a/scripts/Embedded/Examples/FCM_2d_thermal_MPI.jl b/scripts/Embedded/Examples/FCM_2d_thermal_MPI.jl
index b6b5386d..d50013fc 100644
--- a/scripts/Embedded/Examples/FCM_2d_thermal_MPI.jl
+++ b/scripts/Embedded/Examples/FCM_2d_thermal_MPI.jl
@@ -4,7 +4,10 @@ using GridapEmbedded, GridapEmbedded.LevelSetCutters
 
 using GridapTopOpt: StateParamIntegrandWithMeasure
 
-using GridapDistributed, PartitionedArrays, GridapPETSc
+using GridapDistributed, GridapPETSc, GridapSolvers,
+  PartitionedArrays, GridapTopOpt, SparseMatricesCSR
+
+using GridapSolvers: NewtonSolver
 
 function main(mesh_partition,distribute,el_size,path)
   ranks = distribute(LinearIndices((prod(mesh_partition),)))
@@ -53,9 +56,8 @@ function main(mesh_partition,distribute,el_size,path)
     Ωout = DifferentiableTriangulation(Triangulation(cutgeo,PHYSICAL_OUT),V_φ)
     Γ = DifferentiableTriangulation(EmbeddedBoundary(cutgeo),V_φ)
     Γg = GhostSkeleton(cutgeo)
-    n_Γg = get_normal_vector(Γg)
     Ωact = Triangulation(cutgeo,ACTIVE)
-    (; 
+    (;
       :Ωin  => Ωin,
       :dΩin => Measure(Ωin,2*order),
       :Ωout  => Ωout,
@@ -95,8 +97,13 @@ function main(mesh_partition,distribute,el_size,path)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map;analytic_dC=(dVol,))
 
   ## Evolution Method
-  evo = CutFEMEvolve(V_φ,Ωs,dΩ,h;max_steps)
-  reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=ArtificialViscosity(3.0))
+  evo = CutFEMEvolve(V_φ,Ωs,dΩ,h;max_steps,
+    ode_ls = solver,
+    assembler=SparseMatrixAssembler(Tm,Tv,V_φ,V_φ))
+  reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;
+    stabilisation_method=ArtificialViscosity(3.0),
+    assembler=SparseMatrixAssembler(Tm,Tv,V_φ,V_φ),
+    nls = NewtonSolver(LUSolver();maxiter=50,rtol=1.e-14,verbose=i_am_main(ranks)))
   ls_evo = UnfittedFEEvolution(evo,reinit)
   reinit!(ls_evo,φh)
 
@@ -118,9 +125,6 @@ function main(mesh_partition,distribute,el_size,path)
   optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;debug=true,
     γ,verbose=i_am_main(ranks),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) && i_am_main(ranks) && @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])
diff --git a/scripts/Embedded/MWEs/distributed.jl b/scripts/Embedded/MWEs/distributed.jl
index eb10f81d..552e0f78 100644
--- a/scripts/Embedded/MWEs/distributed.jl
+++ b/scripts/Embedded/MWEs/distributed.jl
@@ -5,6 +5,7 @@ using GridapDistributed, PartitionedArrays
 
 using GridapEmbedded
 using GridapEmbedded.LevelSetCutters
+using GridapEmbedded.Distributed
 using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData, Gridap.Adaptivity
 
 using GridapTopOpt: get_subfacet_normal_vector, get_ghost_normal_vector
@@ -39,6 +40,17 @@ cutgeo = cut(model,geo)
 Γ = EmbeddedBoundary(cutgeo)
 Λ = Skeleton(Γ)
 Σ = Boundary(Γ)
+Γg = GhostSkeleton(cutgeo)
+
+n_Γ = get_normal_vector(Γ);
+
+n_S_Λ = get_normal_vector(Λ);
+m_k_Λ = get_conormal_vector(Λ);
+∇ˢφ_Λ = Operation(abs)(n_S_Λ ⋅ ∇(φh).plus);
+
+n_S_Σ = get_normal_vector(Σ);
+m_k_Σ = get_conormal_vector(Σ);
+∇ˢφ_Σ = Operation(abs)(n_S_Σ ⋅ ∇(φh));
 
 ############################################################################################
 
diff --git a/scripts/Embedded/MWEs/distributed_differentiable_trian.jl b/scripts/Embedded/MWEs/distributed_differentiable_trian.jl
index fb1cfd5c..e5bd5f9f 100644
--- a/scripts/Embedded/MWEs/distributed_differentiable_trian.jl
+++ b/scripts/Embedded/MWEs/distributed_differentiable_trian.jl
@@ -1,3 +1,5 @@
+using Test
+
 using GridapTopOpt
 using Gridap
 
@@ -10,14 +12,16 @@ using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData, Gridap.Adaptivity
 using GridapTopOpt: get_subfacet_normal_vector, get_ghost_normal_vector
 using GridapTopOpt: get_conormal_vector, get_tangent_vector
 
+using GridapDistributed: DistributedTriangulation, DistributedDomainContribution
+
 order = 1
 n = 16
 
 parts = (2,2)
 ranks = DebugArray(LinearIndices((prod(parts),)))
 
-_model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),(n,n))
-# _model = CartesianDiscreteModel((0,1,0,1),(n,n))
+# _model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),(n,n))
+_model = CartesianDiscreteModel((0,1,0,1),(n,n))
 
 base_model = UnstructuredDiscreteModel(_model)
 ref_model = refine(base_model, refinement_method = "barycentric")
@@ -32,8 +36,10 @@ fh = interpolate(x->cos(x[1]*x[2]),V_φ)
 geo = DiscreteGeometry(φh,model)
 cutgeo = cut(model,geo)
 
-Ωin = Triangulation(cutgeo,PHYSICAL_IN)
-Ω_AD = DifferentiableTriangulation(Ωin,V_φ)
+# A.1) Volume integral
+
+Ω = Triangulation(cutgeo,PHYSICAL_IN)
+Ω_AD = DifferentiableTriangulation(Ω,V_φ)
 dΩ = Measure(Ω_AD,2*order)
 
 Γ = EmbeddedBoundary(cutgeo)
@@ -46,4 +52,61 @@ dJ_bulk_AD_vec = assemble_vector(dJ_bulk_AD,V_φ)
 dJ_bulk_exact(q) = ∫(-fh*q/(norm ∘ (∇(φh))))dΓ
 dJ_bulk_exact_vec = assemble_vector(dJ_bulk_exact,V_φ)
 
-norm(dJ_bulk_AD_vec - dJ_bulk_exact_vec) < 1e-10
\ No newline at end of file
+@test norm(dJ_bulk_AD_vec - dJ_bulk_exact_vec) < 1e-10
+
+# A.2) Volume integral
+
+g(fh) = ∇(fh)⋅∇(fh)
+J_bulk2(φ) = ∫(g(fh))dΩ
+dJ_bulk_AD2 = gradient(J_bulk2,φh)
+dJ_bulk_AD_vec2 = assemble_vector(dJ_bulk_AD2,V_φ)
+
+dJ_bulk_exact2(q) = ∫(-g(fh)*q/(norm ∘ (∇(φh))))dΓ
+dJ_bulk_exact_vec2 = assemble_vector(dJ_bulk_exact2,V_φ)
+
+@test norm(dJ_bulk_AD_vec2 - dJ_bulk_exact_vec2) < 1e-10
+
+# B.1) Facet integral
+
+Γ = EmbeddedBoundary(cutgeo)
+Γ_AD = DifferentiableTriangulation(Γ,V_φ)
+Λ = Skeleton(Γ)
+Σ = Boundary(Γ)
+
+dΓ = Measure(Γ,2*order);
+dΛ = Measure(Λ,2*order);
+dΣ = Measure(Σ,2*order);
+
+n_Γ = get_normal_vector(Γ);
+
+n_S_Λ = get_normal_vector(Λ);
+m_k_Λ = get_conormal_vector(Λ);
+∇ˢφ_Λ = Operation(abs)(n_S_Λ ⋅ ∇(φh).plus);
+
+n_S_Σ = get_normal_vector(Σ);
+m_k_Σ = get_conormal_vector(Σ);
+∇ˢφ_Σ = Operation(abs)(n_S_Σ ⋅ ∇(φh));
+
+dΓ_AD = Measure(Γ_AD,2*order)
+J_int(φ) = ∫(fh)dΓ_AD
+dJ_int_AD = gradient(J_int,φh)
+dJ_int_AD_vec = assemble_vector(dJ_int_AD,V_φ)
+
+dJ_int_exact(w) = ∫((-n_Γ⋅∇(fh))*w/(norm ∘ (∇(φh))))dΓ +
+                  ∫(-n_S_Λ ⋅ (jump(fh*m_k_Λ) * mean(w) / ∇ˢφ_Λ))dΛ +
+                  ∫(-n_S_Σ ⋅ (fh*m_k_Σ * w / ∇ˢφ_Σ))dΣ
+dJ_int_exact_vec = assemble_vector(dJ_int_exact,V_φ)
+
+@test norm(dJ_int_AD_vec - dJ_int_exact_vec) < 1e-10
+
+# B.2) Facet integral
+J_int2(φ) = ∫(g(fh))dΓ_AD
+dJ_int_AD2 = gradient(J_int2,φh)
+dJ_int_AD_vec2 = assemble_vector(dJ_int_AD2,V_φ)
+
+dJ_int_exact2(w) = ∫((-n_Γ⋅∇(g(fh)))*w/(norm ∘ (∇(φh))))dΓ +
+                  ∫(-n_S_Λ ⋅ (jump(g(fh)*m_k_Λ) * mean(w) / ∇ˢφ_Λ))dΛ +
+                  ∫(-n_S_Σ ⋅ (g(fh)*m_k_Σ * w / ∇ˢφ_Σ))dΣ
+dJ_int_exact_vec2 = assemble_vector(dJ_int_exact2,V_φ)
+
+@test norm(dJ_int_AD_vec2 - dJ_int_exact_vec2) < 1e-10
\ No newline at end of file
diff --git a/src/Embedded/DifferentiableTriangulations.jl b/src/Embedded/DifferentiableTriangulations.jl
index e47ee464..a831e493 100644
--- a/src/Embedded/DifferentiableTriangulations.jl
+++ b/src/Embedded/DifferentiableTriangulations.jl
@@ -4,11 +4,11 @@
 """
   mutable struct DifferentiableTriangulation{Dc,Dp,A<:Triangulation{Dc,Dp}} <: Triangulation{Dc,Dp}
 
-A DifferentiableTriangulation is a wrapper around an embedded triangulation 
-(i.e SubCellTriangulation or SubFacetTriangulation) implementing all the necessary 
+A DifferentiableTriangulation is a wrapper around an embedded triangulation
+(i.e SubCellTriangulation or SubFacetTriangulation) implementing all the necessary
 methods to compute derivatives wrt deformations of the embedded mesh.
 
-To do so, it propagates dual numbers into the geometric maps mapping cut subcells/subfacets 
+To do so, it propagates dual numbers into the geometric maps mapping cut subcells/subfacets
 to the background mesh.
 """
 mutable struct DifferentiableTriangulation{Dc,Dp,A,B} <: Triangulation{Dc,Dp}
@@ -117,7 +117,7 @@ function Geometry.get_glue(ttrian::DifferentiableTriangulation,val::Val{D}) wher
     return glue
   end
 
-  # New reference maps 
+  # New reference maps
   c = ttrian.caches
   cell_values = ttrian.cell_values
   cell_to_rcoords = lazy_map(
@@ -160,6 +160,15 @@ for tdomain in (:ReferenceDomain,:PhysicalDomain)
   end
 end
 
+# function CellData.change_domain(
+#   a::CellField,target_trian::DifferentiableTriangulation{Dc,Dp,A},target_domain::ReferenceDomain
+# ) where {Dc,Dp,A <: Union{SubCellTriangulation,SubFacetTriangulation}}
+#   change_domain(a,get_triangulation(a),DomainStyle(a),target_trian,target_domain)
+#   @assert is_change_possible(strian,ttrian)
+#   b = change_domain(a,$(tdomain)())
+#   return CellData.similar_cell_field(a,CellData.get_data(b),ttrian,$(tdomain)())
+# end
+
 function FESpaces.get_cell_fe_data(fun,f,ttrian::DifferentiableTriangulation)
   FESpaces.get_cell_fe_data(fun,f,ttrian.trian)
 end
@@ -220,9 +229,9 @@ end
 """
     precompute_cut_edge_ids(rcoords,bg_rcoords,edge_list)
 
-Given 
+Given
   - `rcoords`: the node ref coordinates of the cut subcell/subfacet,
-  - `bg_rcoords`: the node ref coordinates of the background cell containing it, 
+  - `bg_rcoords`: the node ref coordinates of the background cell containing it,
   - `edge_list`: the list of nodes defining each edge of the background cell,
 
 this function returns a vector that for each node of the cut subcell/subfacet contains
@@ -368,10 +377,10 @@ end
 # AppendedTriangulation
 #
 # When cutting an embedded domain, we will usually end up with an AppendedTriangulation
-# containing 
+# containing
 #   a) a regular triangulation with the IN/OUT cells
 #   b) a SubCell/SubFacetTriangulation with the CUT cells
-# We only need to propagate the dual numbers to the CUT cells, which is what the 
+# We only need to propagate the dual numbers to the CUT cells, which is what the
 # following implementation does:
 
 DifferentiableTriangulation(trian::Triangulation,fe_space) = trian
@@ -413,6 +422,32 @@ function FESpaces._compute_cell_ids(uh,ttrian::AppendedTriangulation)
   lazy_append(ids_a,ids_b)
 end
 
+#### DistributedTriangulations
+
+function DifferentiableTriangulation(trian::DistributedTriangulation,fe_space)
+  bg_model = get_background_model(trian)
+  trians = map(DifferentiableTriangulation,local_views(trian),local_views(fe_space))
+  return DistributedTriangulation(trians,bg_model)
+end
+
+# TODO: Dispatch more appropriately on trian?
+function FESpaces._change_argument(op,f,trian,uh::DistributedCellField)
+  spaces = map(get_fe_space,local_views(uh))
+  function g(cell_u)
+    cf = DistributedCellField(
+      map(CellField,spaces,cell_u),
+      get_triangulation(uh),
+    )
+    map(update_trian!,local_views(trian),local_views(spaces),local_views(cf))
+    cell_grad = f(cf)
+    map(local_views(trian),local_views(spaces)) do Ω, V
+      update_trian!(Ω,V,nothing)
+    end
+    map(get_contribution,local_views(cell_grad),trian)
+  end
+  g
+end
+
 # TriangulationView
 function DifferentiableTriangulation(
   trian::Geometry.TriangulationView,
@@ -425,18 +460,4 @@ end
 function update_trian!(trian::Geometry.TriangulationView,U,φh)
   update_trian!(trian.parent,U,φh)
   return trian
-end
-
-function FESpaces._change_argument(
-  op,f,trian::Geometry.TriangulationView,uh::SingleFieldFEFunction
-)
-  U = get_fe_space(uh)
-  function g(cell_u)
-    cf = CellField(U,cell_u)
-    update_trian!(trian,U,cf)
-    cell_grad = f(cf)
-    update_trian!(trian,U,nothing) # TODO: experimental
-    get_contribution(cell_grad,trian)
-  end
-  g
 end
\ No newline at end of file
diff --git a/src/Embedded/SubFacetBoundaryTriangulations.jl b/src/Embedded/SubFacetBoundaryTriangulations.jl
index e89d3417..09834263 100644
--- a/src/Embedded/SubFacetBoundaryTriangulations.jl
+++ b/src/Embedded/SubFacetBoundaryTriangulations.jl
@@ -170,17 +170,17 @@ Triangulation containing the interfaces between subfacets. We always have dimens
   - Df = Dc-1 :: Dimension of the cut subfacets
   - Di = Dc-2 :: Dimension of the subfacet interfaces
 
-Description of the different components: 
+Description of the different components:
 
 - `face_trian` :: Original SubFacetTriangulation, built on top of the background mesh.
-- `face_model` :: Subfacet model. Active model for `face_trian`. 
+- `face_model` :: Subfacet model. Active model for `face_trian`.
 - `face_boundary` :: Triangulation of the interfaces between subfacets. It is glued to the `face_model`.
-- `cell_boundary` :: Conceptually the same as `face_boundary`, but it is glued to the 
-                     background mesh cells. Created as a CompositeTriangulation between 
+- `cell_boundary` :: Conceptually the same as `face_boundary`, but it is glued to the
+                     background mesh cells. Created as a CompositeTriangulation between
                       `face_trian` and `face_boundary`.
-- `ghost_boundary` :: Triangulation of the background facets that contain each interface. 
+- `ghost_boundary` :: Triangulation of the background facets that contain each interface.
 
-The "real" triangulation is `cell_boundary`, but we require the other triangulations to 
+The "real" triangulation is `cell_boundary`, but we require the other triangulations to
 perform complex changes of domain.
 """
 struct SubFacetBoundaryTriangulation{Di,Df,Dp} <: Triangulation{Di,Dp}
@@ -279,8 +279,8 @@ function get_interface_sign(trian::SubFacetBoundaryTriangulation)
 end
 
 # TODO: This is only valid when dealing with linear meshes (where normals are constant over facets).
-# If we wanted to use higher-order meshes, we would need to generate the geometric map 
-# going from the facets to the interfaces. 
+# If we wanted to use higher-order meshes, we would need to generate the geometric map
+# going from the facets to the interfaces.
 # However, having a high-order background mesh seems quite silly.
 function get_ghost_facet_normal(
   itrian::CompositeTriangulation{Di,Dc}, # Interface -> BG Cell
@@ -305,8 +305,8 @@ function get_subfacet_facet_normal(
 end
 
 """
-  There is still something sweaty about this... 
-  Why do we apply the sign change but at the same time call `get_edge_tangents` in the 
+  There is still something sweaty about this...
+  Why do we apply the sign change but at the same time call `get_edge_tangents` in the
   creation of conormal vectors in 3D? It's like we are cancelling the sign change...
   There is more to think about here.
 """
@@ -362,6 +362,8 @@ function CellData.get_normal_vector(trian::SubFacetBoundaryTriangulation{Di}) wh
   return n_S * isign
 end
 
+CellData.get_facet_normal(trian::SubFacetBoundaryTriangulation) = get_data(get_normal_vector(trian))
+
 # Tangent vector to the cut interface, t_S = n_S x n_k
 function get_tangent_vector(trian::SubFacetBoundaryTriangulation{Di}) where {Di}
   @notimplementedif Di != 1
@@ -386,9 +388,13 @@ function get_conormal_vector(trian::SubFacetBoundaryTriangulation{Di}) where {Di
   return m_k * isign
 end
 
-# SubFacetSkeletonTriangulation
-
-const SubFacetSkeletonTriangulation{Di,Df,Dp} = SkeletonTriangulation{Di,Dp,SubFacetBoundaryTriangulation{Di,Df,Dp}}
+# SubFacetSkeletonTriangulation & SubFacetBoundaryTriangulationView
+const SubFacetBoundaryTriangulationView{Di,Df,Dp} = TriangulationView{Di,Dp,SubFacetBoundaryTriangulation{Di,Df,Dp}}
+const SubFacetSkeletonTriangulation{Di,Df,Dp} = SkeletonTriangulation{Di,Dp,<:Union{
+  SubFacetBoundaryTriangulation{Di,Df,Dp},
+  SubFacetBoundaryTriangulationView{Di,Df,Dp}
+  }
+}
 
 function Geometry.SkeletonTriangulation(face_trian::SubFacetTriangulation)
   bgmodel = get_background_model(face_trian)
@@ -398,7 +404,7 @@ function Geometry.SkeletonTriangulation(face_trian::SubFacetTriangulation)
   face_skeleton = SkeletonTriangulation(face_model,ghost_mask)
   cell_skeleton = CompositeTriangulation(face_trian,face_skeleton)
   ghost_skeleton = generate_ghost_trian(cell_skeleton,bgmodel)
-  
+
   ctrian_plus = CompositeTriangulation(face_trian,face_skeleton.plus)
   ctrian_minus = CompositeTriangulation(face_trian,face_skeleton.minus)
   isign_plus = get_interface_sign(ctrian_plus,face_trian,ghost_skeleton.plus)
@@ -433,8 +439,18 @@ for func in (:(Gridap.get_normal_vector),:get_tangent_vector)
   end
 end
 
+for func in (:get_tangent_vector,:get_subfacet_normal_vector,:get_ghost_normal_vector,:get_conormal_vector)
+  @eval begin
+    function $func(trian::SubFacetBoundaryTriangulationView)
+      data = CellData.get_data($func(trian.parent))
+      restricted_data = Geometry.restrict(data,trian.cell_to_parent_cell)
+      return GenericCellField(restricted_data,trian,ReferenceDomain())
+    end
+  end
+end
+
 # Distributed
-# Until we merge, we need changes in 
+# Until we merge, we need changes in
 #   - GridapDistributed#master
 #   - GridapEmbedded#distributed
 
diff --git a/src/GridapTopOpt.jl b/src/GridapTopOpt.jl
index ffd5e625..3cdb41e9 100644
--- a/src/GridapTopOpt.jl
+++ b/src/GridapTopOpt.jl
@@ -14,7 +14,7 @@ using Gridap.Helpers, Gridap.Algebra, Gridap.TensorValues
 using Gridap.Geometry, Gridap.CellData, Gridap.Fields, Gridap.Arrays
 using Gridap.ReferenceFEs, Gridap.FESpaces,  Gridap.MultiField, Gridap.Polynomials
 
-using Gridap.Geometry: get_faces, num_nodes
+using Gridap.Geometry: get_faces, num_nodes, TriangulationView
 using Gridap.FESpaces: get_assembly_strategy
 using Gridap.ODEs: ODESolver
 using Gridap: writevtk
diff --git a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl
index b9352a2b..df0862ba 100644
--- a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl
+++ b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl
@@ -1,17 +1,17 @@
 """
     mutable struct CutFEMEvolve{V,M} <: Evolver
 
-CutFEM method for level-set evolution based on method developed by 
+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 
+- `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 
+  max steps `max_steps`, and background mesh skeleton parameters
 - `cache`: Cache for evolver, initially `nothing`.
 
 # Note
@@ -33,7 +33,7 @@ mutable struct CutFEMEvolve{A,B,C} <: Evolver
       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_φ)) where {A,B}
-    Γg = SkeletonTriangulation(get_triangulation(dΩ_bg.quad))
+    Γg = SkeletonTriangulation(get_triangulation(V_φ))
     dΓg = Measure(Γg,2get_order(V_φ))
     n_Γg = get_normal_vector(Γg)
     params = (;γg,h,max_steps,dΓg,n_Γg)
@@ -84,7 +84,7 @@ function solve!(s::CutFEMEvolve,φh,velh,γ,cache::Nothing)
   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,φh)
-  
+
   # March
   march = Base.iterate(ode_sol)
   data, state = march
@@ -110,13 +110,13 @@ function solve!(s::CutFEMEvolve,φh,velh,γ,cache)
 
   ## Update state
   # `get_transient_operator` re-creates the entire TransientLinearFEOperator wrapper.
-  #   We do this so that the first iterate of ODESolution always recomputes the 
+  #   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 
+  # 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 
+  #   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)
diff --git a/test/seq/EmbeddedDifferentiationTests.jl b/test/seq/EmbeddedDifferentiationTests.jl
index a428a722..9b2eb0b7 100644
--- a/test/seq/EmbeddedDifferentiationTests.jl
+++ b/test/seq/EmbeddedDifferentiationTests.jl
@@ -65,7 +65,7 @@ function main(
   # A.1) Volume integral
 
   Ω = Triangulation(cutgeo,PHYSICAL_IN)
-  Ω_AD = DifferentiableTriangulation(Ω)
+  Ω_AD = DifferentiableTriangulation(Ω,V_φ)
   dΩ = Measure(Ω_AD,2*order)
 
   Γ = EmbeddedBoundary(cutgeo)
@@ -95,7 +95,7 @@ function main(
   # B.1) Facet integral
 
   Γ = EmbeddedBoundary(cutgeo)
-  Γ_AD = DifferentiableTriangulation(Γ)
+  Γ_AD = DifferentiableTriangulation(Γ,V_φ)
   Λ = Skeleton(Γ)
   Σ = Boundary(Γ)
 
@@ -118,7 +118,7 @@ function main(
   dJ_int_AD = gradient(J_int,φh)
   dJ_int_AD_vec = assemble_vector(dJ_int_AD,V_φ)
 
-  dJ_int_exact(w) = ∫((-n_Γ⋅∇(fh))*w/(norm ∘ (∇(φh))))dΓ + 
+  dJ_int_exact(w) = ∫((-n_Γ⋅∇(fh))*w/(norm ∘ (∇(φh))))dΓ +
                     ∫(-n_S_Λ ⋅ (jump(fh*m_k_Λ) * mean(w) / ∇ˢφ_Λ))dΛ +
                     ∫(-n_S_Σ ⋅ (fh*m_k_Σ * w / ∇ˢφ_Σ))dΣ
   dJ_int_exact_vec = assemble_vector(dJ_int_exact,V_φ)
@@ -131,7 +131,7 @@ function main(
   dJ_int_AD2 = gradient(J_int2,φh)
   dJ_int_AD_vec2 = assemble_vector(dJ_int_AD2,V_φ)
 
-  dJ_int_exact2(w) = ∫((-n_Γ⋅∇(g(fh)))*w/(norm ∘ (∇(φh))))dΓ + 
+  dJ_int_exact2(w) = ∫((-n_Γ⋅∇(g(fh)))*w/(norm ∘ (∇(φh))))dΓ +
                     ∫(-n_S_Λ ⋅ (jump(g(fh)*m_k_Λ) * mean(w) / ∇ˢφ_Λ))dΛ +
                     ∫(-n_S_Σ ⋅ (g(fh)*m_k_Σ * w / ∇ˢφ_Σ))dΣ
   dJ_int_exact_vec2 = assemble_vector(dJ_int_exact2,V_φ)
@@ -155,7 +155,7 @@ function main(
       ];
       append = false
     )
-  
+
     writevtk(
       Ω, "$(path)_omega"; append = false
     )