diff --git a/Project.toml b/Project.toml
index e048540..f277a42 100755
--- a/Project.toml
+++ b/Project.toml
@@ -24,15 +24,16 @@ SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
 
 [compat]
 BenchmarkTools = "1"
-BlockArrays = "0.16"
+BlockArrays = "0.16,1"
 ChainRulesCore = "1"
 CircularArrays = "1"
 FillArrays = "0.8.4,1"
 Gridap = "0.18"
 GridapDistributed = "0.4"
+GridapEmbedded = "0.9.5"
 GridapPETSc = "0.5"
-GridapSolvers = "0.3"
-JLD2 = "0.4"
+GridapSolvers = "0.4"
+JLD2 = "0.4,0.5"
 MPI = "0.19, 0.20"
 PartitionedArrays = "0.3"
 SparseMatricesCSR = "0.6.6"
diff --git a/src/Embedded/DifferentiableTriangulations.jl b/src/Embedded/DifferentiableTriangulations.jl
new file mode 100644
index 0000000..8780399
--- /dev/null
+++ b/src/Embedded/DifferentiableTriangulations.jl
@@ -0,0 +1,373 @@
+
+# DifferentiableTriangulation
+
+"""
+  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 
+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 the background mesh.
+"""
+mutable struct DifferentiableTriangulation{Dc,Dp,A} <: Triangulation{Dc,Dp}
+  trian :: A
+  cell_values
+  caches
+  function DifferentiableTriangulation(
+    trian :: Triangulation{Dc,Dp},cell_values,caches
+  ) where {Dc,Dp}
+    A = typeof(trian)
+    new{Dc,Dp,A}(trian,cell_values,caches)
+  end
+end
+
+function DifferentiableTriangulation(
+  trian::Union{<:SubCellTriangulation,<:SubFacetTriangulation}
+)
+  caches = precompute_autodiff_caches(trian)
+  return DifferentiableTriangulation(trian,nothing,caches)
+end
+
+(t::DifferentiableTriangulation)(φh) = update_trian!(t,φh)
+
+function update_trian!(trian::DifferentiableTriangulation,φh)
+  trian.cell_values = extract_dualized_cell_values(trian.trian,φh)
+  return trian
+end
+
+function FESpaces._change_argument(
+  op,f,trian::DifferentiableTriangulation,uh::SingleFieldFEFunction
+)
+  U = get_fe_space(uh)
+  function g(cell_u)
+    cf = CellField(U,cell_u)
+    update_trian!(trian,cf)
+    cell_grad = f(cf)
+    get_contribution(cell_grad,trian)
+  end
+  g
+end
+
+function FESpaces._compute_cell_ids(uh,ttrian::DifferentiableTriangulation)
+  FESpaces._compute_cell_ids(uh,ttrian.trian)
+end
+
+function Geometry.get_background_model(t::DifferentiableTriangulation)
+  get_background_model(t.trian)
+end
+
+function Geometry.get_grid(t::DifferentiableTriangulation)
+  get_grid(t.trian)
+end
+
+function Geometry.get_cell_reffe(t::DifferentiableTriangulation)
+  get_cell_reffe(t.trian)
+end
+
+function CellData.get_cell_points(ttrian::DifferentiableTriangulation)
+  if isnothing(ttrian.cell_values)
+    return get_cell_points(ttrian.trian)
+  end
+  c = ttrian.caches
+  cell_values = ttrian.cell_values
+  cell_to_rcoords = lazy_map(
+    DualizeCoordsMap(),c.face_to_rcoords,c.face_to_bgrcoords,
+    cell_values,c.face_to_edges,c.face_to_edge_lists
+  )
+  cell_to_coords = lazy_map(
+    DualizeCoordsMap(),c.face_to_coords,c.face_to_bgcoords,
+    cell_values,c.face_to_edges,c.face_to_edge_lists
+  )
+  return CellPoint(cell_to_rcoords, cell_to_coords, ttrian, ReferenceDomain())
+end
+
+function Geometry.get_cell_map(ttrian::DifferentiableTriangulation)
+  if isnothing(ttrian.cell_values)
+    return get_cell_map(ttrian.trian)
+  end
+  c = ttrian.caches
+  cell_values = ttrian.cell_values
+  cell_to_coords = lazy_map(
+    DualizeCoordsMap(),c.face_to_coords,c.face_to_bgcoords,
+    cell_values,c.face_to_edges,c.face_to_edge_lists
+  )
+  cell_reffe = get_cell_reffe(ttrian)
+  cell_map = compute_cell_maps(cell_to_coords,cell_reffe)
+  return cell_map
+end
+
+function Geometry.get_glue(ttrian::DifferentiableTriangulation,val::Val{D}) where {D}
+  glue = get_glue(ttrian.trian,val)
+  if isnothing(glue) || isnothing(ttrian.cell_values)
+    return glue
+  end
+
+  # New reference maps 
+  c = ttrian.caches
+  cell_values = ttrian.cell_values
+  cell_to_rcoords = lazy_map(
+    DualizeCoordsMap(),c.face_to_rcoords,c.face_to_bgrcoords,
+    cell_values,c.face_to_edges,c.face_to_edge_lists
+  )
+  cell_reffe = get_cell_reffe(ttrian)
+  ref_cell_map = compute_cell_maps(cell_to_rcoords,cell_reffe)
+
+  return FaceToFaceGlue(
+    glue.tface_to_mface,
+    ref_cell_map,
+    glue.mface_to_tface,
+  )
+end
+
+function FESpaces.get_cell_fe_data(fun,f,ttrian::DifferentiableTriangulation)
+  FESpaces.get_cell_fe_data(fun,f,ttrian.trian)
+end
+
+function compute_cell_maps(cell_coords,cell_reffes)
+  cell_shapefuns = lazy_map(get_shapefuns,cell_reffes)
+  default_cell_map = lazy_map(linear_combination,cell_coords,cell_shapefuns)
+  default_cell_grad = lazy_map(∇,default_cell_map)
+  cell_poly = lazy_map(get_polytope,cell_reffes)
+  cell_q0 = lazy_map(p->zero(first(get_vertex_coordinates(p))),cell_poly)
+  origins = lazy_map(evaluate,default_cell_map,cell_q0)
+  gradients = lazy_map(evaluate,default_cell_grad,cell_q0)
+  cell_map = lazy_map(Gridap.Fields.affine_map,gradients,origins)
+  return cell_map
+end
+
+# DualizeCoordsMap
+
+struct DualizeCoordsMap <: Map end
+
+function Arrays.return_cache(
+  k::DualizeCoordsMap,
+  coords::Vector{<:Point{Dp,Tp}},
+  bg_coords::Vector{<:Point{Dp,Tp}},
+  values::Vector{Tv},
+  edges::Vector{Int8},
+  edge_list::Vector{Vector{Int8}}
+) where {Dp,Tp,Tv}
+  T = Point{Dp,Tv}
+  return CachedArray(zeros(T, length(coords)))
+end
+
+function Arrays.evaluate!(
+  cache,
+  k::DualizeCoordsMap,
+  coords::Vector{<:Point{Dp,Tp}},
+  bg_coords::Vector{<:Point{Dp,Tp}},
+  values::Vector{Tv},
+  edges::Vector{Int8},
+  edge_list::Vector{Vector{Int8}}
+) where {Dp,Tp,Tv}
+  setsize!(cache,(length(coords),))
+  new_coords = cache.array
+  for (i,e) in enumerate(edges)
+    if e == -1
+      new_coords[i] = coords[i]
+    else
+      n1, n2 = edge_list[e]
+      q1, q2 = bg_coords[n1], bg_coords[n2]
+      v1, v2 = abs(values[n1]), abs(values[n2])
+      λ = v1/(v1+v2)
+      new_coords[i] = q1 + λ*(q2-q1)
+    end
+  end
+  return new_coords
+end
+
+"""
+    precompute_cut_edge_ids(rcoords,bg_rcoords,edge_list)
+
+Given 
+  - `rcoords`: the node ref coordinates of the cut subcell/subfacet,
+  - `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
+  - `-1` if the node is also a node of the background cell,
+  - the id of the edge containing the node otherwise.
+"""
+function precompute_cut_edge_ids(
+  rcoords::Vector{<:Point{Dp,Tp}},
+  bg_rcoords::Vector{<:Point{Dp,Tp}},
+  edge_list::Vector{<:Vector{<:Integer}}
+) where {Dp,Tp}
+  tol = 10*eps(Tp)
+  edges = Vector{Int8}(undef,length(rcoords))
+  for (i,p) in enumerate(rcoords)
+    if any(q -> norm(q-p) < tol, bg_rcoords)
+      edges[i] = Int8(-1)
+    else
+      e = findfirst(edge -> belongs_to_edge(p,edge,bg_rcoords), edge_list)
+      edges[i] = Int8(e)
+    end
+  end
+  return edges
+end
+
+function get_edge_list(poly::Polytope)
+  ltcell_to_lpoints, simplex = simplexify(poly)
+  simplex_edges = get_faces(simplex,1,0)
+  ltcell_to_edges = map(pts -> map(e -> pts[e], simplex_edges), ltcell_to_lpoints)
+  return collect(Vector{Int8},unique(sort,vcat(ltcell_to_edges...)))
+end
+
+function belongs_to_edge(
+  p::Point{D,T},edge::Vector{<:Integer},bgpts::Vector{Point{D,T}}
+) where {D,T}
+  tol = 10*eps(T)
+  p1, p2 = bgpts[edge]
+  return norm(cross(p-p1,p2-p1)) < tol
+end
+
+function precompute_autodiff_caches(
+  trian::SubCellTriangulation
+)
+  bgmodel = get_background_model(trian)
+  subcells = trian.subcells
+
+  precompute_autodiff_caches(
+    bgmodel,
+    subcells.cell_to_bgcell,
+    subcells.cell_to_points,
+    subcells.point_to_rcoords,
+    subcells.point_to_coords,
+  )
+end
+
+function precompute_autodiff_caches(
+  trian::SubFacetTriangulation
+)
+  bgmodel = get_background_model(trian)
+  subfacets = trian.subfacets
+
+  precompute_autodiff_caches(
+    bgmodel,
+    subfacets.facet_to_bgcell,
+    subfacets.facet_to_points,
+    subfacets.point_to_rcoords,
+    subfacets.point_to_coords,
+  )
+end
+
+function precompute_autodiff_caches(
+  bgmodel,
+  face_to_bgcell,
+  face_to_points,
+  point_to_rcoords,
+  point_to_coords,
+)
+  bg_ctypes = get_cell_type(bgmodel)
+  bgcell_to_polys = expand_cell_data(get_polytopes(bgmodel),bg_ctypes)
+  bgcell_to_coords = get_cell_coordinates(bgmodel)
+  bgcell_to_rcoords = lazy_map(get_vertex_coordinates,bgcell_to_polys)
+
+  face_to_bgcoords = lazy_map(Reindex(bgcell_to_coords),face_to_bgcell)
+  face_to_bgrcoords = lazy_map(Reindex(bgcell_to_rcoords),face_to_bgcell)
+  face_to_rcoords = lazy_map(Broadcasting(Reindex(point_to_rcoords)),face_to_points)
+  face_to_coords = lazy_map(Broadcasting(Reindex(point_to_coords)),face_to_points)
+
+  bgcell_to_edge_lists = lazy_map(get_edge_list,bgcell_to_polys)
+  face_to_edge_lists = lazy_map(Reindex(bgcell_to_edge_lists),face_to_bgcell)
+  face_to_edges = collect(lazy_map(precompute_cut_edge_ids,face_to_rcoords,face_to_bgrcoords,face_to_edge_lists))
+
+  cache = (;
+    face_to_rcoords,
+    face_to_coords,
+    face_to_bgrcoords,
+    face_to_bgcoords,
+    face_to_edges,
+    face_to_edge_lists
+  )
+  return cache
+end
+
+function extract_dualized_cell_values(
+  trian::SubCellTriangulation,
+  φh::CellField,
+)
+  @assert isa(DomainStyle(φh),ReferenceDomain)
+  bgmodel = get_background_model(trian)
+  bgcell_to_values = extract_dualized_cell_values(bgmodel,φh)
+
+  subcells = trian.subcells
+  cell_to_bgcell = subcells.cell_to_bgcell
+  cell_to_values = lazy_map(Reindex(bgcell_to_values),cell_to_bgcell)
+  return cell_to_values
+end
+
+function extract_dualized_cell_values(
+  trian::SubFacetTriangulation,
+  φh::CellField,
+)
+  @assert isa(DomainStyle(φh),ReferenceDomain)
+  bgmodel = get_background_model(trian)
+  bgcell_to_values = extract_dualized_cell_values(bgmodel,φh)
+
+  subfacets = trian.subfacets
+  facet_to_bgcell = subfacets.facet_to_bgcell
+  facet_to_values = lazy_map(Reindex(bgcell_to_values),facet_to_bgcell)
+  return facet_to_values
+end
+
+function extract_dualized_cell_values(
+  bgmodel::DiscreteModel,
+  φh::CellField,
+)
+  @assert isa(DomainStyle(φh),ReferenceDomain)
+  bg_ctypes = get_cell_type(bgmodel)
+  bgcell_to_polys = expand_cell_data(get_polytopes(bgmodel),bg_ctypes)
+  bgcell_to_rcoords = lazy_map(get_vertex_coordinates,bgcell_to_polys)
+  bgcell_to_fields = CellData.get_data(φh)
+  bgcell_to_values = lazy_map(evaluate,bgcell_to_fields,bgcell_to_rcoords)
+  return bgcell_to_values
+end
+
+# AppendedTriangulation
+#
+# When cutting an embedded domain, we will usually end up with an AppendedTriangulation
+# 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 
+# following implementation does:
+
+DifferentiableTriangulation(trian::Triangulation) = trian
+
+function DifferentiableTriangulation(
+  trian::AppendedTriangulation
+)
+  a = DifferentiableTriangulation(trian.a)
+  b = DifferentiableTriangulation(trian.b)
+  return AppendedTriangulation(a,b)
+end
+
+update_trian!(trian::Triangulation,φh) = trian
+
+function update_trian!(trian::AppendedTriangulation,φh)
+  update_trian!(trian.a,φh)
+  update_trian!(trian.b,φh)
+  return trian
+end
+
+function FESpaces._change_argument(
+  op,f,trian::AppendedTriangulation,uh::SingleFieldFEFunction
+)
+  U = get_fe_space(uh)
+  function g(cell_u)
+    cf = CellField(U,cell_u)
+    update_trian!(trian,cf)
+    cell_grad = f(cf)
+    get_contribution(cell_grad,trian)
+  end
+  g
+end
+
+function FESpaces._compute_cell_ids(uh,ttrian::AppendedTriangulation)
+  ids_a = FESpaces._compute_cell_ids(uh,ttrian.a)
+  ids_b = FESpaces._compute_cell_ids(uh,ttrian.b)
+  lazy_append(ids_a,ids_b)
+end
diff --git a/src/Embedded/Embedded.jl b/src/Embedded/Embedded.jl
new file mode 100644
index 0000000..4edefd7
--- /dev/null
+++ b/src/Embedded/Embedded.jl
@@ -0,0 +1,3 @@
+
+include("DifferentiableTriangulations.jl")
+include("SubFacetBoundaryTriangulations.jl")
diff --git a/src/Embedded/SubFacetBoundaryTriangulations.jl b/src/Embedded/SubFacetBoundaryTriangulations.jl
new file mode 100644
index 0000000..fb2dbb0
--- /dev/null
+++ b/src/Embedded/SubFacetBoundaryTriangulations.jl
@@ -0,0 +1,508 @@
+
+using Gridap.Geometry: CompositeTriangulation
+
+Base.round(a::VectorValue{D,T};kwargs...) where {D,T} = VectorValue{D,T}(round.(a.data;kwargs...))
+
+function consistent_facet_to_points(
+  facet_to_points::Table, point_to_coords::Vector
+)
+  f(pt::Number) = round(pt;sigdigits=12)
+  f(id::Integer) = f(point_to_coords[id])
+
+  # Create a list of the unique points composing the facets
+  npts = length(point_to_coords)
+  nfaces = length(facet_to_points)
+  touched = zeros(Bool,npts)
+  for face in 1:nfaces
+    pts = view(facet_to_points,face)
+    touched[pts] .= true
+  end
+  touched_ids = findall(touched)
+  unique_ids = unique(f,touched_ids)
+
+  # Create a mapping from the old point ids to the new ones
+  touched_to_uid = collect(Int32,indexin(f.(touched_ids),f.(unique_ids)))
+  point_to_uid = extend(touched_to_uid,PosNegPartition(touched_ids,npts))
+
+  facet_to_uids = Table(
+    collect(Int32,lazy_map(Reindex(point_to_uid),facet_to_points.data)),
+    facet_to_points.ptrs
+  )
+  return facet_to_uids, unique_ids
+end
+
+function Geometry.get_active_model(
+  trian::SubFacetTriangulation{Dc,Dp}
+) where {Dc,Dp}
+  subgrid = trian.subgrid
+  subfacets = trian.subfacets
+  facet_to_uids, uid_to_point = consistent_facet_to_points(
+    subfacets.facet_to_points,subfacets.point_to_coords
+  )
+  topo = UnstructuredGridTopology(
+    subgrid,facet_to_uids,uid_to_point
+  )
+  return UnstructuredDiscreteModel(subgrid,topo,FaceLabeling(topo))
+end
+
+function generate_ghost_trian(
+  trian::CompositeTriangulation,
+  bgmodel::UnstructuredDiscreteModel{Dc}
+) where {Dc}
+  cell_glue = get_glue(trian,Val(Dc))
+  return generate_ghost_trian(trian,bgmodel,cell_glue)
+end
+
+function generate_ghost_trian(
+  trian::CompositeTriangulation,
+  bgmodel::UnstructuredDiscreteModel{Dc},
+  cell_glue::SkeletonPair{<:FaceToFaceGlue}
+) where {Dc}
+  topo = get_grid_topology(bgmodel)
+  face_to_cell = get_faces(topo,Dc-1,Dc)
+  cell_to_face = get_faces(topo,Dc,Dc-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)
+  m_lcell = ones(Int8,n_bgfaces)
+  for (i,(p_cell, m_cell)) in enumerate(zip(cell_glue.plus.tface_to_mface,cell_glue.minus.tface_to_mface))
+    inter = intersect(cell_to_face[p_cell],cell_to_face[m_cell])
+    @assert length(inter) == 1
+    face = first(inter)
+    ghost_faces[i] = face
+
+    nbors = face_to_cell[ghost_faces[i]]
+    p_lcell[face] = findfirst(==(p_cell),nbors)
+    m_lcell[face] = findfirst(==(m_cell),nbors)
+  end
+
+  plus = BoundaryTriangulation(bgmodel,ghost_faces,p_lcell)
+  minus = BoundaryTriangulation(bgmodel,ghost_faces,m_lcell)
+  return SkeletonTriangulation(plus,minus)
+end
+
+function generate_ghost_trian(
+  trian::CompositeTriangulation,
+  bgmodel::UnstructuredDiscreteModel{Dc},
+  cell_glue::FaceToFaceGlue
+) where {Dc}
+  topo = get_grid_topology(bgmodel)
+  face_to_cell = get_faces(topo,Dc-1,Dc)
+  cell_to_face = get_faces(topo,Dc,Dc-1)
+  is_boundary(f) = isone(length(view(face_to_cell,f)))
+
+  n_faces = num_cells(trian)
+  ghost_faces = zeros(Int32,n_faces)
+  for (i,cell) in enumerate(cell_glue.tface_to_mface)
+    faces = filter(is_boundary,view(cell_to_face,cell))
+    @assert length(faces) == 1 # TODO: This will break if we are in a corner
+    face = first(faces)
+    ghost_faces[i] = face
+  end
+
+  # NOTE: lcell is always 1 for boundary facets
+  return BoundaryTriangulation(bgmodel,ghost_faces)
+end
+
+"""
+    get_ghost_mask(
+      face_trian::SubFacetTriangulation{Df,Dc},
+      face_model::UnstructuredDiscreteModel{Df,Dc}
+    ) where {Df,Dc}
+
+    get_ghost_mask(face_trian::SubFacetTriangulation) = get_ghost_mask(face_trian,get_active_model(face_trian))
+
+  Returns a mask for ghost faces. We define ghost faces as the interfaces between two
+  different cut facets that are located in different background cells.
+
+  The second condition is important: In 3D, some cuts subcells may not be simplices.
+  In this case, we simplexify the subcell. This creates extra cut interfaces that are
+  interior to a background cell. These are not considered ghost faces.
+
+  - In 2D: Dc = 2, Df = 1 -> Ghost faces have dimension 0 (i.e interface points)
+  - In 3D: Dc = 3, Df = 2 -> Ghost faces have dimension 1 (i.e interface edges)
+"""
+function get_ghost_mask(
+  face_trian::SubFacetTriangulation{Df,Dc},
+  face_model::UnstructuredDiscreteModel{Df,Dc}
+) where {Df,Dc}
+  topo = get_grid_topology(face_model)
+  face_to_facets = get_faces(topo,Df-1,Df)
+
+  subfacets = face_trian.subfacets
+  facet_to_bgcell = subfacets.facet_to_bgcell
+
+  n_faces = num_faces(topo,Df-1)
+  face_is_ghost = zeros(Bool,n_faces)
+  for face in 1:n_faces
+    facets = view(face_to_facets,face)
+    is_boundary = isone(length(facets))
+    if !is_boundary
+      @assert length(facets) == 2
+      bgcells = view(facet_to_bgcell,facets)
+      is_ghost = (bgcells[1] != bgcells[2])
+      face_is_ghost[face] = is_ghost
+    end
+  end
+
+  return face_is_ghost
+end
+
+function get_ghost_mask(
+  face_trian::SubFacetTriangulation
+)
+  face_model = get_active_model(face_trian)
+  return get_ghost_mask(face_trian,face_model)
+end
+
+struct SubFacetSkeletonTriangulation{Di,Df,Dp} <: Triangulation{Di,Dp}
+  cell_skeleton::CompositeTriangulation{Di,Dp} # Interface -> BG Cell 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}
+end
+
+function Geometry.SkeletonTriangulation(face_trian::SubFacetTriangulation)
+  bgmodel = get_background_model(face_trian)
+  face_model = get_active_model(face_trian)
+
+  ghost_mask = get_ghost_mask(face_trian,face_model)
+  face_skeleton = SkeletonTriangulation(face_model,ghost_mask)
+  cell_skeleton = CompositeTriangulation(face_trian,face_skeleton)
+  ghost_skeleton = generate_ghost_trian(cell_skeleton,bgmodel)
+  return SubFacetSkeletonTriangulation(cell_skeleton,face_skeleton,ghost_skeleton,face_trian,face_model)
+end
+
+function Geometry.get_background_model(t::SubFacetSkeletonTriangulation)
+  get_background_model(t.cell_skeleton)
+end
+
+function Geometry.get_active_model(t::SubFacetSkeletonTriangulation)
+  get_active_model(t.cell_skeleton)
+end
+
+function Geometry.get_grid(t::SubFacetSkeletonTriangulation)
+  get_grid(t.cell_skeleton)
+end
+
+# Domain changes
+
+function Geometry.get_glue(ttrian::SubFacetSkeletonTriangulation{Di,Df,Dp},::Val{D}) where {D,Di,Df,Dp}
+  get_glue(ttrian.cell_skeleton,Val(D))
+end
+
+function Geometry.is_change_possible(
+  strian::SubFacetTriangulation,ttrian::SubFacetSkeletonTriangulation
+)
+  return strian === ttrian.face_trian
+end
+
+function CellData.change_domain(
+  a::CellField,ttrian::SubFacetSkeletonTriangulation,tdomain::DomainStyle
+)
+  strian = get_triangulation(a)
+  if strian === ttrian
+    # 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)
+  elseif strian === ttrian.face_trian
+    # 3) CellField defined on the cut facets
+    itrian = Triangulation(ttrian.face_model)
+    _a = CellData.similar_cell_field(a,CellData.get_data(a),itrian,DomainStyle(a))
+    b = change_domain(_a,ttrian.face_skeleton,tdomain)
+  else
+    @notimplemented
+  end
+  return CellData.similar_cell_field(b,CellData.get_data(b),ttrian,DomainStyle(b))
+end
+
+function CellData.change_domain(
+  f::CellData.OperationCellField,ttrian::SubFacetSkeletonTriangulation,tdomain::DomainStyle
+)
+  args = map(i->change_domain(i,ttrian,tdomain),f.args)
+  CellData.OperationCellField(f.op,args...)
+end
+
+# Normal vector to the cut facets , n_∂Ω
+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)
+#   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{Di}) where Di
+  n = get_normal_vector(trian.ghost_skeleton)
+  z = zero(VectorValue{Di+1,Float64})
+  plus  = CellField(evaluate(get_data(n.plus),z),trian,ReferenceDomain())
+  minus = CellField(evaluate(get_data(n.minus),z),trian,ReferenceDomain())
+  return SkeletonPair(plus,minus)
+end
+
+# Returns the tangent vector in the PhysicalDomain
+function get_tangent_vector(
+  trian::BoundaryTriangulation{1};
+  ttrian = trian
+)
+  function t(c)
+    @assert length(c) == 2
+    t = c[2] - c[1]
+    return t/norm(t)
+  end
+  data = lazy_map(constant_field,lazy_map(t,get_cell_coordinates(trian)))
+  return GenericCellField(data,ttrian,ReferenceDomain())
+end
+
+function get_tangent_vector(
+  trian::SkeletonTriangulation{1};
+  ttrian = trian
+)
+  plus = get_tangent_vector(trian.plus;ttrian=ttrian)
+  minus = get_tangent_vector(trian.minus;ttrian=ttrian)
+  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_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.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
+  # 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(normalise ∘ cross)(n_S,n_k)
+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
+    m_k = op(n_∂Ω)
+  elseif Di == 1 # 3D
+    t_S = 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
+
+############################################################################################
+# SubFacetBoundaryTriangulation
+
+# TODO: In the future, I imagine unifying both structures, and simply having the Boundary one
+# that can be filtered depending on which interfaces we are aiming at, and the
+# skeleton just being two of those (kinda like gridap does it). For now, let's
+# keep them separate until we figure out all the tangent/normal/conormal stuff.
+
+struct SubFacetBoundaryTriangulation{Di,Df,Dp} <: Triangulation{Di,Dp}
+  cell_boundary::CompositeTriangulation{Di,Dp} # Interface -> BG Cell
+  face_boundary::BoundaryTriangulation{Di,Dp}  # Interface -> Cut Facet
+  ghost_boundary::BoundaryTriangulation{Df,Dp} # Ghost Facet -> BG Cell
+  face_trian::SubFacetTriangulation{Df,Dp}     # Cut Facet -> BG Cell
+  face_model::UnstructuredDiscreteModel{Df,Dp}
+end
+
+function Geometry.BoundaryTriangulation(face_trian::SubFacetTriangulation)
+  bgmodel = get_background_model(face_trian)
+  face_model = get_active_model(face_trian)
+
+  face_boundary = BoundaryTriangulation(face_model)
+  cell_boundary = CompositeTriangulation(face_trian,face_boundary)
+  ghost_boundary = generate_ghost_trian(cell_boundary,bgmodel)
+  return SubFacetBoundaryTriangulation(cell_boundary,face_boundary,ghost_boundary,face_trian,face_model)
+end
+
+function Geometry.get_background_model(t::SubFacetBoundaryTriangulation)
+  get_background_model(t.cell_boundary)
+end
+
+function Geometry.get_active_model(t::SubFacetBoundaryTriangulation)
+  get_active_model(t.cell_boundary)
+end
+
+function Geometry.get_grid(t::SubFacetBoundaryTriangulation)
+  get_grid(t.cell_boundary)
+end
+
+# Domain changes
+
+function Geometry.get_glue(ttrian::SubFacetBoundaryTriangulation{Di,Df,Dp},::Val{D}) where {D,Di,Df,Dp}
+  get_glue(ttrian.cell_boundary,Val(D))
+end
+
+function Geometry.is_change_possible(
+  strian::SubFacetTriangulation,ttrian::SubFacetBoundaryTriangulation
+)
+  return strian === ttrian.face_trian
+end
+
+function CellData.change_domain(
+  a::CellField,ttrian::SubFacetBoundaryTriangulation,tdomain::DomainStyle
+)
+  strian = get_triangulation(a)
+  if strian === ttrian
+    # 1) CellField defined on the skeleton
+    return change_domain(a,DomainStyle(a),tdomain)
+  end
+
+  if is_change_possible(strian,ttrian.cell_boundary)
+    # 2) CellField defined on the bgmodel
+    b = change_domain(a,ttrian.cell_boundary,tdomain)
+  elseif strian === ttrian.face_trian
+    # 3) CellField defined on the cut facets
+    itrian = Triangulation(ttrian.face_model)
+    _a = CellData.similar_cell_field(a,CellData.get_data(a),itrian,DomainStyle(a))
+    b = change_domain(_a,ttrian.face_boundary,tdomain)
+  else
+    @notimplemented
+  end
+  return CellData.similar_cell_field(b,CellData.get_data(b),ttrian,DomainStyle(b))
+end
+
+function CellData.change_domain(
+  f::CellData.OperationCellField,ttrian::SubFacetBoundaryTriangulation,tdomain::DomainStyle
+)
+  args = map(i->change_domain(i,ttrian,tdomain),f.args)
+  CellData.OperationCellField(f.op,args...)
+end
+
+# Normal vector to the cut facets , n_∂Ω
+function get_subfacet_normal_vector(trian::SubFacetBoundaryTriangulation)
+  n_∂Ω = get_normal_vector(trian.face_trian)
+  return change_domain(n_∂Ω,trian,ReferenceDomain())
+end
+
+function get_ghost_normal_vector(trian::SubFacetBoundaryTriangulation{Di}) where Di
+  n = get_normal_vector(trian.ghost_boundary)
+  z = zero(VectorValue{Di+1,Float64})
+  return CellField(evaluate(get_data(n.plus),z),trian,ReferenceDomain())
+end
+
+# Normal vector to the cut interface, n_S
+function CellData.get_normal_vector(trian::SubFacetBoundaryTriangulation{Di}) where {Di}
+  if Di == 0
+    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.face_boundary;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
+  # We pick the plus side. Consistency is given by later
+  # computing t_S as t_S = n_S x n_k
+  return n_S
+end
+
+# 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
+  n_S = get_normal_vector(trian)
+  n_k = get_ghost_normal_vector(trian)
+  return Operation(normalise ∘ cross)(n_S,n_k)
+end
+
+# Conormal vectors, m_k = t_S x n_∂Ω
+function get_conormal_vector(trian::SubFacetBoundaryTriangulation{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 = 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
+
+function Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:CellField})
+  plus = k(a.plus)
+  minus = k(a.minus)
+  SkeletonPair(plus,minus)
+end
+
+function Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:CellField},b::SkeletonPair{<:CellField})
+  plus = k(a.plus,b.plus)
+  minus = k(a.minus,b.minus)
+  SkeletonPair(plus,minus)
+end
+
+import Gridap.TensorValues: inner, outer
+import LinearAlgebra: dot
+import Base: abs, *, +, -, /
+
+for op in (:/,)
+  @eval begin
+    ($op)(a::CellField,b::SkeletonPair{<:CellField}) = Operation($op)(a,b)
+    ($op)(a::SkeletonPair{<:CellField},b::CellField) = Operation($op)(a,b)
+  end
+end
+
+for op in (:outer,:*,:dot,:/)
+  @eval begin
+    ($op)(a::SkeletonPair{<:CellField},b::SkeletonPair{<:CellField}) = Operation($op)(a,b)
+  end
+end
+
+function CellData.change_domain(a::SkeletonPair, ::ReferenceDomain, ::PhysicalDomain)
+  plus = change_domain(a.plus,ReferenceDomain(),PhysicalDomain())
+  minus = change_domain(a.minus,ReferenceDomain(),PhysicalDomain())
+  return SkeletonPair(plus,minus)
+end
diff --git a/src/GridapTopOpt.jl b/src/GridapTopOpt.jl
index d0e654b..c19784f 100644
--- a/src/GridapTopOpt.jl
+++ b/src/GridapTopOpt.jl
@@ -4,17 +4,17 @@ using GridapPETSc, GridapPETSc.PETSC
 using GridapPETSc: PetscScalar, PetscInt, PETSC,  @check_error_code
 
 using MPI
