diff --git a/scripts/Embedded/MWEs/jordi_embedded.jl b/scripts/Embedded/MWEs/jordi_embedded.jl
index e1ab5ab..9213c21 100644
--- a/scripts/Embedded/MWEs/jordi_embedded.jl
+++ b/scripts/Embedded/MWEs/jordi_embedded.jl
@@ -1,5 +1,4 @@
 using GridapTopOpt
-
 using Gridap
 
 using GridapEmbedded
@@ -23,12 +22,12 @@ V_φ = TestFESpace(model,reffe_scalar)
 # ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,1/(10n),0)
 
 # φh = interpolate(x->max(abs(x[1]-0.5),abs(x[2]-0.5))-0.25,V_φ) # Square
- φh = interpolate(x->((x[1]-0.5)^N+(x[2]-0.5)^N)^(1/N)-0.25,V_φ) # Curved corner example
+#φh = interpolate(x->((x[1]-0.5)^N+(x[2]-0.5)^N)^(1/N)-0.25,V_φ) # Curved corner example
 #φh = interpolate(x->abs(x[1]-0.5)+abs(x[2]-0.5)-0.25-0/n/10,V_φ) # Diamond
 # φh = interpolate(x->(1+0.25)abs(x[1]-0.5)+0abs(x[2]-0.5)-0.25,V_φ) # Straight lines with scaling
 # φh = interpolate(x->abs(x[1]-0.5)+0abs(x[2]-0.5)-0.25/(1+0.25),V_φ) # Straight lines without scaling
 # φh = interpolate(x->tan(-pi/4)*(x[1]-0.5)+(x[2]-0.5),V_φ) # Angled interface
-# φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.303,V_φ) # Circle
+φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)-0.303,V_φ) # Circle
 
 x_φ = get_free_dof_values(φh)
 # @assert ~any(isapprox(0.0;atol=10^-10),x_φ)
@@ -37,9 +36,6 @@ if any(isapprox(0.0;atol=10^-10),x_φ)
   idx = findall(isapprox(0.0;atol=10^-10),x_φ)
   x_φ[idx] .+= 10eps()
 end
-any(x -> x < 0,x_φ)
-any(x -> x > 0,x_φ)
-# reinit!(ls_evo,φh,0.5)
 
 geo = DiscreteGeometry(φh,model)
 cutgeo = cut(model,geo)
@@ -50,8 +46,6 @@ cutgeo = cut(model,geo)
 diff_Ωin = DifferentiableTriangulation(Ωin)
 diff_Ωout = DifferentiableTriangulation(Ωout)
 
-oh = interpolate(1.0,V_φ)
-
 dΩin = Measure(diff_Ωin,3*order)
 j_in(φ) = ∫(1)dΩin
 dj_in = gradient(j_in,φh)
@@ -59,7 +53,7 @@ dj_vec_in = assemble_vector(dj_in,V_φ)
 norm(dj_vec_in)
 
 dΩout = Measure(diff_Ωout,3*order)
-j_out(φ) = ∫(oh)dΩout
+j_out(φ) = ∫(1)dΩout
 dj_out = gradient(j_out,φh)
 dj_vec_out = -assemble_vector(dj_out,V_φ)
 norm(dj_vec_out)
@@ -78,6 +72,21 @@ dΓ2 = Measure(diff_Γ,3*order)
 jΓ(φ) = ∫(1)dΓ2
 djΓ = gradient(jΓ,φh)
 
+############################################################################################
+
+include("../SubFacetSkeletons.jl")
+
+Γ = EmbeddedBoundary(cutgeo)
+Λ = Skeleton(Γ)
+
+n_∂Ω = get_subfacet_normal_vector(Λ)
+n_k = get_ghost_normal_vector(Λ)
+#t_S = get_tangent_vector(Λ)
+n_S = get_normal_vector(Λ)
+m_k = get_conormal_vector(Λ)
+
+############################################################################################
+
 writevtk(
   Ω,"results/test",
   cellfields=["φh"=>φh,"∇φh"=>∇(φh),"dj_in"=>FEFunction(V_φ,dj_vec_in),"dj_expected"=>FEFunction(V_φ,dj_exp_vec),"dj_out"=>FEFunction(V_φ,dj_vec_out)],
@@ -92,5 +101,4 @@ writevtk(
 )
 writevtk(
   Γ,"results/gamma"
-)
-
+)
\ No newline at end of file
diff --git a/scripts/Embedded/SubFacetSkeletonTriangulation.jl b/scripts/Embedded/SubFacetSkeletonTriangulation.jl
index 3104451..b6c26cc 100644
--- a/scripts/Embedded/SubFacetSkeletonTriangulation.jl
+++ b/scripts/Embedded/SubFacetSkeletonTriangulation.jl
@@ -1,6 +1,6 @@
 using Gridap
 using Gridap.Geometry, Gridap.Adaptivity, Gridap.Algebra, Gridap.Arrays, Gridap.FESpaces
