diff --git a/scripts/Embedded/MWEs/jordi_embedded.jl b/scripts/Embedded/MWEs/jordi_embedded.jl
index 68e987f..4071b65 100644
--- a/scripts/Embedded/MWEs/jordi_embedded.jl
+++ b/scripts/Embedded/MWEs/jordi_embedded.jl
@@ -3,7 +3,7 @@ using Gridap
 
 using GridapEmbedded
 using GridapEmbedded.LevelSetCutters
-using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData
+using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData, Gridap.Adaptivity
 import Gridap.Geometry: get_node_coordinates, collect1d
 
 include("../differentiable_trians.jl")
@@ -12,14 +12,18 @@ order = 1
 n = 2
 N = 8
 
-model = CartesianDiscreteModel((0,1,0,1),(n,n))
-model = simplexify(model)
+# _model = CartesianDiscreteModel((0,1,0,1),(n,n))
+_model = CartesianDiscreteModel((0,1,0,1,0,1),(n,n,n))
+cd = Gridap.Geometry.get_cartesian_descriptor(_model)
+base_model = UnstructuredDiscreteModel(_model)
+ref_model = refine(base_model, refinement_method = "barycentric")
+model = ref_model.model
+# model = simplexify(model)
 Ω = Triangulation(model)
 dΩ = Measure(Ω,2*order)
 
 reffe_scalar = ReferenceFE(lagrangian,Float64,order)
 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
@@ -27,7 +31,11 @@ V_φ = TestFESpace(model,reffe_scalar)
 # φ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
+φh = interpolate(x->sqrt((x[1]-0.5)^2+(x[2]-0.5)^2+(x[3]-0.5)^2)-0.303,V_φ) # Sphere
+
+fh = interpolate(x->cos(x[1]*x[2]),V_φ)
 
 x_φ = get_free_dof_values(φh)
 # @assert ~any(isapprox(0.0;atol=10^-10),x_φ)
@@ -47,31 +55,26 @@ diff_Ωin = DifferentiableTriangulation(Ωin)
 diff_Ωout = DifferentiableTriangulation(Ωout)
 
 dΩin = Measure(diff_Ωin,3*order)
-j_in(φ) = ∫(1)dΩin
+j_in(φ) = ∫(fh)dΩin
 dj_in = gradient(j_in,φh)
 dj_vec_in = assemble_vector(dj_in,V_φ)
 norm(dj_vec_in)
 
 dΩout = Measure(diff_Ωout,3*order)
-j_out(φ) = ∫(1)dΩout
+j_out(φ) = ∫(fh)dΩout
 dj_out = gradient(j_out,φh)
 dj_vec_out = -assemble_vector(dj_out,V_φ)
 norm(dj_vec_out)
 
 Γ = EmbeddedBoundary(cutgeo)
 dΓ = Measure(Γ,3*order)
-dj_expected(q) = ∫(-q)dΓ
+dj_expected(q) = ∫(-fh*q/(norm ∘ (∇(φh))))dΓ
 dj_exp_vec = assemble_vector(dj_expected,V_φ)
 norm(dj_exp_vec)
 
 norm(dj_vec_in-dj_exp_vec)
 norm(dj_vec_out-dj_exp_vec)
 
-diff_Γ = DifferentiableTriangulation(Γ)
-dΓ2 = Measure(diff_Γ,3*order)
-jΓ(φ) = ∫(1)dΓ2
-djΓ = gradient(jΓ,φh)
-
 ############################################################################################
 
 include("../SubFacetSkeletons.jl")
@@ -87,35 +90,59 @@ dΛ = Measure(Λ,2*order)
 
 n_∂Ω = get_subfacet_normal_vector(Λ)
 n_k = get_ghost_normal_vector(Λ)
-#t_S = get_tangent_vector(Λ)
+# t_S = get_tangent_vector(Λ)
 n_S = get_normal_vector(Λ)
 m_k = get_conormal_vector(Λ)
 