-using BlockArrays
-using CircularArrays
+using BlockArrays, CircularArrays, FillArrays
 using LinearAlgebra
 using ChainRulesCore
 using DelimitedFiles, Printf
 
 using Gridap
 using Gridap.Helpers, Gridap.Algebra, Gridap.TensorValues
-using Gridap.Geometry, Gridap.CellData, Gridap.Fields
+using Gridap.Geometry, Gridap.CellData, Gridap.Fields, Gridap.Arrays
 using Gridap.ReferenceFEs, Gridap.FESpaces,  Gridap.MultiField
-using Gridap.Geometry: get_faces
+
+using Gridap.Geometry: get_faces, num_nodes
 using Gridap.FESpaces: get_assembly_strategy
 using Gridap: writevtk
 
@@ -31,6 +31,10 @@ using GridapSolvers
 using GridapSolvers.LinearSolvers, GridapSolvers.NonlinearSolvers, GridapSolvers.BlockSolvers
 using GridapSolvers.SolverInterfaces: SolverVerboseLevel, SOLVER_VERBOSE_NONE, SOLVER_VERBOSE_LOW, SOLVER_VERBOSE_HIGH
 
+using GridapEmbedded
+using GridapEmbedded.LevelSetCutters, GridapEmbedded.Interfaces
+using GridapEmbedded.Interfaces: SubFacetData, SubCellTriangulation, SubFacetTriangulation
+
 using JLD2: save_object, load_object, jldsave
 
 import Base: +
@@ -46,6 +50,11 @@ export get_state
 export evaluate_functionals!
 export evaluate_derivatives!
 
+include("Embedded/Embedded.jl")
+export DifferentiableTriangulation
+export SubFacetBoundaryTriangulation
+export SubFacetSkeletonTriangulation
+
 include("Utilities.jl")
 export SmoothErsatzMaterialInterpolation
 export update_labels!