diff --git a/scripts/Benchmarks/benchmark.jl b/scripts/Benchmarks/benchmark.jl
index c31d395..13900d6 100644
--- a/scripts/Benchmarks/benchmark.jl
+++ b/scripts/Benchmarks/benchmark.jl
@@ -76,23 +76,23 @@ function nl_elast(mesh_partition,ranks,el_size,order,verbose)
   ## Volume change
   J(F) = sqrt(det(C(F)))
   ## Residual
-  res(u,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*((dE ∘ (∇(v),∇(u))) ⊙ (S ∘ ∇(u))) )*dΩ - ∫(g⋅v)dΓ_N
+  res(u,v,φ) = ∫( (I ∘ φ)*((dE ∘ (∇(v),∇(u))) ⊙ (S ∘ ∇(u))) )*dΩ - ∫(g⋅v)dΓ_N
   Tm = SparseMatrixCSR{0,PetscScalar,PetscInt}
   Tv = Vector{PetscScalar}
   lin_solver = ElasticitySolver(V)
   nl_solver = NewtonSolver(lin_solver;maxiter=50,rtol=10^-8,verbose)
   state_map = NonlinearFEStateMap(
-    res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N;
+    res,U,V,V_φ,U_reg,φh;
     assem_U = SparseMatrixAssembler(Tm,Tv,U,V),
     assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U),
     assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg),
     nls = nl_solver, adjoint_ls = lin_solver
   )
   # Objective and constraints
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((dE ∘ (∇(u),∇(u))) ⊙ (S ∘ ∇(u))))dΩ
+  J(u,φ) = ∫((I ∘ φ)*((dE ∘ (∇(u),∇(u))) ⊙ (S ∘ ∇(u))))dΩ
   vol_D = sum(∫(1)dΩ)
-  C1(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
-  dC1(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  C1(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+  dC1(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
   pcfs = PDEConstrainedFunctionals(J,[C1],state_map,analytic_dC=[dC1])
   # Velocity extension
   α = 4max_steps*γ*maximum(get_el_Δ(model))
@@ -152,24 +152,24 @@ function therm(mesh_partition,ranks,el_size,order,verbose)
   interp = SmoothErsatzMaterialInterpolation(η = 2*maximum(get_el_Δ(model)))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
   # Weak formulation
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(g*v)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
+  l(v,φ) = ∫(g*v)dΓ_N
   Tm = SparseMatrixCSR{0,PetscScalar,PetscInt}
   Tv = Vector{PetscScalar}
   solver = PETScLinearSolver()
   state_map = AffineFEStateMap(
-    a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N;
+    a,l,U,V,V_φ,U_reg,φh;
     assem_U = SparseMatrixAssembler(Tm,Tv,U,V),
     assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U),
     assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg),
     ls = solver,adjoint_ls = solver
   )
   # Objective and constraints
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
+  dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
   vol_D = sum(∫(1)dΩ)
-  C1(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
-  dC1(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  C1(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+  dC1(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
   pcfs = PDEConstrainedFunctionals(J,[C1],state_map,
     analytic_dJ=dJ,analytic_dC=[dC1])
   # Velocity extension
@@ -228,24 +228,24 @@ function elast(mesh_partition,ranks,el_size,order,verbose)
   interp = SmoothErsatzMaterialInterpolation(η = 2*maximum(get_el_Δ(model)))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
   # Weak formulation
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v⋅g)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ
+  l(v,φ) = ∫(v⋅g)dΓ_N
   Tm = SparseMatrixCSR{0,PetscScalar,PetscInt}
   Tv = Vector{PetscScalar}
   solver = ElasticitySolver(V)
   state_map = AffineFEStateMap(
-    a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N;
+    a,l,U,V,V_φ,U_reg,φh;
     assem_U = SparseMatrixAssembler(Tm,Tv,U,V),
     assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U),
     assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg),
     ls = solver,adjoint_ls = solver
   )
   # Objective and constraints
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ
+  dJ(q,u,φ) = ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
   vol_D = sum(∫(1)dΩ)
-  C1(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
-  dC1(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  C1(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+  dC1(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
   pcfs = PDEConstrainedFunctionals(J,[C1],state_map,
     analytic_dJ=dJ,analytic_dC=[dC1])
   # Velocity extension
@@ -320,15 +320,15 @@ function inverter_HPM(mesh_partition,ranks,el_size,order,verbose)
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_in,dΓ_out) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out
-  l(v,φ,dΩ,dΓ_in,dΓ_out) = ∫(v⋅g)dΓ_in
+  a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out
+  l(v,φ) = ∫(v⋅g)dΓ_in
 
   ## Optimisation functionals
   e₁ = VectorValue(1,0,0)
-  J(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in
-  Vol(u,φ,dΩ,dΓ_in,dΓ_out) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_in,dΓ_out) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
-  UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out
+  J(u,φ) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  UΓ_out(u,φ) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out
 
   ## Finite difference solver
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
@@ -339,7 +339,7 @@ function inverter_HPM(mesh_partition,ranks,el_size,order,verbose)
   solver = ElasticitySolver(V)
 
   state_map = AffineFEStateMap(
-    a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_in,dΓ_out;
+    a,l,U,V,V_φ,U_reg,φh;
     assem_U = SparseMatrixAssembler(Tm,Tv,U,V),
     assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U),
     assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg),
diff --git a/scripts/MPI/elastic_compliance_ALM.jl b/scripts/MPI/elastic_compliance_ALM.jl
index 2cdafcc..ee18819 100644
--- a/scripts/MPI/elastic_compliance_ALM.jl
+++ b/scripts/MPI/elastic_compliance_ALM.jl
@@ -65,14 +65,14 @@ function main(mesh_partition,distribute,el_size,path)
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v⋅g)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ
+  l(v,φ) = ∫(v⋅g)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ
+  dJ(q,u,φ) = ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
@@ -83,7 +83,7 @@ function main(mesh_partition,distribute,el_size,path)
   solver = ElasticitySolver(V)
 
   state_map = AffineFEStateMap(
-    a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N;
+    a,l,U,V,V_φ,U_reg,φh;
     assem_U = SparseMatrixAssembler(Tm,Tv,U,V),
     assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U),
     assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg),
diff --git a/scripts/MPI/inverse_homenisation_ALM.jl b/scripts/MPI/inverse_homenisation_ALM.jl
index 5946f58..0e36f2f 100644
--- a/scripts/MPI/inverse_homenisation_ALM.jl
+++ b/scripts/MPI/inverse_homenisation_ALM.jl
@@ -62,16 +62,16 @@ function main(mesh_partition,distribute,el_size,path)
         TensorValue(0.,0.,0.,1.),     # ϵᵢⱼ⁽²²⁾≡ϵᵢⱼ⁽²⁾
         TensorValue(0.,1/2,1/2,0.))   # ϵᵢⱼ⁽¹²⁾≡ϵᵢⱼ⁽³⁾
 
-  a(u,v,φ,dΩ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v))dΩ
-  l = [(v,φ,dΩ) -> ∫(-(I ∘ φ) * C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3]
+  a(u,v,φ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v))dΩ
+  l = [(v,φ) -> ∫(-(I ∘ φ) * C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3]
 
   ## Optimisation functionals
-  Cᴴ(r,s,u,φ,dΩ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ
-  dCᴴ(r,s,q,u,φ,dΩ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
-  κ(u,φ,dΩ) = -1/4*(Cᴴ(1,1,u,φ,dΩ)+Cᴴ(2,2,u,φ,dΩ)+2*Cᴴ(1,2,u,φ,dΩ))
-  dκ(q,u,φ,dΩ) = -1/4*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+2*dCᴴ(1,2,q,u,φ,dΩ))
-  Vol(u,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  Cᴴ(r,s,u,φ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ
+  dCᴴ(r,s,q,u,φ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  κ(u,φ) = -1/4*(Cᴴ(1,1,u,φ)+Cᴴ(2,2,u,φ)+2*Cᴴ(1,2,u,φ))
+  dκ(q,u,φ) = -1/4*(dCᴴ(1,1,q,u,φ)+dCᴴ(2,2,q,u,φ)+2*dCᴴ(1,2,q,u,φ))
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
@@ -82,7 +82,7 @@ function main(mesh_partition,distribute,el_size,path)
   solver = ElasticitySolver(V)
 
   state_map = RepeatingAffineFEStateMap(
-    3,a,l,U,V,V_φ,U_reg,φh,dΩ;
+    3,a,l,U,V,V_φ,U_reg,φh;
     assem_U = SparseMatrixAssembler(Tm,Tv,U,V),
     assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U),
     assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg),
diff --git a/scripts/MPI/nonlinear_thermal_compliance_ALM.jl b/scripts/MPI/nonlinear_thermal_compliance_ALM.jl
index beccd0e..0991d79 100644
--- a/scripts/MPI/nonlinear_thermal_compliance_ALM.jl
+++ b/scripts/MPI/nonlinear_thermal_compliance_ALM.jl
@@ -68,18 +68,18 @@ function main(mesh_partition,distribute,el_size,path)
   κ0 = 1
   ξ = -1
   κ(u) = κ0*(exp(ξ*u))
-  res(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N
+  res(u,v,φ) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(u))dΩ
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(u))dΩ
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/scripts/MPI/thermal_compliance_ALM.jl b/scripts/MPI/thermal_compliance_ALM.jl
index 9c78fbb..1e9ea3d 100644
--- a/scripts/MPI/thermal_compliance_ALM.jl
+++ b/scripts/MPI/thermal_compliance_ALM.jl
@@ -63,20 +63,20 @@ function main(mesh_partition,distribute,path)
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
+  l(v,φ) = ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
+  dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/scripts/MPI/thermal_compliance_ALM_PETSc.jl b/scripts/MPI/thermal_compliance_ALM_PETSc.jl
index 9248276..c5fb28e 100644
--- a/scripts/MPI/thermal_compliance_ALM_PETSc.jl
+++ b/scripts/MPI/thermal_compliance_ALM_PETSc.jl
@@ -63,14 +63,14 @@ function main(mesh_partition,distribute,el_size,path)
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
+  l(v,φ) = ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
+  dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
@@ -81,7 +81,7 @@ function main(mesh_partition,distribute,el_size,path)
   solver = PETScLinearSolver()
 
   state_map = AffineFEStateMap(
-    a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N;
+    a,l,U,V,V_φ,U_reg,φh;
     assem_U = SparseMatrixAssembler(Tm,Tv,U,V),
     assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U),
     assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg),
diff --git a/scripts/Serial/elastic_compliance_ALM.jl b/scripts/Serial/elastic_compliance_ALM.jl
index 9f0b702..1e400ee 100644
--- a/scripts/Serial/elastic_compliance_ALM.jl
+++ b/scripts/Serial/elastic_compliance_ALM.jl
@@ -60,20 +60,20 @@ function main(path="./results/elastic_compliance_ALM/")
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v⋅g)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ
+  l(v,φ) = ∫(v⋅g)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ
+  dJ(q,u,φ) = ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/scripts/Serial/elastic_compliance_HPM.jl b/scripts/Serial/elastic_compliance_HPM.jl
index 9f787c3..eea2c88 100644
--- a/scripts/Serial/elastic_compliance_HPM.jl
+++ b/scripts/Serial/elastic_compliance_HPM.jl
@@ -61,20 +61,20 @@ function main(path="./results/elastic_compliance_HPM/")
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v⋅g)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ
+  l(v,φ) = ∫(v⋅g)dΓ_N
 
   ## Optimisation functionals
-  J = (u,φ,dΩ,dΓ_N) -> ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ
-  dJ = (q,u,φ,dΩ,dΓ_N) -> ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
-  Vol = (u,φ,dΩ,dΓ_N) -> ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol = (q,u,φ,dΩ,dΓ_N) -> ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J = (u,φ) -> ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ
+  dJ = (q,u,φ) -> ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  Vol = (u,φ) -> ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol = (q,u,φ) -> ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/scripts/Serial/elastic_compliance_LM.jl b/scripts/Serial/elastic_compliance_LM.jl
index b09e501..e4c090d 100644
--- a/scripts/Serial/elastic_compliance_LM.jl
+++ b/scripts/Serial/elastic_compliance_LM.jl
@@ -58,19 +58,19 @@ function main(path="./results/elastic_compliance_LM/")
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v⋅g)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ
+  l(v,φ) = ∫(v⋅g)dΓ_N
 
   ## Optimisation functionals
   ξ = 2