-fh = interpolate(x->1,V_φ)
+fh = interpolate(x->cos(x[1]),V_φ)
 ∇ˢφ = Operation(abs)(n_S ⋅ ∇(φh).plus)
-dJ2(w) = ∫(n_S ⋅ (jump(fh*m_k) * mean(w) / ∇ˢφ))dΛ
+_n = get_normal_vector(Γ)
+dJ2(w) = ∫((-_n⋅∇(fh))*w/(norm ∘ (∇(φh))))dΓ + ∫(-n_S ⋅ (jump(fh*m_k) * mean(w) / ∇ˢφ))dΛ
 dj2 = assemble_vector(dJ2,V_φ)
 
 diff_Γ = DifferentiableTriangulation(Γ)
 dΓ_AD = Measure(diff_Γ,2*order)
 J2(φ) = ∫(fh)dΓ_AD
 dJ2_AD = gradient(J2,φh)
-dj2_AD = assemble_vector(dJ2_AD,V_φ)
+dj2_AD = assemble_vector(dJ2_AD,V_φ) # TODO: Why is this +ve but AD on volume is -ve???
+
+dj2 - dj2_AD
 
 ############################################################################################
 
+writevtk(
+  Λ,
+  "results/GammaSkel",
+  cellfields=[
+    "n_∂Ω.plus" => n_∂Ω.plus,"n_∂Ω.minus" => n_∂Ω.minus,
+    "n_k.plus" => n_k.plus,"n_k.minus" => n_k.minus,
+    # "t_S.plus" => t_S.plus,"t_S.minus" => t_S.minus,
+    "n_S" => n_S,
+    "m_k.plus" => m_k.plus,"m_k.minus" => m_k.minus,
+    "∇ˢφ"=>∇ˢφ,
+    "∇φh_Γs_plus"=>∇(φh).plus,"∇φh_Γs_minus"=>∇(φh).minus,
+    "jump(fh*m_k)"=>jump(fh*m_k)
+  ];
+  append=false
+)
+
+bgcell_to_inoutcut = compute_bgcell_to_inoutcut(model,geo)
 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)],
-  celldata=["inoutcut"=>bgcell_to_inoutcut]
+  cellfields=["φh"=>φh,"∇φh"=>∇(φh),"dj2_in"=>FEFunction(V_φ,dj2_AD),"dj_expected"=>FEFunction(V_φ,dj2)],
+  celldata=["inoutcut"=>bgcell_to_inoutcut];
+  append=false
 )
 
 writevtk(
-  Triangulation(cutgeo,PHYSICAL_IN),"results/trian_in"
+  Triangulation(cutgeo,PHYSICAL_IN),"results/trian_in";
+  append=false
 )
 writevtk(
-  Triangulation(cutgeo,PHYSICAL_OUT),"results/trian_out"
+  Triangulation(cutgeo,PHYSICAL_OUT),"results/trian_out";
+  append=false
 )
 writevtk(
-  Γ,"results/gamma"
+  Γ,"results/gamma";
+  append=false
 )