-  Gridap.CellData, Gridap.Fields, Gridap.Helpers, Gridap.ReferenceFEs, Gridap.Polynomials
+using Gridap.CellData, Gridap.Fields, Gridap.Helpers, Gridap.ReferenceFEs, Gridap.Polynomials
 
 using GridapEmbedded
 using GridapEmbedded.Interfaces
@@ -18,12 +18,14 @@ function SubFacetSkeletonTriangulation(trian::SubFacetTriangulation{Dc,Dp}) wher
   if Dp == 2
     Γ_facet_to_points = map(Reindex(subfacets.point_to_coords),subfacets.facet_to_points)
     Γ_pts = unique(x -> round(x;sigdigits=12),reduce(vcat,Γ_facet_to_points))
-    Γ_facet_to_pts = convert(Vector{Vector{Int32}},map(v->indexin(round.(v;sigdigits=12),
-      round.(Γ_pts;sigdigits=12)),Γ_facet_to_points))
+    Γ_facet_to_pts = Table(convert(
+      Vector{Vector{Int32}},
+      map(v->indexin(round.(v;sigdigits=12),round.(Γ_pts;sigdigits=12)),Γ_facet_to_points)
+    ))
 
-    grid_top = GridTopology(Γ.subgrid,Table(Γ_facet_to_pts),IdentityVector(length(Γ_pts)))
+    grid_top = GridTopology(Γ.subgrid,Γ_facet_to_pts,IdentityVector(length(Γ_pts)))
     grid = UnstructuredGrid(
-      Γ_pts,Table(Γ_facet_to_pts),
+      Γ_pts,Γ_facet_to_pts,
       Γ.subgrid.reffes,
       Γ.subgrid.cell_types,
       Γ.subgrid.orientation_style,
@@ -96,7 +98,7 @@ function orient(a::VectorValue{2,T},b::VectorValue{2,T}) where T
 end
 
 order = 1
-n = 51
+n = 10
 N = 4
 _model = CartesianDiscreteModel((0,1,0,1),(n,n))
 cd = Gridap.Geometry.get_cartesian_descriptor(_model)
@@ -139,6 +141,9 @@ dΓ = Measure(Γ,2*order)
 dΓs = Measure(Γs,2)
 fh = interpolate(x->1,V_φ)
 
+# Jordi: I think this is equivalent to the _normal_vector & _orthogonal_vector methods 
+# in GridapEmbedded/LevelSetCutters/LookupTables - line 333/350
+# There it is generalised to 3D and 4D.
 _2d_cross(n) = VectorValue(n[2],-n[1]);
 
 Γg_to_Γs = Γg_to_Γs_perm(Γs,Γg)
@@ -215,21 +220,25 @@ writevtk(
   "n_∂Ω_plus"=>n_∂Ω_plus,"n_∂Ω_minus"=>n_∂Ω_minus,
   "ns"=>nˢ,
   "∇φh_Γs_plus"=>∇φh_Γs.plus,"∇φh_Γs_minus"=>∇φh_Γs.minus,
-  "m.minus"=>m.minus,"m.plus"=>m.plus]
+  "m.minus"=>m.minus,"m.plus"=>m.plus];
+  append=false
 )
 bgcell_to_inoutcut = compute_bgcell_to_inoutcut(model,geo)