-  J = (u,φ,dΩ,dΓ_N) -> ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)) + ξ*(ρ ∘ φ))dΩ
-  dJ = (q,u,φ,dΩ,dΓ_N) -> ∫((- ξ + C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  J = (u,φ) -> ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)) + ξ*(ρ ∘ φ))dΩ
+  dJ = (q,u,φ) -> ∫((- ξ + C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,state_map,analytic_dJ=dJ)
 
   ## Hilbertian extension-regularisation problems
diff --git a/scripts/Serial/hyperelastic_compliance_ALM.jl b/scripts/Serial/hyperelastic_compliance_ALM.jl
index dce96bd..42ef224 100644
--- a/scripts/Serial/hyperelastic_compliance_ALM.jl
+++ b/scripts/Serial/hyperelastic_compliance_ALM.jl
@@ -62,18 +62,18 @@ function main(path="./results/hyperelastic_compliance_ALM/")
   E(F) = 0.5*( F' ⋅ F - one(F) )
   Σ(∇u) = λ*tr(E(F(∇u)))*one(∇u)+2*μ*E(F(∇u))
   T(∇u) = F(∇u) ⋅ Σ(∇u)
-  res(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((T ∘ ∇(u)) ⊙ ∇(v)))*dΩ - ∫(g⋅v)dΓ_N
+  res(u,v,φ) = ∫((I ∘ φ)*((T ∘ ∇(u)) ⊙ ∇(v)))*dΩ - ∫(g⋅v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((T ∘ ∇(u)) ⊙ ∇(u)))dΩ
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*((T ∘ ∇(u)) ⊙ ∇(u)))dΩ
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl b/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl
index eca9d1c..8927fcd 100644
--- a/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl
+++ b/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl
@@ -76,18 +76,18 @@ function main(path="./results/hyperelastic_compliance_neohook_ALM/")
 
   # Cauchy stress tensor and residual
   σ(∇u) = (1.0/J(F(∇u)))*F(∇u)⋅S(∇u)⋅(F(∇u))'
-  res(u,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*((dE∘(∇(v),∇(u))) ⊙ (S∘∇(u))) )*dΩ - ∫(g⋅v)dΓ_N
+  res(u,v,φ) = ∫( (I ∘ φ)*((dE∘(∇(v),∇(u))) ⊙ (S∘∇(u))) )*dΩ - ∫(g⋅v)dΓ_N
 
   ## Optimisation functionals
-  Obj(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((dE∘(∇(u),∇(u))) ⊙ (S∘∇(u))))dΩ
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  Obj(u,φ) = ∫((I ∘ φ)*((dE∘(∇(u),∇(u))) ⊙ (S∘∇(u))))dΩ
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(Obj,[Vol],state_map,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/scripts/Serial/inverse_homenisation_ALM.jl b/scripts/Serial/inverse_homenisation_ALM.jl
index 972d432..0689e11 100644
--- a/scripts/Serial/inverse_homenisation_ALM.jl
+++ b/scripts/Serial/inverse_homenisation_ALM.jl
@@ -58,22 +58,22 @@ function main(path="./results/inverse_homenisation_ALM/")
         TensorValue(0.,0.,0.,1.),     # ϵᵢⱼ⁽²²⁾≡ϵᵢⱼ⁽²⁾
         TensorValue(0.,1/2,1/2,0.))   # ϵᵢⱼ⁽¹²⁾≡ϵᵢⱼ⁽³⁾
 
-  a(u,v,φ,dΩ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v) )dΩ
-  l = [(v,φ,dΩ) -> ∫(-(I ∘ φ)* C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3]
+  a(u,v,φ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v) )dΩ
+  l = [(v,φ) -> ∫(-(I ∘ φ)* C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3]
 
   ## Optimisation functionals
-  Cᴴ(r,s,u,φ,dΩ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ
-  dCᴴ(r,s,q,u,φ,dΩ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
-  κ(u,φ,dΩ) = -1/4*(Cᴴ(1,1,u,φ,dΩ)+Cᴴ(2,2,u,φ,dΩ)+2*Cᴴ(1,2,u,φ,dΩ))
-  dκ(q,u,φ,dΩ) = -1/4*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+2*dCᴴ(1,2,q,u,φ,dΩ))
-  Vol(u,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
-  dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  Cᴴ(r,s,u,φ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ
+  dCᴴ(r,s,q,u,φ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  κ(u,φ) = -1/4*(Cᴴ(1,1,u,φ)+Cᴴ(2,2,u,φ)+2*Cᴴ(1,2,u,φ))
+  dκ(q,u,φ) = -1/4*(dCᴴ(1,1,q,u,φ)+dCᴴ(2,2,q,u,φ)+2*dCᴴ(1,2,q,u,φ))
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh,dΩ)
+  state_map = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/scripts/Serial/inverter_ALM.jl b/scripts/Serial/inverter_ALM.jl
index ad2d878..d1db290 100644
--- a/scripts/Serial/inverter_ALM.jl
+++ b/scripts/Serial/inverter_ALM.jl
@@ -74,21 +74,21 @@ function main(path="./results/inverter_ALM/")
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_in,dΓ_out) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out
-  l(v,φ,dΩ,dΓ_in,dΓ_out) = ∫(v⋅g)dΓ_in
+  a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out
+  l(v,φ) = ∫(v⋅g)dΓ_in
 
   ## Optimisation functionals
   e₁ = VectorValue(1,0)
-  J(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in
-  Vol(u,φ,dΩ,dΓ_in,dΓ_out) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
-  dVol(q,u,φ,dΩ,dΓ_in,dΓ_out) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
-  UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out
+  J(u,φ) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  UΓ_out(u,φ) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_in,dΓ_out)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol,UΓ_out],state_map,analytic_dC=[dVol,nothing])
 
   ## Hilbertian extension-regularisation problems
diff --git a/scripts/Serial/inverter_HPM.jl b/scripts/Serial/inverter_HPM.jl
index fbe1312..cc20629 100644
--- a/scripts/Serial/inverter_HPM.jl
+++ b/scripts/Serial/inverter_HPM.jl
@@ -75,21 +75,21 @@ function main(path="./results/inverter_HPM/")
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_in,dΓ_out) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out
-  l(v,φ,dΩ,dΓ_in,dΓ_out) = ∫(v⋅g)dΓ_in
+  a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out
+  l(v,φ) = ∫(v⋅g)dΓ_in
 
   ## Optimisation functionals
   e₁ = VectorValue(1,0)
-  J(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in
-  Vol(u,φ,dΩ,dΓ_in,dΓ_out) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
-  dVol(q,u,φ,dΩ,dΓ_in,dΓ_out) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
-  UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out
+  J(u,φ) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  UΓ_out(u,φ) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_in,dΓ_out)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol,UΓ_out],state_map,analytic_dC=[dVol,nothing])
 
   ## Hilbertian extension-regularisation problems
diff --git a/scripts/Serial/nonlinear_thermal_compliance_ALM.jl b/scripts/Serial/nonlinear_thermal_compliance_ALM.jl
index 2bf3248..b69a31b 100644
--- a/scripts/Serial/nonlinear_thermal_compliance_ALM.jl
+++ b/scripts/Serial/nonlinear_thermal_compliance_ALM.jl
@@ -65,18 +65,18 @@ function main(path="./results/nonlinear_thermal_compliance_ALM/")
   κ0 = 1
   ξ = -1
   κ(u) = κ0*(exp(ξ*u))
-  res(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N
+  res(u,v,φ) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(u))dΩ
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(u))dΩ
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/scripts/Serial/thermal_compliance_ALM.jl b/scripts/Serial/thermal_compliance_ALM.jl
index aaf45b2..a013cfa 100644
--- a/scripts/Serial/thermal_compliance_ALM.jl
+++ b/scripts/Serial/thermal_compliance_ALM.jl
@@ -61,20 +61,20 @@ function main(path="./results/thermal_compliance_ALM/")
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
+  l(v,φ) = ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
+  dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl b/scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl
index 51f8e1a..4720998 100644
--- a/scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl
+++ b/scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl
@@ -67,20 +67,20 @@ function main(path="./results/thermal_compliance_ALM_higherorderlsf/")
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
+  l(v,φ) = ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
+  dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/src/ChainRules.jl b/src/ChainRules.jl
index 4661c6b..bcedbab 100644
--- a/src/ChainRules.jl
+++ b/src/ChainRules.jl
@@ -1,93 +1,38 @@
-abstract type AbstractIntegrandWithMeasure end
-
-"""
-    struct IntegrandWithMeasure{A,B<:Tuple}
-
-A wrapper to enable serial or parallel partial differentation of an
-integral `F` using `Gridap.gradient`. This is required to allow automatic
-differentation with `DistributedMeasure`.
-
-# Properties
-- `F  :: A`: A function that returns a `DomainContribution` or `DistributedDomainContribution`.
-- `dΩ :: B`: A tuple of measures.
-"""
-struct IntegrandWithMeasure{A,B<:Tuple} <: AbstractIntegrandWithMeasure
-  F  :: A
-  dΩ :: B
-end
-
-"""
-    (F::IntegrandWithMeasure)(args...)
-
-Evaluate `F.F` given arguments `args`.
 """