\ No newline at end of file
diff --git a/scripts/Embedded/SubFacetSkeletons.jl b/scripts/Embedded/SubFacetSkeletons.jl
index 075ed8f..9222549 100644
--- a/scripts/Embedded/SubFacetSkeletons.jl
+++ b/scripts/Embedded/SubFacetSkeletons.jl
@@ -61,7 +61,7 @@ function generate_ghost_trian(
   face_to_cell = get_faces(topo,Dc-1,Dc)
   cell_to_face = get_faces(topo,Dc,Dc-1)
 
-  n_bgfaces = num_faces(bgmodel,1)
+  n_bgfaces = num_faces(bgmodel,Dc-1)
   n_faces = num_cells(trian)
   ghost_faces = zeros(Int32,n_faces)
   p_lcell = ones(Int8,n_bgfaces)
@@ -131,7 +131,7 @@ function CellData.change_domain(
     # 1) CellField defined on the skeleton
     return change_domain(a,DomainStyle(a),tdomain)
   end
-  
+
   if is_change_possible(strian,ttrian.cell_skeleton)
     # 2) CellField defined on the bgmodel
     b = change_domain(a,ttrian.cell_skeleton,tdomain)
@@ -154,23 +154,44 @@ function CellData.change_domain(
 end
 
 # Normal vector to the cut facets , n_∂Ω
-function get_subfacet_normal_vector(trian::SubFacetSkeletonTriangulation) 
+function get_subfacet_normal_vector(trian::SubFacetSkeletonTriangulation{0})
   n_∂Ω = get_normal_vector(trian.face_trian)
   plus = change_domain(n_∂Ω.plus,trian,ReferenceDomain())
   minus = change_domain(-n_∂Ω.minus,trian,ReferenceDomain())
   return SkeletonPair(plus,minus)
 end
+function get_subfacet_normal_vector(trian::SubFacetSkeletonTriangulation{1})
+  n_∂Ω = get_normal_vector(trian.face_trian)
+  plus = change_domain(n_∂Ω.plus,trian,ReferenceDomain())
+  minus = change_domain(n_∂Ω.minus,trian,ReferenceDomain()) # This needs to be postive for mk to be correctly oriented in 3d
+  return SkeletonPair(plus,minus)
+end
 
 # Normal vector to the ghost facets, n_k
-function get_ghost_normal_vector(trian::SubFacetSkeletonTriangulation)
+# function get_ghost_normal_vector(trian::SubFacetSkeletonTriangulation)
+#   n = get_normal_vector(trian.ghost_skeleton)
+#   plus  = GenericCellField(CellData.get_data(n.plus),trian,ReferenceDomain())
+#   minus = GenericCellField(CellData.get_data(n.minus),trian,ReferenceDomain())
+#   return SkeletonPair(plus,minus)
+# end
+# The above fails as fields defined on SubFacetSkeletonTriangulation have 1 input in 3d points
+#   but fields defined on ghost skeleton have 2 inputs. Should change_domain fix this?
+
+function get_ghost_normal_vector(trian::SubFacetSkeletonTriangulation{0})
   n = get_normal_vector(trian.ghost_skeleton)
-  plus  = GenericCellField(CellData.get_data(n.plus),trian,ReferenceDomain())
-  minus = GenericCellField(CellData.get_data(n.minus),trian,ReferenceDomain())
+  plus  = CellField(evaluate(get_data(n.plus),Point(0,)),trian,ReferenceDomain())
+  minus = CellField(evaluate(get_data(n.minus),Point(0,)),trian,ReferenceDomain())
+  return SkeletonPair(plus,minus)
+end
+function get_ghost_normal_vector(trian::SubFacetSkeletonTriangulation{1})
+  n = get_normal_vector(trian.ghost_skeleton)
+  plus  = CellField(evaluate(get_data(n.plus),Point(0,0)),trian,ReferenceDomain())
+  minus = CellField(evaluate(get_data(n.minus),Point(0,0)),trian,ReferenceDomain())
   return SkeletonPair(plus,minus)
 end
 
 """
-# Returns a consistent tangent vector in the ReferenceDomain. Consistency 
+# Returns a consistent tangent vector in the ReferenceDomain. Consistency
 # is achieved by choosing it's direction as going from the node with the
 # smallest id to the one with the largest id.
 function get_tangent_vector(
@@ -214,14 +235,26 @@ function get_tangent_vector(
   return SkeletonPair(plus,minus)
 end
 
+_2d_cross(n) = VectorValue(-n[2],n[1]);
+function normalise(v) # <- Double check RE if this neccessary?
+  m = sqrt(inner(v,v))
+  if m < eps()
+    return zero(v)
+  else
+    return v/m
+  end
+end
+
 # Normal vector to the cut interface, n_S
 function CellData.get_normal_vector(trian::SubFacetSkeletonTriangulation{Di}) where {Di}
   if Di == 0
-    n_S = get_tangent_vector(trian.ghost_skeleton;ttrian=trian)
+    n_k = get_ghost_normal_vector(trian)
+    n_S = Operation(_2d_cross)(n_k) # nS = nk x tS and tS = ±e₃ in 2D
+    # n_S = get_tangent_vector(trian.ghost_skeleton;ttrian=trian) # Not oriented correctly in 2d. # TODO: understand why!?
   elseif Di == 1
     n_k = get_ghost_normal_vector(trian)
-    t_S = get_tangent_vector(trian.cell_skeleton;ttrian=trian)
-    n_S = Operation(cross)(n_k,t_S) # nk = tS x nS -> nS = nk x tS (eq 6.25)
+    t_S = get_tangent_vector(trian.face_skeleton;ttrian=trian)
+    n_S = Operation(normalise ∘ cross)(n_k,t_S) # nk = tS x nS -> nS = nk x tS (eq 6.25)
   else
     @notimplemented
   end
@@ -229,20 +262,40 @@ function CellData.get_normal_vector(trian::SubFacetSkeletonTriangulation{Di}) wh
   # computing t_S as t_S = n_S x n_k
   return n_S.plus
 end
+# function alt_get_normal_vector(trian::SubFacetSkeletonTriangulation{Di}) where {Di}
+#   if Di == 0
+#     n_k = alt_get_ghost_normal_vector(trian)
+#     n_S = Operation(_2d_cross)(n_k) # nS = nk x tS and tS = ±e₃ in 2D
+#   elseif Di == 1
+#     n_k = alt_get_ghost_normal_vector(trian)
+#     t_S = alt_get_tangent_vector(trian)
+#     n_S = Operation(cross)(n_k,t_S) # nk = tS x nS -> nS = nk x tS (eq 6.25)
+#   else
+#     @notimplemented
+#   end
+#   # We pick the plus side. Consistency is given by later
+#   # computing t_S as t_S = n_S x n_k
+#   return n_S.plus
+# end
 
 # Tangent vector to the cut interface, t_S = n_S x n_k
 function get_tangent_vector(trian::SubFacetSkeletonTriangulation{Di}) where {Di}
   @notimplementedif Di != 1
   n_S = get_normal_vector(trian)
   n_k = get_ghost_normal_vector(trian)
-  return Operation(cross)(n_S,n_k)
+  return Operation(normalise ∘ cross)(n_S,n_k)
 end
+# function alt_get_tangent_vector(trian::SubFacetSkeletonTriangulation{1})
+#   n_∂Ω = get_subfacet_normal_vector(trian)
+#   n_k = get_ghost_normal_vector(trian)
+#   return Operation(normalise ∘ cross)(n_k,n_∂Ω) # t_S = n_k × n_∂Ω # <- need to show this if we wanted to actually use it!
+# end
 
 # Conormal vectors, m_k = t_S x n_∂Ω
 function get_conormal_vector(trian::SubFacetSkeletonTriangulation{Di}) where {Di}
   op = Operation(GridapEmbedded.LevelSetCutters._normal_vector)
   n_∂Ω = get_subfacet_normal_vector(trian)
-  if Di == 0 # 2D 
+  if Di == 0 # 2D
     m_k = op(n_∂Ω)
   elseif Di == 1 # 3D
     t_S = get_tangent_vector(trian)
@@ -252,6 +305,19 @@ function get_conormal_vector(trian::SubFacetSkeletonTriangulation{Di}) where {Di
   end
   return m_k
 end
+# function alt_get_conormal_vector(trian::SubFacetSkeletonTriangulation{Di}) where {Di}
+#   op = Operation(GridapEmbedded.LevelSetCutters._normal_vector)
+#   n_∂Ω = get_subfacet_normal_vector(trian)
+#   if Di == 0 # 2D
+#     m_k = op(n_∂Ω)
+#   elseif Di == 1 # 3D
+#     t_S = alt_get_tangent_vector(trian)
+#     m_k = op(t_S,n_∂Ω) # m_k = t_S x n_∂Ω (eq 6.26)
+#   else
+#     @notimplemented
+#   end
+#   return m_k
+# end
 
 # This will go to Gridap