-  writevtk(
+writevtk(
   Ω,
   "results/Background",
   cellfields=["φh"=>∇(φh)]#,"dj"=>FEFunction(V_φ,djΓ_exp_vec),"dj_AD"=>FEFunction(V_φ,djΓ_vec_out)]
   ,
-  celldata=["inoutcut"=>bgcell_to_inoutcut]
-  )
+  celldata=["inoutcut"=>bgcell_to_inoutcut];
+  append=false
+)
 writevtk(
   Γ,
-  "results/Boundary"
+  "results/Boundary";
+  append=false
 )
 writevtk(
   Γg,
-  "results/GhostSkel"
-)
\ No newline at end of file
+  "results/GhostSkel";
+  append=false
+)
diff --git a/scripts/Embedded/SubFacetSkeletons.jl b/scripts/Embedded/SubFacetSkeletons.jl
index f3a8e22..7f943f3 100644
--- a/scripts/Embedded/SubFacetSkeletons.jl
+++ b/scripts/Embedded/SubFacetSkeletons.jl
@@ -62,7 +62,7 @@ function generate_ghost_trian(
   cell_to_face = get_faces(topo,Dc,Dc-1)
 
   n_bgfaces = num_faces(bgmodel,1)
-  n_faces = num_cells(cell_skeleton)
+  n_faces = num_cells(trian)
   ghost_faces = zeros(Int32,n_faces)
   p_lcell = ones(Int8,n_bgfaces)
   m_lcell = ones(Int8,n_bgfaces)
@@ -81,9 +81,10 @@ function generate_ghost_trian(
   minus = BoundaryTriangulation(bgmodel,ghost_faces,m_lcell)
   return SkeletonTriangulation(plus,minus)
 end
+
 struct SubFacetSkeletonTriangulation{Di,Df,Dp} <: Triangulation{Di,Dp}
   cell_skeleton::CompositeTriangulation{Di,Dp} # Interface -> BG Cell pair
-  face_skeleton::SkeletonTriangulation{Df,Dp}  # Interface -> Cut Facet pair
+  face_skeleton::SkeletonTriangulation{Di,Dp}  # Interface -> Cut Facet pair
   ghost_skeleton::SkeletonTriangulation{Df,Dp} # Ghost Facet -> BG Cell pair
   face_trian::SubFacetTriangulation{Df,Dp}     # Cut Facet -> BG Cell
   face_model::UnstructuredDiscreteModel{Df,Dp}
@@ -128,9 +129,9 @@ end
 # Normal vector to the ghost facets, n_k
 function get_ghost_normal_vector(trian::SubFacetSkeletonTriangulation)
   n_out = get_normal_vector(trian.ghost_skeleton)
-  n_in  = (-1) * n_out
-  plus  =  GenericCellField(CellData.get_data(n_in.plus),trian.cell_skeleton.plus,ReferenceDomain())
-  minus = GenericCellField(CellData.get_data(n_in.minus),trian.cell_skeleton.minus,ReferenceDomain())
+  n_in  = SkeletonPair((-1) * n_out.plus, (-1) * n_out.minus)
+  plus  = GenericCellField(CellData.get_data(n_in.plus),trian,ReferenceDomain())
+  minus = GenericCellField(CellData.get_data(n_in.minus),trian,ReferenceDomain())
   return SkeletonPair(plus,minus)
 end
 
@@ -144,15 +145,14 @@ function get_tangent_vector(trian::SkeletonTriangulation{1})
     t = c[2] - c[1]
     return t/norm(t)
   end
-  edge_coords = get_cell_coordinates(trian)
-  data = lazy_map(t,CellData.get_data(edge_coords))
+  data = lazy_map(t,get_cell_coordinates(trian))
   plus = GenericCellField(data,trian.plus,PhysicalDomain())
   minus = GenericCellField(data,trian.minus,PhysicalDomain())
   return SkeletonPair(plus,minus)
 end
 
 # Normal vector to the cut interface, n_S
-function get_normal_vector(trian::SubFacetSkeletonTriangulation{Di}) where {Di}
+function CellData.get_normal_vector(trian::SubFacetSkeletonTriangulation{Di}) where {Di}
   if Di == 0
     n_S = get_tangent_vector(trian.ghost_skeleton)
   elseif Di == 1
@@ -175,7 +175,7 @@ end
 
 # Conormal vectors, m_k = t_S x n_∂Ω
 function get_conormal_vector(trian::SubFacetSkeletonTriangulation{Di}) where {Di}
-  op = Operation(GridapEmbedded._normal_vector)
+  op = Operation(GridapEmbedded.LevelSetCutters._normal_vector)
   n_∂Ω = get_subfacet_normal_vector(trian)
   if Di == 0 # 2D 
     m_k = op(n_∂Ω)
@@ -187,3 +187,9 @@ function get_conormal_vector(trian::SubFacetSkeletonTriangulation{Di}) where {Di
   end
   return m_k
 end
+
+function Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:CellField})
+  plus = k(a.plus)
+  minus = k(a.minus)
+  SkeletonPair(plus,minus)
+end
\ No newline at end of file