-(F::IntegrandWithMeasure)(args...) = F.F(args...,F.dΩ...)
+    Gridap.gradient(F,uh::Vector,K::Int)
 
-"""
-    Gridap.gradient(F::IntegrandWithMeasure,uh::Vector,K::Int)
-
-Given an an `IntegrandWithMeasure` `F` and a vector of `FEFunctions` `uh` (excluding measures)
-evaluate the partial derivative of `F.F` with respect to `uh[K]`.
+Given a function `F` that returns a DomainContribution when called, and a vector of 
+`FEFunctions` `uh`, evaluate the partial derivative of `F` with respect to `uh[K]`.
 
 # Example
 
 Suppose `uh` and `φh` are FEFunctions with measures `dΩ` and `dΓ_N`.
 Then the partial derivative of a function `J` wrt to `φh` is computed via
 ````
-J(u,φ,dΩ,dΓ_N) = ∫(f(u,φ))dΩ + ∫(g(u,φ))dΓ_N
-J_iwm = IntegrandWithMeasure(J,(dΩ,dΓ_N))
-∂J∂φh = ∇(J_iwm,[uh,φh],2)
+J(u,φ) = ∫(f(u,φ))dΩ + ∫(g(u,φ))dΓ_N
+∂J∂φh = ∇(J,[uh,φh],2)
 ````
 where `f` and `g` are user defined.
 """
-function Gridap.gradient(F::IntegrandWithMeasure,uh::Vector{<:FEFunction},K::Int)
+function Gridap.gradient(F,uh::Vector{<:CellField},K::Int)
   @check 0 < K <= length(uh)
-  _f(uk) = F.F(uh[1:K-1]...,uk,uh[K+1:end]...,F.dΩ...)
+  _f(uk) = F(uh[1:K-1]...,uk,uh[K+1:end]...)
   return Gridap.gradient(_f,uh[K])
 end
 
-function Gridap.gradient(F::IntegrandWithMeasure,uh::Vector,K::Int)
-  @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
-  contribs = map(local_measures,local_fields) do dΩ,lf
-    # TODO: Remove second term below, this is a fix for the problem discussed in
-    #  https://github.com/zjwegert/GridapTopOpt.jl/issues/46
-    _f = u -> F.F(lf[1:K-1]...,u,lf[K+1:end]...,dΩ...) #+ ∑(∫(0)dΩ[i] for i = 1:length(dΩ))
-    return Gridap.Fields.gradient(_f,lf[K])
-  end
-  return DistributedDomainContribution(contribs)
-end
-
-Gridap.gradient(F::IntegrandWithMeasure,uh) = Gridap.gradient(F,[uh],1)
-
 """
-    Gridap.jacobian(F::IntegrandWithMeasure,uh::Vector,K::Int)
+    Gridap.jacobian(F,uh::Vector,K::Int)
 
-Given an an `IntegrandWithMeasure` `F` and a vector of `FEFunctions` or `CellField` `uh`
-(excluding measures) evaluate the Jacobian `F.F` with respect to `uh[K]`.
+Given a function `F` that returns a DomainContribution when called, and a 
+vector of `FEFunctions` or `CellField` `uh`, evaluate the Jacobian
+`F` with respect to `uh[K]`.
 """
-function Gridap.jacobian(F::IntegrandWithMeasure,uh::Vector{<:Union{FEFunction,CellField}},K::Int)
+function Gridap.jacobian(F,uh::Vector{<:CellField},K::Int)
   @check 0 < K <= length(uh)
-  _f(uk) = F.F(uh[1:K-1]...,uk,uh[K+1:end]...,F.dΩ...)
+  _f(uk) = F(uh[1:K-1]...,uk,uh[K+1:end]...)
   return Gridap.jacobian(_f,uh[K])
 end
 
-DistributedFields = Union{DistributedCellField,DistributedMultiFieldCellField}
-
-function Gridap.jacobian(F::IntegrandWithMeasure,uh::Vector{<:DistributedFields},K::Int)
-  @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
-  contribs = map(local_measures,local_fields) do dΩ,lf
-    _f = u -> F.F(lf[1:K-1]...,u,lf[K+1:end]...,dΩ...)
-    return Gridap.jacobian(_f,lf[K])
-  end
-  return DistributedDomainContribution(contribs)
-end
-
-Gridap.jacobian(F::IntegrandWithMeasure,uh) = Gridap.jacobian(F,[uh],1)
-
 function GridapDistributed.to_parray_of_arrays(a::NTuple{N,T}) where {N,T<:DebugArray}
   indices = linear_indices(first(a))
   map(indices) do i
