diff --git a/scripts/Embedded/2d_thermal.jl b/scripts/Embedded/2d_thermal.jl
index 5b6f8f90..c0a21ce1 100644
--- a/scripts/Embedded/2d_thermal.jl
+++ b/scripts/Embedded/2d_thermal.jl
@@ -4,7 +4,7 @@ using Gridap,GridapTopOpt
 include("embedded_measures.jl")
 
 function main()
-  path="./results/CellFEM_thermal_compliance_ALM/"
+  path="./results/UnfittedFEM_thermal_compliance_ALM/"
   n = 200
   order = 1
   γ = 0.1
diff --git a/scripts/Embedded/2d_thermal_optimised.jl b/scripts/Embedded/2d_thermal_optimised.jl
index f456c6e1..24dc88c4 100644
--- a/scripts/Embedded/2d_thermal_optimised.jl
+++ b/scripts/Embedded/2d_thermal_optimised.jl
@@ -1,10 +1,10 @@
 using Pkg; Pkg.activate()
 
-using Gridap,GridapTopOpt
+using Gridap,GridapTopOpt,GridapEmbedded
 include("embedded_measures.jl")
 
 function main()
-  path="./results/CellFEM_thermal_compliance_ALM_Optimised/"
+  path="./results/UnfittedFEM_thermal_compliance_ALM_Optimised_test/"
   n = 200
   order = 1
   γ = 0.1
@@ -76,7 +76,7 @@ function main()
   ## Optimiser
   rm(path,force=true,recursive=true)
   mkpath(path)
-  optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;
+  optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;reinit_mod=5,
     γ,γ_reinit,verbose=true,constraint_names=[:Vol])
   for (it,uh,φh) in optimiser
     dΩ1,_,dΓ = get_meas(φh)