@@ -107,21 +52,21 @@ function GridapDistributed.to_parray_of_arrays(a::NTuple{N,T}) where {N,T<:MPIAr
 end
 
 """
-    struct StateParamIntegrandWithMeasure{A<:IntegrandWithMeasure,B,C,D}
+    struct StateParamIntegrandWithMeasure{A,B,C,D}
 
-A wrapper to handle partial differentation of an [`IntegrandWithMeasure`](@ref)
+A wrapper to handle partial differentation of a function F
 of a specific form (see below) in a `ChainRules.jl` compatible way with caching.
 
 # Assumptions
 
-We assume that we have a `IntegrandWithMeasure` of the following form:
+We assume that we have a function F of the following form:
 
-`F(u,φ,dΩ₁,dΩ₂,...) = ∫(f(u,φ))dΩ₁ + ∫(g(u,φ))dΩ₂ + ...,`.
+`F(u,φ) = ∫(f(u,φ))dΩ₁ + ∫(g(u,φ))dΩ₂ + ...,`.
 
 where `u` and `φ` are each expected to inherit from `Union{FEFunction,MultiFieldFEFunction}`
 or the GridapDistributed equivalent.
 """
-struct StateParamIntegrandWithMeasure{A<:AbstractIntegrandWithMeasure,B,C,D}
+struct StateParamIntegrandWithMeasure{A,B,C,D}
   F       :: A
   spaces  :: B
   assems  :: C
@@ -129,17 +74,15 @@ struct StateParamIntegrandWithMeasure{A<:AbstractIntegrandWithMeasure,B,C,D}
 end
 
 """
-    StateParamIntegrandWithMeasure(F::IntegrandWithMeasure,U::FESpace,V_φ::FESpace,
+    StateParamIntegrandWithMeasure(F,U::FESpace,V_φ::FESpace,
     U_reg::FESpace,assem_U::Assembler,assem_deriv::Assembler)
 
 Create an instance of `StateParamIntegrandWithMeasure`.
 """
 function StateParamIntegrandWithMeasure(
-  F::AbstractIntegrandWithMeasure,
-  U::FESpace,V_φ::FESpace,U_reg::FESpace,
+  F,U::FESpace,V_φ::FESpace,U_reg::FESpace,
   assem_U::Assembler,assem_deriv::Assembler
 )
-  # TODO: Find alternative to this later, is a fix for grad of embedded FEs
   φ₀, u₀ = interpolate(x->-sqrt((x[1]-1/2)^2+(x[2]-1/2)^2)+0.2,V_φ), zero(U)
   ∂j∂u_vecdata = collect_cell_vector(U,∇(F,[u₀,φ₀],1))
   ∂j∂φ_vecdata = collect_cell_vector(U_reg,∇(F,[u₀,φ₀],2))
@@ -215,13 +158,6 @@ Return the solution/state `u` to the FE problem.
 """
 get_state(::AbstractFEStateMap) = @abstractmethod
 
-"""
-    get_measure(m::AbstractFEStateMap)
-
-Return the measures associated with the FE problem.
-"""
-get_measure(::AbstractFEStateMap) = @abstractmethod
-
 """
     get_spaces(m::AbstractFEStateMap)
 
@@ -400,19 +336,7 @@ function ChainRulesCore.rrule(φ_to_u::AbstractFEStateMap,φ::AbstractVector)
   return ChainRulesCore.rrule(φ_to_u,φh)
 end
 
-"""
-    StateParamIntegrandWithMeasure(f::Function,φ_to_u::AbstractFEStateMap)
-
-Create an instance of `StateParamIntegrandWithMeasure` given a `f` and
-`φ_to_u`.
-"""
-function StateParamIntegrandWithMeasure(f::Function,φ_to_u::AbstractFEStateMap)
-  dΩ = get_measure(φ_to_u)
-  F  = IntegrandWithMeasure(f,dΩ)
-  StateParamIntegrandWithMeasure(F,φ_to_u)
-end
-
-function StateParamIntegrandWithMeasure(F::AbstractIntegrandWithMeasure,φ_to_u::AbstractFEStateMap)
+function StateParamIntegrandWithMeasure(F::Function,φ_to_u::AbstractFEStateMap)
   U,V,V_φ,U_reg = φ_to_u.spaces
   assem_deriv = get_deriv_assembler(φ_to_u)
   assem_U = get_pde_assembler(φ_to_u)
@@ -427,8 +351,8 @@ element operators `AffineFEOperator`.
 
 # Parameters
 
-- `biform::A`: `IntegrandWithMeasure` defining the bilinear form.
-- `liform::B`: `IntegrandWithMeasure` defining the linear form.
+- `biform::A`: `Function` defining the bilinear form.
+- `liform::B`: `Function` defining the linear form.
 - `spaces::C`: `Tuple` of finite element spaces.
 - `plb_caches::D`: A cache for the pullback operator.
 - `fwd_caches::E`: A cache for the forward problem.
@@ -445,7 +369,7 @@ struct AffineFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap
   @doc """
       AffineFEStateMap(
         a::Function,l::Function,
-        U,V,V_φ,U_reg,φh,dΩ...;
+        U,V,V_φ,U_reg,φh;
         assem_U = SparseMatrixAssembler(U,V),
         assem_adjoint = SparseMatrixAssembler(V,U),
         assem_deriv = SparseMatrixAssembler(U_reg,U_reg),
@@ -460,8 +384,8 @@ struct AffineFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap
   Optional arguments enable specification of assemblers and linear solvers.
   """
     function AffineFEStateMap(
-      biform::AbstractIntegrandWithMeasure,liform::AbstractIntegrandWithMeasure,
-      U,V,V_φ,U_reg,φh,dΩ...;
+      biform::Function,liform::Function,
+      U,V,V_φ,U_reg,φh;
       assem_U = SparseMatrixAssembler(U,V),
       assem_adjoint = SparseMatrixAssembler(V,U),
       assem_deriv = SparseMatrixAssembler(U_reg,U_reg),
@@ -497,15 +421,8 @@ struct AffineFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap
   end
 end
 
-function AffineFEStateMap(a::Function,l::Function,U,V,V_φ,U_reg,φh,dΩ...;kwargs...)
-  biform = IntegrandWithMeasure(a,dΩ)
-  liform = IntegrandWithMeasure(l,dΩ)
-  AffineFEStateMap(biform,liform,U,V,V_φ,U_reg,φh,dΩ;kwargs...)
-end
-
 # Getters
 get_state(m::AffineFEStateMap) = FEFunction(get_trial_space(m),m.fwd_caches[4])
-get_measure(m::AffineFEStateMap) = m.biform.dΩ
 get_spaces(m::AffineFEStateMap) = m.spaces
 get_assemblers(m::AffineFEStateMap) = (m.fwd_caches[6],m.plb_caches[2],m.adj_caches[4])
 
@@ -549,7 +466,7 @@ element operators.
 
 # Parameters
 
-- `res::A`: an `IntegrandWithMeasure` defining the residual of the problem.
+- `res::A`: a `Function` defining the residual of the problem.
 - `jac::B`: a `Function` defining Jacobian of the residual.
 - `spaces::C`: `Tuple` of finite element spaces.
 - `plb_caches::D`: A cache for the pullback operator.
@@ -566,7 +483,7 @@ struct NonlinearFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap
 
   @doc """
       NonlinearFEStateMap(
-        res::Function,U,V,V_φ,U_reg,φh,dΩ...;
+        res::Function,U,V,V_φ,U_reg,φh;
         assem_U = SparseMatrixAssembler(U,V),
         assem_adjoint = SparseMatrixAssembler(V,U),
         assem_deriv = SparseMatrixAssembler(U_reg,U_reg),
@@ -581,7 +498,7 @@ struct NonlinearFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap
   Optional arguments enable specification of assemblers, nonlinear solver, and adjoint (linear) solver.
   """
   function NonlinearFEStateMap(
-    res::AbstractIntegrandWithMeasure,jac,U,V,V_φ,U_reg,φh,dΩ...;
+    res::Function,jac::Function,U,V,V_φ,U_reg,φh;
     assem_U = SparseMatrixAssembler(U,V),
     assem_adjoint = SparseMatrixAssembler(V,U),
     assem_deriv = SparseMatrixAssembler(U_reg,U_reg),
@@ -616,23 +533,14 @@ struct NonlinearFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap
   end
 end
 
-function NonlinearFEStateMap(res::AbstractIntegrandWithMeasure,U,V,V_φ,U_reg,φh,dΩ...;jac::Nothing=nothing,kwargs...)
-  jacf = (u,du,v,φh) -> jacobian(res,[u,v,φh],1)
-  NonlinearFEStateMap(res,jacf,U,V,V_φ,U_reg,φh,dΩ;kwargs...)
-end
-
-function NonlinearFEStateMap(res::Function,U,V,V_φ,U_reg,φh,dΩ...;jac=nothing,kwargs...)
-  _res = IntegrandWithMeasure(res,dΩ)
+function NonlinearFEStateMap(res::Function,U,V,V_φ,U_reg,φh;jac=nothing,kwargs...)
   if isnothing(jac)
-    _jac = (u,du,v,φh) -> jacobian(_res,[u,v,φh],1)
-  else
-    _jac = IntegrandWithMeasure(jac,dΩ)
+    jac = (u,du,v,φh) -> jacobian(res,[u,v,φh],1)
   end
-  NonlinearFEStateMap(_res,_jac,U,V,V_φ,U_reg,φh,dΩ;kwargs...)
+  NonlinearFEStateMap(res,jac,U,V,V_φ,U_reg,φh;kwargs...)
 end
 
 get_state(m::NonlinearFEStateMap) = FEFunction(get_trial_space(m),m.fwd_caches[3])
-get_measure(m::NonlinearFEStateMap) = m.res.dΩ
 get_spaces(m::NonlinearFEStateMap) = m.spaces
 get_assemblers(m::NonlinearFEStateMap) = (m.fwd_caches[4],m.plb_caches[2],m.adj_caches[4])
 
@@ -676,8 +584,8 @@ a single bilinear form.
 
 # Parameters
 
-- `biform`: `IntegrandWithMeasure` defining the bilinear form.
-- `liform`: A vector of `IntegrandWithMeasure` defining the linear forms.
+- `biform`: A `Function` defining the bilinear form.
+- `liform`: A vector of `Function` defining the linear forms.
 - `spaces`: Repeated finite element spaces.
 - `spaces_0`: Original finite element spaces that are being repeated.
 - `plb_caches`: A cache for the pullback operator.
@@ -696,7 +604,7 @@ struct RepeatingAffineFEStateMap{A,B,C,D,E,F,G} <: AbstractFEStateMap
   @doc """
       RepeatingAffineFEStateMap(
         nblocks::Int,a::Function,l::Vector{<:Function},
-        U0,V0,V_φ,U_reg,φh,dΩ...;
+        U0,V0,V_φ,U_reg,φh;
         assem_U = SparseMatrixAssembler(U0,V0),
         assem_adjoint = SparseMatrixAssembler(V0,U0),
         assem_deriv = SparseMatrixAssembler(U_reg,U_reg),
@@ -717,8 +625,8 @@ struct RepeatingAffineFEStateMap{A,B,C,D,E,F,G} <: AbstractFEStateMap
     where each field corresponds to an entry in the vector of linear forms
   """
   function RepeatingAffineFEStateMap(
-    nblocks::Int,biform::AbstractIntegrandWithMeasure,liforms::Vector{<:AbstractIntegrandWithMeasure},
-    U0,V0,V_φ,U_reg,φh,dΩ...;
+    nblocks::Int,biform::Function,liforms::Vector{<:Function},
+    U0,V0,V_φ,U_reg,φh;
     assem_U = SparseMatrixAssembler(U0,V0),
     assem_adjoint = SparseMatrixAssembler(V0,U0),
     assem_deriv = SparseMatrixAssembler(U_reg,U_reg),
@@ -766,14 +674,6 @@ struct RepeatingAffineFEStateMap{A,B,C,D,E,F,G} <: AbstractFEStateMap
   end
 end
 
-function RepeatingAffineFEStateMap(
-    nblocks::Int,a::Function,l::Vector{<:Function},
-    U0,V0,V_φ,U_reg,φh,dΩ...;kwargs...)
-  biform = IntegrandWithMeasure(a,dΩ)
-  liforms = map(li -> IntegrandWithMeasure(li,dΩ),l)
-  RepeatingAffineFEStateMap(nblocks::Int,biform,liforms,U0,V0,V_φ,U_reg,φh,dΩ...;kwargs...)
-end
-
 const MultiFieldSpaceTypes = Union{<:MultiField.MultiFieldFESpace,<:GridapDistributed.DistributedMultiFieldFESpace}
 
 function repeat_spaces(nblocks::Integer,U0::FESpace,V0::FESpace)
@@ -840,7 +740,6 @@ function repeated_allocate_in_domain(nblocks::Integer,M::AbstractBlockMatrix)
 end
 
 get_state(m::RepeatingAffineFEStateMap) = FEFunction(get_trial_space(m),m.fwd_caches[4])
-get_measure(m::RepeatingAffineFEStateMap) = m.biform.dΩ
 get_spaces(m::RepeatingAffineFEStateMap) = m.spaces
 get_assemblers(m::RepeatingAffineFEStateMap) = (m.fwd_caches[8],m.plb_caches[2],m.adj_caches[4])
 
@@ -1046,7 +945,6 @@ function Fields.evaluate!(pcf::PDEConstrainedFunctionals,φh)
 
   U_reg = get_deriv_space(pcf.state_map)
   deriv_assem = get_deriv_assembler(pcf.state_map)
-  dΩ = get_measure(pcf.state_map)
 
   ## Foward problem
   u, u_pullback = rrule(pcf.state_map,φh)
@@ -1062,13 +960,6 @@ function Fields.evaluate!(pcf::PDEConstrainedFunctionals,φh)
     return j_val
   end
   function ∇!(F::StateParamIntegrandWithMeasure,dF,dF_analytic::Function)
-    # Analytic shape derivative
-    j_val = F(uh,φh)
-    _dF(q) = dF_analytic(q,uh,φh,dΩ...)
-    assemble_vector!(_dF,dF,deriv_assem,U_reg)
-    return j_val
-  end
-  function ∇!(F::StateParamIntegrandWithMeasure,dF,dF_analytic::AbstractIntegrandWithMeasure)
     # Analytic shape derivative
     j_val = F(uh,φh)
     _dF(q) = dF_analytic(q,uh,φh)
@@ -1089,10 +980,6 @@ end
 
 # IO
 
-function Base.show(io::IO,object::IntegrandWithMeasure)
-  print(io,"$(nameof(typeof(object)))")
-end
-
 function Base.show(io::IO,object::StateParamIntegrandWithMeasure)
   print(io,"$(nameof(typeof(object)))")
 end
diff --git a/src/GridapTopOpt.jl b/src/GridapTopOpt.jl
index 19eaa8d..ba410c0 100644
--- a/src/GridapTopOpt.jl
+++ b/src/GridapTopOpt.jl
@@ -23,7 +23,8 @@ using GridapDistributed
 using GridapDistributed: DistributedDiscreteModel, DistributedTriangulation,
   DistributedFESpace, DistributedDomainContribution, to_parray_of_arrays,
   allocate_in_domain, DistributedCellField, DistributedMultiFieldCellField,
-  DistributedMultiFieldFEBasis, BlockPMatrix, BlockPVector, change_ghost
+  DistributedMultiFieldFEBasis, BlockPMatrix, BlockPVector, change_ghost,
+  gradient, jacobian
 
 using PartitionedArrays
 using PartitionedArrays: getany, tuple_of_arrays, matching_ghost_indices
diff --git a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl
index c50c372..d8710f6 100644
--- a/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl
+++ b/src/LevelSetEvolution/UnfittedEvolution/CutFEMEvolve.jl
@@ -6,62 +6,79 @@ Burman et al. (2017). DOI: `10.1016/j.cma.2017.09.005`.
 
 # Parameters
 - `ode_solver::ODESolver`: ODE solver 
-- `model::A`: FE model
+- `Ωs::B`: `EmbeddedCollection` holding updatable triangulation and measures from GridapEmbedded
+- `dΩ_bg::C`: Measure for integration
 - `space::B`: Level-set FE space
-- `dΩ::C`: Measure for integration
 - `assembler::Assembler`: FE assembler
-- `params::D`: Tuple of stabilisation parameter `Γg`, mesh size `h`, 
-  max steps `max_steps`, and FE space `order`
+- `params::D`: Tuple of stabilisation parameter `Γg`, mesh size `h`, and
+  max steps `max_steps`
 - `cache`: Cache for evolver, initially `nothing`.
 
 # Note
-The stepsize `dt = 0.1` in `RungeKutta` is a place-holder and is updated using
-the `γ` passed to `solve!`.
+- The stepsize `dt = 0.1` in `RungeKutta` is a place-holder and is updated using
+  the `γ` passed to `solve!`.
+- We expect the EmbeddedCollection `Ωs` to contain `:F_act` and `:dF_act`. If 
+  this is not available we add it to the recipe list in `Ωs` and a warning will appear. 
 """
-mutable struct CutFEMEvolve{A,B,C,D} <: Evolver
+mutable struct CutFEMEvolve{A,B,C} <: Evolver
   ode_solver::ODESolver
-  model::A
+  Ωs::EmbeddedCollection
+  dΩ_bg::A
   space::B
-  dΩ::C
   assembler::Assembler
-  params::D
+  params::C
   cache
-  function CutFEMEvolve(model::A,V_φ::B,dΩ::C,h::Real;
+  function CutFEMEvolve(V_φ::B,Ωs::EmbeddedCollection,dΩ_bg::A,h::Real;
       max_steps=10,
       Γg = 0.1,
       ode_ls = LUSolver(),
       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,C}
-    params = (;Γg,h,max_steps,order=get_order(V_φ))
-    new{A,B,C,typeof(params)}(ode_solver,model,V_φ,dΩ,assembler,params,nothing)
+      assembler=SparseMatrixAssembler(V_φ,V_φ)) where {A,B}
+    if !(:F_act ∈ keys(Ωs.objects))
+      @warn "Expected triangulation ':F_act' not found in the 
+      EmbeddedCollection. This and the corresponding measure ':dF_act' have been
+      added to the recipe list. 
+      
+      Ensure that you are not using ':F_act' under a different
+      name to avoid additional computation for cutting."
+      function F_act_recipe(cutgeo)
+        Ω_act = Triangulation(cutgeo,ACTIVE)
+        F_act = SkeletonTriangulation(Ω_act)
+        (; 
+          :F_act => F_act,
+          :dF_act => Measure(F_act,2get_order(V_φ))
+        )
+      end
+      add_recipe!(Ωs,F_act_recipe)
+    end
+    params = (;Γg,h,max_steps)
+    new{A,B,typeof(params)}(ode_solver,Ωs,dΩ_bg,V_φ,assembler,params,nothing)
   end
 end
 
 get_ode_solver(s::CutFEMEvolve) = s.ode_solver
 get_assembler(s::CutFEMEvolve) = s.assembler
 get_space(s::CutFEMEvolve) = s.space
-get_model(s::CutFEMEvolve) = s.model
-get_measure(s::CutFEMEvolve) = s.dΩ
+get_embedded_collection(s::CutFEMEvolve) = s.Ωs
+get_measure(s::CutFEMEvolve) = s.dΩ_bg
 get_params(s::CutFEMEvolve) = s.params
 get_cache(s::CutFEMEvolve) = s.cache
 
 function get_transient_operator(φh,velh,s::CutFEMEvolve)
-  model, V_φ, dΩ, assembler, params = s.model, s.space, s.dΩ, s.assembler, s.params
-  Γg, h, order = params.Γg, params.h, params.order
-
-  geo = DiscreteGeometry(φh,model)
-  cutgeo = cut(model,geo)
-  Ω_act = Triangulation(cutgeo,ACTIVE)
-  F_act = SkeletonTriangulation(Ω_act)
-  dF_act = Measure(F_act,2*order)
-  n = get_normal_vector(F_act)
+  Ωs, V_φ, dΩ_bg, assembler, params = s.Ωs, s.space, s.dΩ_bg, s.assembler, s.params
+  Γg, h = params.Γg, params.h
   ϵ = 1e-20
 
+  update_collection!(Ωs,φh)
+  F_act = Ωs.F_act
+  dF_act = Ωs.dF_act
+  n = get_normal_vector(F_act)
+
   β = velh*∇(φh)/(ϵ + norm ∘ ∇(φh))
-  stiffness(t,u,v) = ∫((β ⋅ ∇(u)) * v)dΩ + ∫(Γg*h^2*jump(∇(u) ⋅ n)*jump(∇(v) ⋅ n))dF_act
-  mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ
-  forcing(t,v) = ∫(0v)dΩ
+  stiffness(t,u,v) = ∫((β ⋅ ∇(u)) * v)dΩ_bg + ∫(Γg*h^2*jump(∇(u) ⋅ n)*jump(∇(v) ⋅ n))dF_act
+  mass(t, ∂ₜu, v) = ∫(∂ₜu * v)dΩ_bg
+  forcing(t,v) = ∫(0v)dΩ_bg
   Ut_φ = TransientTrialFESpace(V_φ)
   ode_op = TransientLinearFEOperator((stiffness,mass),forcing,Ut_φ,V_φ;
     constant_forms=(false,true),assembler)
diff --git a/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl b/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl
index fb1e4ab..4bda0a3 100644
--- a/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl
+++ b/src/LevelSetEvolution/UnfittedEvolution/StabilisedReinit.jl
@@ -11,32 +11,49 @@ the artifical viscosity term with an interior jump penalty term.
 # Parameters
 - `nls::NonlinearSolver`: Nonlinear solver for solving the reinitialisation equation
 - `stabilisation_method::A`: A `StabilisationMethod` method for stabilising the problem 
-- `model::B`: FE model
+- `Ωs::B`: `EmbeddedCollection` holding updatable triangulation and measures from GridapEmbedded
+- `dΩ_bg::D`: Background mesh measure
 - `space::C`: FE space for level-set function
-- `dΩ::D`: Background mesh measure
 - `assembler::Assembler`: Assembler for LS FE space 
-- `params::E`: Tuple of Nitsche parameter `γd`, mesh size `h`, 
-  and FE space `order`
+- `params::E`: Tuple of Nitsche parameter `γd` and mesh size `h`
 - `cache`: Cache for evolver, initially `nothing`.
 
+# Note
+- We expect the EmbeddedCollection `Ωs` to contain `:dΓ`. If this is not 
+  available we add it to the recipe list in `Ωs` and a warning will appear. 
 """
-mutable struct StabilisedReinit{A,B,C,D,E} <: Reinitialiser
+mutable struct StabilisedReinit{A,B,C,D} <: Reinitialiser
   nls::NonlinearSolver
   stabilisation_method::A
-  model::B
+  Ωs::EmbeddedCollection
+  dΩ_bg::B
   space::C
-  dΩ::D
   assembler::Assembler
-  params::E
+  params::D
   cache
-  function StabilisedReinit(model::B,V_φ::C,dΩ::D,h::Real;
+  function StabilisedReinit(V_φ::C,Ωs::EmbeddedCollection,dΩ_bg::B,h::Real;
       γd = 20.0,
       nls = NewtonSolver(LUSolver();maxiter=50,rtol=1.e-14,verbose=true),
       assembler=SparseMatrixAssembler(V_φ,V_φ),
-      stabilisation_method::A = InteriorPenalty(V_φ)) where {A,B,C,D}
-    params = (;γd,h,order=get_order(V_φ))
-    new{A,B,C,D,typeof(params)}(nls,stabilisation_method,model,
-      V_φ,dΩ,assembler,params,nothing)
+      stabilisation_method::A = InteriorPenalty(V_φ)) where {A,B,C}
+    
+    if !(:dΓ ∈ keys(Ωs.objects))
+      @warn "Expected measure ':dΓ' not found in the 
+      EmbeddedCollection. This has been added to the recipe list. 
+      
+      Ensure that you are not using ':dΓ' under a different
+      name to avoid additional computation for cutting."
+      function dΓ_recipe(cutgeo)
+        Γ = EmbeddedBoundary(cutgeo)
+        (; 
+          :dΓ => Measure(Γ,2get_order(V_φ))
+        )
+      end
+      add_recipe!(Ωs,dΓ_recipe)
+    end
+    params = (;γd,h)
+    new{A,B,C,typeof(params)}(nls,stabilisation_method,Ωs,
+    dΩ_bg,V_φ,assembler,params,nothing)
   end
 end
 
@@ -44,8 +61,8 @@ get_nls(s::StabilisedReinit) = s.nls
 get_stabilisation_method(s::StabilisedReinit) = s.stabilisation_method
 get_assembler(s::StabilisedReinit) = s.assembler
 get_space(s::StabilisedReinit) = s.space
-get_model(s::StabilisedReinit) = s.model
-get_measure(s::StabilisedReinit) = s.dΩ
+get_embedded_collection(s::StabilisedReinit) = s.Ωs
+get_measure(s::StabilisedReinit) = s.dΩ_bg
 get_params(s::StabilisedReinit) = s.params
 get_cache(s::StabilisedReinit) = s.cache
 
@@ -78,23 +95,20 @@ struct ArtificialViscosity <: StabilisationMethod
 end
 
 function get_residual_and_jacobian(s::StabilisedReinit{ArtificialViscosity},φh)
-  model, dΩ, = s.model, s.dΩ
-  γd, h, order = s.params
+  Ωs, dΩ_bg, = s.Ωs, s.dΩ_bg
+  γd, h = s.params
   ca = s.stabilisation_method.stabilisation_coefficent
   ϵ = 1e-20
 
-  geo = DiscreteGeometry(φh,model)
-  cutgeo = cut(model,geo)
-  Γ = EmbeddedBoundary(cutgeo)
-  dΓ = Measure(Γ,2*order)
+  update_collection!(Ωs,φh)
+  dΓ = Ωs.dΓ
 
   W(u,∇u) = sign(u) * ∇u / (ϵ + norm(∇u))
   V(w) = ca*h*(sqrt ∘ ( w ⋅ w ))
-  a(w,u,v) = ∫(v*(W ∘ (w,∇(w))) ⋅ ∇(u) + V(W ∘ (w,∇(w)))*∇(u) ⋅ ∇(v))dΩ + ∫((γd/h)*v*u)dΓ
-  b(w,v) = ∫((sign ∘ w)*v)dΩ
+  a(w,u,v) = ∫(v*(W ∘ (w,∇(w))) ⋅ ∇(u) + V(W ∘ (w,∇(w)))*∇(u) ⋅ ∇(v))dΩ_bg + ∫((γd/h)*v*u)dΓ
+  b(w,v) = ∫((sign ∘ w)*v)dΩ_bg
   res(u,v) = a(u,u,v) - b(u,v)
   jac(u,du,v) = a(u,du,v)
-
   return res,jac
 end
 
@@ -108,20 +122,18 @@ function InteriorPenalty(V_φ::FESpace)
 end
 
 function get_residual_and_jacobian(s::StabilisedReinit{InteriorPenalty},φh)
-  model, dΩ, = s.model, s.dΩ
-  γd, h, order = s.params
+  Ωs, dΩ_bg, = s.Ωs, s.dΩ_bg
+  γd, h = s.params
   ϵ = 1e-20
   dΛ = s.stabilisation_method.dΛ
 
-  geo = DiscreteGeometry(φh,model)
-  cutgeo = cut(model,geo)
-  Γ = EmbeddedBoundary(cutgeo)
-  dΓ = Measure(Γ,2*order)
+  update_collection!(Ωs,φh)
+  dΓ = Ωs.dΓ
 
   W(u,∇u) = sign(u) * ∇u / (ϵ + norm(∇u))
   V(w) = ca*h*(sqrt ∘ ( w ⋅ w ))
-  a(w,u,v) = ∫(v*(W ∘ (w,∇(w))) ⋅ ∇(u))dΩ + ∫(h^2*jump(∇(u)) ⋅ jump(∇(v)))dΛ + ∫((γd/h)*v*u)dΓ
-  b(w,v) = ∫((sign ∘ w)*v)dΩ
+  a(w,u,v) = ∫(v*(W ∘ (w,∇(w))) ⋅ ∇(u))dΩ_bg + ∫(h^2*jump(∇(u)) ⋅ jump(∇(v)))dΛ + ∫((γd/h)*v*u)dΓ
+  b(w,v) = ∫((sign ∘ w)*v)dΩ_bg
   res(u,v) = a(u,u,v) - b(u,v)
   jac(u,du,v) = a(u,du,v)
 
diff --git a/src/Optimisers/HilbertianProjection.jl b/src/Optimisers/HilbertianProjection.jl
index 3d436f2..9b57612 100644
--- a/src/Optimisers/HilbertianProjection.jl
+++ b/src/Optimisers/HilbertianProjection.jl
@@ -352,7 +352,7 @@ function Base.iterate(m::HilbertianProjection,state)
   U_reg = get_deriv_space(m.problem.state_map)
   V_φ   = get_aux_space(m.problem.state_map)
   interpolate!(FEFunction(U_reg,θ),vel,V_φ)
-  J, C, dJ, dC, γ = _linesearch!(m,state)
+  J, C, dJ, dC = _linesearch!(m,state,γ)
 
   ## Hilbertian extension-regularisation
   project!(m.vel_ext,dJ)
@@ -367,8 +367,8 @@ function Base.iterate(m::HilbertianProjection,state)
   return vars, state
 end
 
-function _linesearch!(m::HilbertianProjection{WithAutoDiff},state)
-  it, J, C, θ, dJ, dC, uh, φh, vel, φ_tmp, γ, os_it = state
+function _linesearch!(m::HilbertianProjection{WithAutoDiff},state,γ)
+  it, J, C, θ, dJ, dC, uh, φh, vel, φ_tmp, _, os_it = state
 
   params = m.params; history = m.history
   ls_enabled = params.ls_enabled; reinit_mod = params.reinit_mod
@@ -409,8 +409,8 @@ function _linesearch!(m::HilbertianProjection{WithAutoDiff},state)
   return J, C, dJ, dC, γ
 end
 
-function _linesearch!(m::HilbertianProjection{NoAutoDiff},state)
-  it, J, C, θ, dJ, dC, uh, φh, vel, φ_tmp, γ, os_it = state
+function _linesearch!(m::HilbertianProjection{NoAutoDiff},state,γ)
+  it, J, C, θ, dJ, dC, uh, φh, vel, φ_tmp, _, os_it = state
 
   params = m.params; history = m.history
   ls_enabled = params.ls_enabled; reinit_mod = params.reinit_mod
diff --git a/test/mpi/InverseHomogenisationALMTests.jl b/test/mpi/InverseHomogenisationALMTests.jl
index dcac064..f124c22 100644
--- a/test/mpi/InverseHomogenisationALMTests.jl
+++ b/test/mpi/InverseHomogenisationALMTests.jl
@@ -61,22 +61,22 @@ function main(distribute,mesh_partition;AD)
         TensorValue(0.,0.,0.,1.),     # ϵᵢⱼ⁽²²⁾≡ϵᵢⱼ⁽²⁾
         TensorValue(0.,1/2,1/2,0.))   # ϵᵢⱼ⁽¹²⁾≡ϵᵢⱼ⁽³⁾
 
-  a(u,v,φ,dΩ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v) )dΩ
-  l = [(v,φ,dΩ) -> ∫(-(I ∘ φ)* C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3]
+  a(u,v,φ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v) )dΩ
+  l = [(v,φ) -> ∫(-(I ∘ φ)* C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3]
 
   ## Optimisation functionals
-  Cᴴ(r,s,u,φ,dΩ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ
-  dCᴴ(r,s,q,u,φ,dΩ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
-  κ(u,φ,dΩ) = -1/4*(Cᴴ(1,1,u,φ,dΩ)+Cᴴ(2,2,u,φ,dΩ)+2*Cᴴ(1,2,u,φ,dΩ))
-  dκ(q,u,φ,dΩ) = -1/4*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+2*dCᴴ(1,2,q,u,φ,dΩ))
-  Vol(u,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
-  dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  Cᴴ(r,s,u,φ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ
+  dCᴴ(r,s,q,u,φ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  κ(u,φ) = -1/4*(Cᴴ(1,1,u,φ)+Cᴴ(2,2,u,φ)+2*Cᴴ(1,2,u,φ))
+  dκ(q,u,φ) = -1/4*(dCᴴ(1,1,q,u,φ)+dCᴴ(2,2,q,u,φ)+2*dCᴴ(1,2,q,u,φ))
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh,dΩ)
+  state_map = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh)
   pcfs = if AD
     PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol])
   else
diff --git a/test/mpi/InverterHPMTests.jl b/test/mpi/InverterHPMTests.jl
index f4c614d..0450dfa 100644
--- a/test/mpi/InverterHPMTests.jl
+++ b/test/mpi/InverterHPMTests.jl
@@ -78,21 +78,21 @@ function main(distribute,mesh_partition)
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_in,dΓ_out) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out
-  l(v,φ,dΩ,dΓ_in,dΓ_out) = ∫(v⋅g)dΓ_in
+  a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out
+  l(v,φ) = ∫(v⋅g)dΓ_in
 
   ## Optimisation functionals
   e₁ = VectorValue(1,0)
-  J(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in
-  Vol(u,φ,dΩ,dΓ_in,dΓ_out) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
-  dVol(q,u,φ,dΩ,dΓ_in,dΓ_out) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
-  UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out
+  J(u,φ) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  UΓ_out(u,φ) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_in,dΓ_out)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol,UΓ_out],state_map,analytic_dC=[dVol,nothing])
 
   ## Hilbertian extension-regularisation problems
diff --git a/test/mpi/NeohookAnalyticJacALMTests.jl b/test/mpi/NeohookAnalyticJacALMTests.jl
index df8041a..bf9936c 100644
--- a/test/mpi/NeohookAnalyticJacALMTests.jl
+++ b/test/mpi/NeohookAnalyticJacALMTests.jl
@@ -85,21 +85,21 @@ function main(distribute,mesh_partition)
 
   # Cauchy stress tensor and residual
   σ(∇u) = (1.0/J(F(∇u)))*F(∇u)⋅S(∇u)⋅(F(∇u))'
-  res(u,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*((dE∘(∇(v),∇(u))) ⊙ (S∘∇(u))) )*dΩ - ∫(g⋅v)dΓ_N
-  jac_mat(u,du,v,φ,dΩ,dΓ_N) =  ∫( (I ∘ φ)*(dE∘(∇(v),∇(u))) ⊙ (dS∘(∇(du),∇(u))) )*dΩ
-  jac_geo(u,du,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*∇(v) ⊙ ( (S∘∇(u))⋅∇(du) ) )*dΩ
-  jac(u,du,v,φ,dΩ,dΓ_N) = jac_mat(u,du,v,φ,dΩ,dΓ_N) + jac_geo(u,du,v,φ,dΩ,dΓ_N)
+  res(u,v,φ) = ∫( (I ∘ φ)*((dE∘(∇(v),∇(u))) ⊙ (S∘∇(u))) )*dΩ - ∫(g⋅v)dΓ_N
+  jac_mat(u,du,v,φ) =  ∫( (I ∘ φ)*(dE∘(∇(v),∇(u))) ⊙ (dS∘(∇(du),∇(u))) )*dΩ
+  jac_geo(u,du,v,φ) = ∫( (I ∘ φ)*∇(v) ⊙ ( (S∘∇(u))⋅∇(du) ) )*dΩ
+  jac(u,du,v,φ) = jac_mat(u,du,v,φ) + jac_geo(u,du,v,φ)
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((dE∘(∇(u),∇(u))) ⊙ (S∘∇(u))))dΩ
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*((dE∘(∇(u),∇(u))) ⊙ (S∘∇(u))))dΩ
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N;jac)
+  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh;jac)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/test/mpi/NonlinearThermalComplianceALMTests.jl b/test/mpi/NonlinearThermalComplianceALMTests.jl
index 1b63a03..134a4cb 100644
--- a/test/mpi/NonlinearThermalComplianceALMTests.jl
+++ b/test/mpi/NonlinearThermalComplianceALMTests.jl
@@ -68,18 +68,18 @@ function main(distribute,mesh_partition)
   κ0 = 1
   ξ = -1
   κ(u) = κ0*(exp(ξ*u))