diff --git a/scripts/Embedded/2d_thermal_optimised_MPI.jl b/scripts/Embedded/2d_thermal_optimised_MPI.jl
new file mode 100644
index 00000000..47b77839
--- /dev/null
+++ b/scripts/Embedded/2d_thermal_optimised_MPI.jl
@@ -0,0 +1,100 @@
+using Gridap,GridapTopOpt,GridapEmbedded
+using GridapDistributed, GridapPETSc, PartitionedArrays
+
+include("embedded_measures.jl")
+
+function main(parts,distribute)
+  ranks = distribute(LinearIndices((prod(parts),)))
+
+  path="./results/UnfittedFEM_thermal_compliance_ALM_Optimised_MPI_test/"
+  n = 200
+  order = 1
+  γ = 0.1
+  γ_reinit = 0.5
+  max_steps = floor(Int,order*minimum(n)/10)
+  tol = 1/(5*order^2)/minimum(n)
+  κ = 1
+  vf = 0.4
+  α_coeff = 4max_steps*γ
+  iter_mod = 1
+  i_am_main(ranks) && rm(path,force=true,recursive=true)
+  i_am_main(ranks) && mkpath(path)
+
+  model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),(n,n));
+  el_Δ = get_el_Δ(model)
+  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())
+  update_labels!(1,model,f_Γ_D,"Gamma_D")
+  update_labels!(2,model,f_Γ_N,"Gamma_N")
+
+  ## 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Ω)
+
+  ## Spaces
+  reffe_scalar = ReferenceFE(lagrangian,Float64,order)
+  V = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_D"])
+  U = TrialFESpace(V,0.0)
+  V_φ = TestFESpace(model,reffe_scalar)
+  V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_N"])
+  U_reg = TrialFESpace(V_reg,0)
+
+  φh = interpolate(initial_lsf(8,0.2),V_φ)
+  embedded_meas = EmbeddedMeasureCache(φh,V_φ)
+  update_meas(φ) = update_embedded_measures!(φ,embedded_meas)
+  get_meas(φ) = get_embedded_measures(φ,embedded_meas)
+
+  ## Weak form
+  a(u,v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(v))dΩ1 + ∫(10^-6*∇(u)⋅∇(v))dΩ2
+  l(v,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(v)dΓ_N
+
+  ## Optimisation functionals
+  J(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(u))dΩ1
+  dJ(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(∇(u)⋅∇(u)*q)dΓ;
+  Vol(u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(1/vol_D)dΩ1 - ∫(vf)dΩ;
+  dVol(q,u,φ,dΩ,dΓ_N,dΩ1,dΩ2,dΓ) = ∫(-1/vol_D*q)dΓ
+
+  ## IntegrandWithEmbeddedMeasure
+  a_iem = IntegrandWithEmbeddedMeasure(a,(dΩ,dΓ_N),update_meas)
+  l_iem = IntegrandWithEmbeddedMeasure(l,(dΩ,dΓ_N),get_meas)
+  J_iem = IntegrandWithEmbeddedMeasure(J,(dΩ,dΓ_N),get_meas)
+  dJ_iem = IntegrandWithEmbeddedMeasure(dJ,(dΩ,dΓ_N),get_meas)
+  Vol_iem = IntegrandWithEmbeddedMeasure(Vol,(dΩ,dΓ_N),get_meas)
+  dVol_iem = IntegrandWithEmbeddedMeasure(dVol,(dΩ,dΓ_N),get_meas)
+
+  ## Evolution Method
+  ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
+
+  ## Setup solver and FE operators
+  state_map = AffineFEStateMap(a_iem,l_iem,U,V,V_φ,U_reg,φh,(dΩ,dΓ_N))
+  pcfs = PDEConstrainedFunctionals(J_iem,[Vol_iem],state_map,analytic_dJ=dJ_iem,analytic_dC=[dVol_iem])
+
+  ## Hilbertian extension-regularisation problems
+  α = α_coeff*maximum(el_Δ)
+  a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ;
+  vel_ext = VelocityExtension(a_hilb,U_reg,V_reg)
+
+  ## Optimiser
+  optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh;reinit_mod=5,
+    γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
+  for (it,uh,φh) in optimiser
+    _Ω1,_,_Γ = get_embedded_triangulations(φh,embedded_meas)
+    if iszero(it % iter_mod)
+      writevtk(_Ω1,path*"Omega_out$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+      writevtk(_Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(_Γ)])
+    end
+    write_history(path*"/history.txt",optimiser.history;ranks)
+  end
+  it = get_history(optimiser).niter; uh = get_state(pcfs)
+  _Ω1,_,_Γ = get_embedded_triangulations(φh,embedded_meas)
+  writevtk(_Ω1,path*"Omega_out$it",cellfields=["φ"=>φh,"|∇(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
+  writevtk(_Γ,path*"Gamma_out$it",cellfields=["normal"=>get_normal_vector(_Γ)])
+end
+
+with_mpi() do distribute
+  parts = (3,3);
+  main(parts,distribute)
+end
\ No newline at end of file
diff --git a/scripts/Embedded/embedded_measures.jl b/scripts/Embedded/embedded_measures.jl
index d7ef4773..1de0c874 100644
--- a/scripts/Embedded/embedded_measures.jl
+++ b/scripts/Embedded/embedded_measures.jl
@@ -1,49 +1,80 @@
-using GridapEmbedded
+using GridapTopOpt: to_parray_of_arrays
+using GridapDistributed
+using GridapDistributed: DistributedDiscreteModel, DistributedDomainContribution
 using GridapEmbedded.LevelSetCutters
-using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData
+# using GridapEmbedded.LevelSetCutters: DiscreteGeometry
+using GridapEmbedded.Distributed: DistributedDiscreteGeometry
+using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData, Gridap.Helpers
 import Gridap.Geometry: get_node_coordinates, collect1d
 
-function cv_to_dof(cv,V)
-	fv=zeros(eltype(eltype(cv)),num_free_dofs(V))
-	gather_free_values!(fv,V,cv)
-end
-
-function field_to_cv(uh::FEFunction)
-	get_cell_dof_values(uh)
-end
-
-function field_to_cv(cf::CellField)
-	cv=cf.cell_field.args[1]
-end
-
-function get_geo_params(ϕₕ::FEFunction,Vbg)
+# function cv_to_dof(cv,V)
+# 	fv=zeros(eltype(eltype(cv)),num_free_dofs(V))
+# 	gather_free_values!(fv,V,cv)
+# end
+
+# function field_to_cv(uh::FEFunction)
+# 	get_cell_dof_values(uh)
+# end
+
+# function field_to_cv(cf::CellField)
+# 	cf.cell_field.args[1]
+# end
+
+# function get_geo_params(φh::FEFunction,Vbg)
+#   Ωbg = get_triangulation(Vbg)
+#   bgmodel = get_background_model(Ωbg)
+#   point_to_coords = collect1d(get_node_coordinates(bgmodel))
+#   ls_to_point_to_value_unmasked = field_to_cv(φh)
+#   p0 = cv_to_dof(ls_to_point_to_value_unmasked,Vbg)
+#   geo1 = DiscreteGeometry(p0,point_to_coords)
+#   geo2 = DiscreteGeometry(-1*p0,point_to_coords,name="")
+#   get_geo_params(geo1,geo2,bgmodel)
+# end
+
+# function get_geo_params(φh::CellField,Vbg)
+#   Ωbg = get_triangulation(Vbg)
+#   bgmodel = get_background_model(Ωbg)
+#   point_to_coords = collect1d(get_node_coordinates(bgmodel))
+#   ls_to_point_to_value_unmasked = field_to_cv(φh)
+#   p0 = cv_to_dof(ls_to_point_to_value_unmasked,Vbg)
+#   geo1 = DiscreteGeometry(p0,point_to_coords)
+#   geo2 = DiscreteGeometry(-1*p0,point_to_coords,name="")
+#   get_geo_params(geo1,geo2,bgmodel)
+# end
+
+# function get_geo_params(φ::AbstractVector,Vbg)
+# 	Ωbg = get_triangulation(Vbg)
+# 	bgmodel = get_background_model(Ωbg)
+# 	point_to_coords = collect1d(get_node_coordinates(bgmodel))
+# 	geo1 = DiscreteGeometry(φ,point_to_coords,name="")
+# 	geo2 = DiscreteGeometry(-φ,point_to_coords,name="")
+# 	get_geo_params(geo1,geo2,bgmodel)
+# end
+
+# The above is more efficent for serial problems but does not work for periodic problems or distirbuted
+#  problems. This is subject to change.
+
+_DiscreteGeometry(φ,model::DistributedDiscreteModel;name::String="") =
+  DistributedDiscreteGeometry(φ,model;name)
+
+_DiscreteGeometry(φh::CellField,model::CartesianDiscreteModel;name::String="") =
+  DiscreteGeometry(φh(collect1d(get_node_coordinates(model))),collect1d(get_node_coordinates(model));name)
+
+function get_geo_params(ϕh::CellField,Vbg)#::GridapDistributed.DistributedFESpace)
   Ωbg = get_triangulation(Vbg)
   bgmodel = get_background_model(Ωbg)
-  point_to_coords = collect1d(get_node_coordinates(bgmodel))
-  ls_to_point_to_value_unmasked = field_to_cv(ϕₕ)
-  p0 = cv_to_dof(ls_to_point_to_value_unmasked,Vbg)
-  geo1 = DiscreteGeometry(p0,point_to_coords)
-  geo2 = DiscreteGeometry(-1*p0,point_to_coords,name="")
+  ϕhminus = FEFunction(Vbg,-get_free_dof_values(ϕh))
+  geo1 = _DiscreteGeometry(ϕh,bgmodel)
+  geo2 = _DiscreteGeometry(ϕhminus,bgmodel,name="")
   get_geo_params(geo1,geo2,bgmodel)
 end
 
-function get_geo_params(ϕₕ::CellField,Vbg)
-  Ωbg = get_triangulation(Vbg)
-  bgmodel = get_background_model(Ωbg)
-  point_to_coords = collect1d(get_node_coordinates(bgmodel))
-  ls_to_point_to_value_unmasked = field_to_cv(ϕₕ)
-  p0 = cv_to_dof(ls_to_point_to_value_unmasked,Vbg)
-  geo1 = DiscreteGeometry(p0,point_to_coords)
-  geo2 = DiscreteGeometry(-1*p0,point_to_coords,name="")
-  get_geo_params(geo1,geo2,bgmodel)
-end
-
-function get_geo_params(ϕ::AbstractVector,Vbg)
+function get_geo_params(φ::AbstractVector,Vbg)#::GridapDistributed.DistributedFESpace)
 	Ωbg = get_triangulation(Vbg)
 	bgmodel = get_background_model(Ωbg)
-	point_to_coords = collect1d(get_node_coordinates(bgmodel))
-	geo1 = DiscreteGeometry(ϕ,point_to_coords,name="")
-	geo2 = DiscreteGeometry(-ϕ,point_to_coords,name="")
+	ϕh = FEFunction(Vbg,φ); ϕhminus = FEFunction(Vbg,-φ)
+	geo1 = _DiscreteGeometry(ϕh,bgmodel,name="")
+	geo2 = _DiscreteGeometry(ϕhminus,bgmodel,name="")
 	get_geo_params(geo1,geo2,bgmodel)
 end
 
@@ -60,7 +91,7 @@ function get_geo_params(geo1,geo2,bgmodel)
 	dΩ1 = Measure(Ω1,degree)
 	dΩ2 = Measure(Ω2,degree)
 	dΓ = Measure(Γ,degree)
-	(dΩ1,dΩ2,dΓ)
+	(Ω1,Ω2,Γ),(dΩ1,dΩ2,dΓ)
 end
 
 #######################
@@ -68,6 +99,7 @@ end
 # For problems where the derivatives are known, we only want to update measures once
 mutable struct EmbeddedMeasureCache
   const space
+  trians
   measures
 
   function EmbeddedMeasureCache(φ,space)
@@ -77,7 +109,9 @@ mutable struct EmbeddedMeasureCache
 end
 
 function update_embedded_measures!(φ,s::EmbeddedMeasureCache)
-  s.measures = get_geo_params(φ,s.space)
+  trians, measures = get_geo_params(φ,s.space)
+  s.measures = measures
+  s.trians = trians
   return s.measures
 end
 
@@ -85,6 +119,10 @@ function get_embedded_measures(φ,s::EmbeddedMeasureCache)
   return s.measures
 end
 
+function get_embedded_triangulations(φ,s::EmbeddedMeasureCache)
+  return s.trians
+end
+
 #######################
 
 import GridapTopOpt: AbstractIntegrandWithMeasure
@@ -98,14 +136,43 @@ end
 (F::IntegrandWithEmbeddedMeasure)(args...) = F.F(args...,F.dΩ...,F.get_embedded_dΩ(args[end])...)
 
 function Gridap.gradient(F::IntegrandWithEmbeddedMeasure,uh::Vector{<:FEFunction},K::Int)
-  # @check 0 < K <= length(uh)
+  @warn "Automatic differentation is currently disabled for `IntegrandWithEmbeddedMeasure` types" maxlog=1
+  @check 0 < K <= length(uh)
   _f(uk) = if K == length(uh)
     F.F(uh[1:K-1]...,uk,uh[K+1:end]...,F.dΩ...,F.get_embedded_dΩ(uk)...)
   else
     F.F(uh[1:K-1]...,uk,uh[K+1:end]...,F.dΩ...,F.get_embedded_dΩ(uh[end])...)
   end
-  # return Gridap.gradient(_f,uh[K])
-  return Gridap.gradient(u -> ∑(∫(0)F.dΩ[i] for i = 1:length(F.dΩ)),uh[K]) # AD is currently disabled
+  # return Gridap.gradient(_f,uh[K]) # AD is currently disabled due to error (under investigation)
+  return Gridap.gradient(u -> ∑(∫(0)F.dΩ[i] for i = 1:length(F.dΩ)),uh[K])
+end
+
+# This doesn't currently work, we need a nice way to differentiate local_embedded_measures
+#  like the above.
+function Gridap.gradient(F::IntegrandWithEmbeddedMeasure,uh::Vector,K::Int)
+  @warn "Automatic differentation is currently disabled for `IntegrandWithEmbeddedMeasure` types" maxlog=1
+  @check 0 < K <= length(uh)
+  local_fields = map(local_views,uh) |> to_parray_of_arrays
+  local_measures = map(local_views,F.dΩ) |> to_parray_of_arrays
+
+  # if K == length(uh)
+  #   # Not sure how to do this just yet...
+  # else
+  #   local_embedded_measures = map(local_views,F.get_embedded_dΩ(uh[end])) |> to_parray_of_arrays
+  #   contribs = map(local_measures,local_embedded_measures,local_fields) do dΩ,dΩ_embedded,lf
+  #     _f = u -> F.F(lf[1:K-1]...,u,lf[K+1:end]...,dΩ...,dΩ_embedded...)
+  #     return Gridap.Fields.gradient(_f,lf[K])
+  #   end
+  # end
+
+  # Placeholder
+  local_embedded_measures = map(local_views,F.get_embedded_dΩ(uh[end])) |> to_parray_of_arrays
+  contribs = map(local_measures,local_embedded_measures,local_fields) do dΩ,dΩ_embedded,lf
+    _f = u -> ∑(∫(0)dΩ[i] for i = 1:length(dΩ))
+    return Gridap.Fields.gradient(_f,lf[K])
+  end
+
+  return DistributedDomainContribution(contribs)
 end
 
 Gridap.gradient(F::IntegrandWithEmbeddedMeasure,uh) = Gridap.gradient(F,[uh],1)
\ No newline at end of file