-  res(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N
+  res(u,v,φ) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(u))dΩ
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(u))dΩ
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/test/mpi/PZMultiFieldRepeatingStateTests.jl b/test/mpi/PZMultiFieldRepeatingStateTests.jl
index dbe6f25..c083b7f 100644
--- a/test/mpi/PZMultiFieldRepeatingStateTests.jl
+++ b/test/mpi/PZMultiFieldRepeatingStateTests.jl
@@ -67,21 +67,21 @@ function main(distribute,mesh_partition;AD)
   Eⁱ = (VectorValue(1.0,0.0,),
         VectorValue(0.0,1.0))
 
-  a((u,ϕ),(v,q),φ,dΩ) = ∫((I ∘ φ) * (1/k0*((C ⊙ ε(u)) ⊙ ε(v)) -
+  a((u,ϕ),(v,q),φ) = ∫((I ∘ φ) * (1/k0*((C ⊙ ε(u)) ⊙ ε(v)) -
                                     1/α0*((-∇(ϕ) ⋅ e) ⊙ ε(v)) +
                                     -1/α0*((e ⋅² ε(u)) ⋅ -∇(q)) +
                                     -γ0/β0*((κ ⋅ -∇(ϕ)) ⋅ -∇(q))) )dΩ;
 
-  l_ε = [((v,q),φ,dΩ) -> ∫(((I ∘ φ) * (-C ⊙ εᴹ[i] ⊙ ε(v) + k0/α0*(e ⋅² εᴹ[i]) ⋅ -∇(q))))dΩ for i = 1:3];
-  l_E = [((v,q),φ,dΩ) -> ∫((I ∘ φ) * ((Eⁱ[i] ⋅ e ⊙ ε(v) + k0/α0*(κ ⋅ Eⁱ[i]) ⋅ -∇(q))))dΩ for i = 1:2];
+  l_ε = [((v,q),φ) -> ∫(((I ∘ φ) * (-C ⊙ εᴹ[i] ⊙ ε(v) + k0/α0*(e ⋅² εᴹ[i]) ⋅ -∇(q))))dΩ for i = 1:3];
+  l_E = [((v,q),φ) -> ∫((I ∘ φ) * ((Eⁱ[i] ⋅ e ⊙ ε(v) + k0/α0*(κ ⋅ Eⁱ[i]) ⋅ -∇(q))))dΩ for i = 1:2];
   l = [l_ε; l_E]
 
-  function Cᴴ(r,s,uϕ,φ,dΩ)
+  function Cᴴ(r,s,uϕ,φ)
     u_s = uϕ[2s-1]; ϕ_s = uϕ[2s]
     ∫(1/k0 * (I ∘ φ) * (((C ⊙ (1/k0*ε(u_s) + εᴹ[s])) ⊙ εᴹ[r]) - ((-1/α0*∇(ϕ_s) ⋅ e) ⊙ εᴹ[r])))dΩ;
   end
 
-  function DCᴴ(r,s,q,uϕ,φ,dΩ)
+  function DCᴴ(r,s,q,uϕ,φ)
     u_r = uϕ[2r-1]; ϕ_r = uϕ[2r]
     u_s = uϕ[2s-1]; ϕ_s = uϕ[2s]
     ∫(- 1/k0 * q * (
@@ -93,19 +93,19 @@ function main(distribute,mesh_partition;AD)
     )dΩ;
   end
 
-  Bᴴ(uϕ,φ,dΩ) = 1/4*(Cᴴ(1,1,uϕ,φ,dΩ)+Cᴴ(2,2,uϕ,φ,dΩ)+2*Cᴴ(1,2,uϕ,φ,dΩ))
-  DBᴴ(q,uϕ,φ,dΩ) = 1/4*(DCᴴ(1,1,q,uϕ,φ,dΩ)+DCᴴ(2,2,q,uϕ,φ,dΩ)+2*DCᴴ(1,2,q,uϕ,φ,dΩ))
+  Bᴴ(uϕ,φ) = 1/4*(Cᴴ(1,1,uϕ,φ)+Cᴴ(2,2,uϕ,φ)+2*Cᴴ(1,2,uϕ,φ))
+  DBᴴ(q,uϕ,φ) = 1/4*(DCᴴ(1,1,q,uϕ,φ)+DCᴴ(2,2,q,uϕ,φ)+2*DCᴴ(1,2,q,uϕ,φ))
 
-  J(uϕ,φ,dΩ) = -1*Bᴴ(uϕ,φ,dΩ)
-  DJ(q,uϕ,φ,dΩ) = -1*DBᴴ(q,uϕ,φ,dΩ)
-  C1(uϕ,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  DC1(q,uϕ,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(uϕ,φ) = -1*Bᴴ(uϕ,φ)
+  DJ(q,uϕ,φ) = -1*DBᴴ(q,uϕ,φ)
+  C1(uϕ,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  DC1(q,uϕ,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   stencil = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = RepeatingAffineFEStateMap(5,a,l,UP,VQ,V_φ,U_reg,φh,dΩ)
+  state_map = RepeatingAffineFEStateMap(5,a,l,UP,VQ,V_φ,U_reg,φh)
   pcfs = if AD
     PDEConstrainedFunctionals(J,[C1],state_map;analytic_dJ=DJ,analytic_dC=[DC1])
   else
diff --git a/test/mpi/ThermalComplianceALMTests.jl b/test/mpi/ThermalComplianceALMTests.jl
index 9b7e0a6..dee3d38 100644
--- a/test/mpi/ThermalComplianceALMTests.jl
+++ b/test/mpi/ThermalComplianceALMTests.jl
@@ -63,20 +63,20 @@ function main(distribute,mesh_partition;order,AD)
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
+  l(v,φ) = ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
+  dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = if AD
     PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol])
   else
@@ -157,20 +157,20 @@ function main_3d(distribute,mesh_partition,;order)
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
+  l(v,φ) = ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
+  dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/test/mpi/ThermalComplianceALMTestsPETSc.jl b/test/mpi/ThermalComplianceALMTestsPETSc.jl
index db1b6d4..b3832bd 100644
--- a/test/mpi/ThermalComplianceALMTestsPETSc.jl
+++ b/test/mpi/ThermalComplianceALMTestsPETSc.jl
@@ -63,14 +63,14 @@ function main(distribute,mesh_partition;order,AD)
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
+  l(v,φ) = ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
+  dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
@@ -81,7 +81,7 @@ function main(distribute,mesh_partition;order,AD)
   solver = PETScLinearSolver()
 
   state_map = AffineFEStateMap(
-    a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N;
+    a,l,U,V,V_φ,U_reg,φh;
     assem_U = SparseMatrixAssembler(Tm,Tv,U,V),
     assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U),
     assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg),
@@ -171,14 +171,14 @@ function main_3d(distribute,mesh_partition,;order)
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
+  l(v,φ) = ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
+  dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
@@ -189,7 +189,7 @@ function main_3d(distribute,mesh_partition,;order)
   solver = PETScLinearSolver()
 
   state_map = AffineFEStateMap(
-    a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N;
+    a,l,U,V,V_φ,U_reg,φh;
     assem_U = SparseMatrixAssembler(Tm,Tv,U,V),
     assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U),
     assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg),
diff --git a/test/seq/InverseHomogenisationALMTests.jl b/test/seq/InverseHomogenisationALMTests.jl
index b0b79d4..3d1a264 100644
--- a/test/seq/InverseHomogenisationALMTests.jl
+++ b/test/seq/InverseHomogenisationALMTests.jl
@@ -59,22 +59,22 @@ function main(;AD)
         TensorValue(0.,0.,0.,1.),     # ϵᵢⱼ⁽²²⁾≡ϵᵢⱼ⁽²⁾
         TensorValue(0.,1/2,1/2,0.))   # ϵᵢⱼ⁽¹²⁾≡ϵᵢⱼ⁽³⁾
 
-  a(u,v,φ,dΩ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v) )dΩ
-  l = [(v,φ,dΩ) -> ∫(-(I ∘ φ)* C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3]
+  a(u,v,φ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v) )dΩ
+  l = [(v,φ) -> ∫(-(I ∘ φ)* C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3]
 
   ## Optimisation functionals
-  Cᴴ(r,s,u,φ,dΩ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ
-  dCᴴ(r,s,q,u,φ,dΩ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
-  κ(u,φ,dΩ) = -1/4*(Cᴴ(1,1,u,φ,dΩ)+Cᴴ(2,2,u,φ,dΩ)+2*Cᴴ(1,2,u,φ,dΩ))
-  dκ(q,u,φ,dΩ) = -1/4*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+2*dCᴴ(1,2,q,u,φ,dΩ))
-  Vol(u,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
-  dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  Cᴴ(r,s,u,φ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ
+  dCᴴ(r,s,q,u,φ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  κ(u,φ) = -1/4*(Cᴴ(1,1,u,φ)+Cᴴ(2,2,u,φ)+2*Cᴴ(1,2,u,φ))
+  dκ(q,u,φ) = -1/4*(dCᴴ(1,1,q,u,φ)+dCᴴ(2,2,q,u,φ)+2*dCᴴ(1,2,q,u,φ))
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh,dΩ)
+  state_map = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh)
   pcfs = if AD
     PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol])
   else
diff --git a/test/seq/InverterHPMTests.jl b/test/seq/InverterHPMTests.jl
index 06fffac..1511cb8 100644
--- a/test/seq/InverterHPMTests.jl
+++ b/test/seq/InverterHPMTests.jl
@@ -76,21 +76,21 @@ function main()
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_in,dΓ_out) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out
-  l(v,φ,dΩ,dΓ_in,dΓ_out) = ∫(v⋅g)dΓ_in
+  a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out
+  l(v,φ) = ∫(v⋅g)dΓ_in
 
   ## Optimisation functionals
   e₁ = VectorValue(1,0)
-  J(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in
-  Vol(u,φ,dΩ,dΓ_in,dΓ_out) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
-  dVol(q,u,φ,dΩ,dΓ_in,dΓ_out) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
-  UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out
+  J(u,φ) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  UΓ_out(u,φ) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_in,dΓ_out)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol,UΓ_out],state_map,analytic_dC=[dVol,nothing])
 
   ## Hilbertian extension-regularisation problems
diff --git a/test/seq/NeohookAnalyticJacALMTests.jl b/test/seq/NeohookAnalyticJacALMTests.jl
index 8ef5d70..e4b8bc0 100644
--- a/test/seq/NeohookAnalyticJacALMTests.jl
+++ b/test/seq/NeohookAnalyticJacALMTests.jl
@@ -83,21 +83,21 @@ function main()
 
   # Cauchy stress tensor and residual
   σ(∇u) = (1.0/J(F(∇u)))*F(∇u)⋅S(∇u)⋅(F(∇u))'
-  res(u,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*((dE∘(∇(v),∇(u))) ⊙ (S∘∇(u))) )*dΩ - ∫(g⋅v)dΓ_N
-  jac_mat(u,du,v,φ,dΩ,dΓ_N) =  ∫( (I ∘ φ)*(dE∘(∇(v),∇(u))) ⊙ (dS∘(∇(du),∇(u))) )*dΩ
-  jac_geo(u,du,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*∇(v) ⊙ ( (S∘∇(u))⋅∇(du) ) )*dΩ
-  jac(u,du,v,φ,dΩ,dΓ_N) = jac_mat(u,du,v,φ,dΩ,dΓ_N) + jac_geo(u,du,v,φ,dΩ,dΓ_N)
+  res(u,v,φ) = ∫( (I ∘ φ)*((dE∘(∇(v),∇(u))) ⊙ (S∘∇(u))) )*dΩ - ∫(g⋅v)dΓ_N
+  jac_mat(u,du,v,φ) =  ∫( (I ∘ φ)*(dE∘(∇(v),∇(u))) ⊙ (dS∘(∇(du),∇(u))) )*dΩ
+  jac_geo(u,du,v,φ) = ∫( (I ∘ φ)*∇(v) ⊙ ( (S∘∇(u))⋅∇(du) ) )*dΩ
+  jac(u,du,v,φ) = jac_mat(u,du,v,φ) + jac_geo(u,du,v,φ)
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((dE∘(∇(u),∇(u))) ⊙ (S∘∇(u))))dΩ
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*((dE∘(∇(u),∇(u))) ⊙ (S∘∇(u))))dΩ
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N;jac)
+  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh;jac)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/test/seq/NonlinearThermalComplianceALMTests.jl b/test/seq/NonlinearThermalComplianceALMTests.jl
index c480870..c50f61c 100644
--- a/test/seq/NonlinearThermalComplianceALMTests.jl
+++ b/test/seq/NonlinearThermalComplianceALMTests.jl
@@ -66,18 +66,18 @@ function main()
   κ0 = 1
   ξ = -1
   κ(u) = κ0*(exp(ξ*u))
-  res(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N
+  res(u,v,φ) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(u))dΩ
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(u))dΩ
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/test/seq/PZMultiFieldRepeatingStateTests.jl b/test/seq/PZMultiFieldRepeatingStateTests.jl
index 5864d9f..8ed6c07 100644
--- a/test/seq/PZMultiFieldRepeatingStateTests.jl
+++ b/test/seq/PZMultiFieldRepeatingStateTests.jl
@@ -72,21 +72,21 @@ function main(;AD,use_mfs=false)
   Eⁱ = (VectorValue(1.0,0.0,),
         VectorValue(0.0,1.0))
 
-  a((u,ϕ),(v,q),φ,dΩ) = ∫((I ∘ φ) * (1/k0*((C ⊙ ε(u)) ⊙ ε(v)) -
+  a((u,ϕ),(v,q),φ) = ∫((I ∘ φ) * (1/k0*((C ⊙ ε(u)) ⊙ ε(v)) -
                                     1/α0*((-∇(ϕ) ⋅ e) ⊙ ε(v)) +
                                     -1/α0*((e ⋅² ε(u)) ⋅ -∇(q)) +
                                     -γ0/β0*((κ ⋅ -∇(ϕ)) ⋅ -∇(q))) )dΩ;
 
-  l_ε = [((v,q),φ,dΩ) -> ∫(((I ∘ φ) * (-C ⊙ εᴹ[i] ⊙ ε(v) + k0/α0*(e ⋅² εᴹ[i]) ⋅ -∇(q))))dΩ for i = 1:3];
-  l_E = [((v,q),φ,dΩ) -> ∫((I ∘ φ) * ((Eⁱ[i] ⋅ e ⊙ ε(v) + k0/α0*(κ ⋅ Eⁱ[i]) ⋅ -∇(q))))dΩ for i = 1:2];
+  l_ε = [((v,q),φ) -> ∫(((I ∘ φ) * (-C ⊙ εᴹ[i] ⊙ ε(v) + k0/α0*(e ⋅² εᴹ[i]) ⋅ -∇(q))))dΩ for i = 1:3];
+  l_E = [((v,q),φ) -> ∫((I ∘ φ) * ((Eⁱ[i] ⋅ e ⊙ ε(v) + k0/α0*(κ ⋅ Eⁱ[i]) ⋅ -∇(q))))dΩ for i = 1:2];
   l = [l_ε; l_E]
 
-  function Cᴴ(r,s,uϕ,φ,dΩ)
+  function Cᴴ(r,s,uϕ,φ)
     u_s = uϕ[2s-1]; ϕ_s = uϕ[2s]
     ∫(1/k0 * (I ∘ φ) * (((C ⊙ (1/k0*ε(u_s) + εᴹ[s])) ⊙ εᴹ[r]) - ((-1/α0*∇(ϕ_s) ⋅ e) ⊙ εᴹ[r])))dΩ;
   end
 
-  function DCᴴ(r,s,q,uϕ,φ,dΩ)
+  function DCᴴ(r,s,q,uϕ,φ)
     u_r = uϕ[2r-1]; ϕ_r = uϕ[2r]
     u_s = uϕ[2s-1]; ϕ_s = uϕ[2s]
     ∫(- 1/k0 * q * (
@@ -98,19 +98,19 @@ function main(;AD,use_mfs=false)
     )dΩ;
   end
 
-  Bᴴ(uϕ,φ,dΩ) = 1/4*(Cᴴ(1,1,uϕ,φ,dΩ)+Cᴴ(2,2,uϕ,φ,dΩ)+2*Cᴴ(1,2,uϕ,φ,dΩ))
-  DBᴴ(q,uϕ,φ,dΩ) = 1/4*(DCᴴ(1,1,q,uϕ,φ,dΩ)+DCᴴ(2,2,q,uϕ,φ,dΩ)+2*DCᴴ(1,2,q,uϕ,φ,dΩ))
+  Bᴴ(uϕ,φ) = 1/4*(Cᴴ(1,1,uϕ,φ)+Cᴴ(2,2,uϕ,φ)+2*Cᴴ(1,2,uϕ,φ))
+  DBᴴ(q,uϕ,φ) = 1/4*(DCᴴ(1,1,q,uϕ,φ)+DCᴴ(2,2,q,uϕ,φ)+2*DCᴴ(1,2,q,uϕ,φ))
 
-  J(uϕ,φ,dΩ) = -1*Bᴴ(uϕ,φ,dΩ)
-  DJ(q,uϕ,φ,dΩ) = -1*DBᴴ(q,uϕ,φ,dΩ)
-  C1(uϕ,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  DC1(q,uϕ,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(uϕ,φ) = -1*Bᴴ(uϕ,φ)
+  DJ(q,uϕ,φ) = -1*DBᴴ(q,uϕ,φ)
+  C1(uϕ,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  DC1(q,uϕ,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   stencil = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = RepeatingAffineFEStateMap(5,a,l,UP,VQ,V_φ,U_reg,φh,dΩ)
+  state_map = RepeatingAffineFEStateMap(5,a,l,UP,VQ,V_φ,U_reg,φh)
   pcfs = if AD
     PDEConstrainedFunctionals(J,[C1],state_map;analytic_dJ=DJ,analytic_dC=[DC1])
   else
diff --git a/test/seq/ThermalComplianceALMTests.jl b/test/seq/ThermalComplianceALMTests.jl
index fcf3704..33591d4 100644
--- a/test/seq/ThermalComplianceALMTests.jl
+++ b/test/seq/ThermalComplianceALMTests.jl
@@ -61,20 +61,20 @@ function main(;order,AD)
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
+  l(v,φ) = ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
+  dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = if AD
     PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol])
   else
@@ -154,20 +154,20 @@ function main_3d(;order)
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
+  l(v,φ) = ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
+  dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol])
 
   ## Hilbertian extension-regularisation problems
diff --git a/test/seq/ThermalComplianceHPMTests.jl b/test/seq/ThermalComplianceHPMTests.jl
index 87fce14..2ebe2a5 100644
--- a/test/seq/ThermalComplianceHPMTests.jl
+++ b/test/seq/ThermalComplianceHPMTests.jl
@@ -61,20 +61,20 @@ function main(;order,AD_case)
   interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ))
   I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ
 
-  a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
-  l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N
+  a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ
+  l(v,φ) = ∫(v)dΓ_N
 
   ## Optimisation functionals
-  J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
-  dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
-  Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
-  dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
+  J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
+  dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ;
+  Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ;
+  dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
 
   ## Finite difference solver and level set function
   ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps)
 
   ## Setup solver and FE operators
-  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
+  state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh)
   pcfs = if AD_case == :no_ad
     PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol])
   elseif AD_case == :with_ad
diff --git a/test/seq/UnfittedEvolutionTests/ArtificialViscosityStabilisedReinitTest.jl b/test/seq/UnfittedEvolutionTests/ArtificialViscosityStabilisedReinitTest.jl
index 34bcb62..12bf6c8 100644
--- a/test/seq/UnfittedEvolutionTests/ArtificialViscosityStabilisedReinitTest.jl
+++ b/test/seq/UnfittedEvolutionTests/ArtificialViscosityStabilisedReinitTest.jl
@@ -22,8 +22,16 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-(x[1]-0.5)^2-(x[2]-0.5)^2+0.25^2,V_φ)
 φh0 = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ)
 
-ls_evo = CutFEMEvolve(model,V_φ,dΩ,h)
-ls_reinit = StabilisedReinit(model,V_φ,dΩ,h;stabilisation_method=ArtificialViscosity(1.5h))
+Ωs = EmbeddedCollection(model,φh) do cutgeo
+  Γ = EmbeddedBoundary(cutgeo)
+  (; 
+    :Γ => Γ,
+    :dΓ => Measure(Γ,2*order)
+  )
+end
+
+ls_evo = CutFEMEvolve(V_φ,Ωs,dΩ,h)
+ls_reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=ArtificialViscosity(1.5h))
 evo = UnfittedFEEvolution(ls_evo,ls_reinit)
 reinit!(evo,φh);
 
diff --git a/test/seq/UnfittedEvolutionTests/CutFEMEvolveTest.jl b/test/seq/UnfittedEvolutionTests/CutFEMEvolveTest.jl
index a73ee3d..1fe19f0 100644
--- a/test/seq/UnfittedEvolutionTests/CutFEMEvolveTest.jl
+++ b/test/seq/UnfittedEvolutionTests/CutFEMEvolveTest.jl
@@ -3,6 +3,7 @@ using Test
 
 using GridapTopOpt
 using Gridap, Gridap.Geometry, Gridap.Adaptivity
+using GridapEmbedded
 
 order = 1
 n = 201
@@ -20,8 +21,16 @@ V_φ = TestFESpace(model,reffe_scalar)
 
 φh = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ)
 
-ls_evo = CutFEMEvolve(model,V_φ,dΩ,h)
-ls_reinit = StabilisedReinit(model,V_φ,dΩ,h)
+Ωs = EmbeddedCollection(model,φh) do cutgeo
+  Γ = EmbeddedBoundary(cutgeo)
+  (; 
+    :Γ  => Γ,
+    :dΓ => Measure(Γ,2*order),
+  )
+end
+
+ls_evo = CutFEMEvolve(V_φ,Ωs,dΩ,h)
+ls_reinit = StabilisedReinit(V_φ,Ωs,dΩ,h)
 evo = UnfittedFEEvolution(ls_evo,ls_reinit)
 
 φ0 = copy(get_free_dof_values(φh))
diff --git a/test/seq/UnfittedEvolutionTests/InteriorPenaltyStabilisedReinitTest.jl b/test/seq/UnfittedEvolutionTests/InteriorPenaltyStabilisedReinitTest.jl
index c914fc3..f35a7b9 100644
--- a/test/seq/UnfittedEvolutionTests/InteriorPenaltyStabilisedReinitTest.jl
+++ b/test/seq/UnfittedEvolutionTests/InteriorPenaltyStabilisedReinitTest.jl
@@ -22,8 +22,16 @@ V_φ = TestFESpace(model,reffe_scalar)
 φh = interpolate(x->-(x[1]-0.5)^2-(x[2]-0.5)^2+0.25^2,V_φ)
 φh0 = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ)
 
-ls_evo = CutFEMEvolve(model,V_φ,dΩ,h)
-ls_reinit = StabilisedReinit(model,V_φ,dΩ,h;stabilisation_method=InteriorPenalty(V_φ))
+Ωs = EmbeddedCollection(model,φh) do cutgeo
+  Γ = EmbeddedBoundary(cutgeo)
+  (; 
+    :Γ => Γ,
+    :dΓ => Measure(Γ,2*order)
+  )
+end
+
+ls_evo = CutFEMEvolve(V_φ,Ωs,dΩ,h)
+ls_reinit = StabilisedReinit(V_φ,Ωs,dΩ,h;stabilisation_method=InteriorPenalty(V_φ))
 evo = UnfittedFEEvolution(ls_evo,ls_reinit)
 reinit!(evo,φh);