diff --git a/Manifest.toml b/Manifest.toml
index 59c6cd3f..e5c90b75 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -148,7 +148,7 @@ version = "0.10.20"
 
 [[Gridap]]
 deps = ["AbstractTrees", "BSON", "BlockArrays", "Combinatorics", "DocStringExtensions", "FastGaussQuadrature", "FileIO", "FillArrays", "ForwardDiff", "JLD2", "JSON", "LineSearches", "LinearAlgebra", "NLsolve", "NearestNeighbors", "QuadGK", "Random", "SparseArrays", "SparseMatricesCSR", "StaticArrays", "Test", "WriteVTK"]
-git-tree-sha1 = "423ab291db7df5f86d796ab890e22532649161b6"
+git-tree-sha1 = "c58a6aade499e243bddc152f6acdc17d45886a94"
 repo-rev = "gridap_distributed"
 repo-url = "https://github.com/gridap/Gridap.jl.git"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
diff --git a/src/CellData.jl b/src/CellData.jl
index 25132e4e..a1b4794f 100644
--- a/src/CellData.jl
+++ b/src/CellData.jl
@@ -267,6 +267,53 @@ function CellData.get_cell_points(a::DistributedTriangulation)
 end
 
 function CellData.get_normal_vector(a::DistributedTriangulation)
-  DistributedCellField(map_parts(get_normal_vector,a.trians))
+  fields = map_parts(get_normal_vector,a.trians)
+  DistributedCellField(fields)
+end
+
+# Skeleton related
+
+function DistributedCellField(a::AbstractPData{<:SkeletonPair})
+  plus, minus = map_parts(s->(s.plus,s.minus),a)
+  dplus = DistributedCellField(plus)
+  dminus = DistributedCellField(minus)
+  SkeletonPair(dplus,dminus)
+end
+
+function Base.getproperty(x::DistributedCellField, sym::Symbol)
+  if sym in (:⁺,:plus)
+    DistributedCellField(map_parts(i->i.plus,x.fields))
+  elseif sym in (:⁻, :minus)
+    DistributedCellField(map_parts(i->i.minus,x.fields))
+  else
+    getfield(x, sym)
+  end
+end
+
+function Base.propertynames(x::DistributedCellField, private::Bool=false)
+  (fieldnames(typeof(x))...,:⁺,:plus,:⁻,:minus)
 end
 
+for op in (:outer,:*,:dot)
+  @eval begin
+    ($op)(a::DistributedCellField,b::SkeletonPair{<:DistributedCellField}) = Operation($op)(a,b)
+    ($op)(a::SkeletonPair{<:DistributedCellField},b::DistributedCellField) = Operation($op)(a,b)
+  end
+end
+
+function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField,b::SkeletonPair{<:DistributedCellField})
+  plus = k(a.plus,b.plus)
+  minus = k(a.minus,b.minus)
+  SkeletonPair(plus,minus)
+end
+
+function Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:DistributedCellField},b::DistributedCellField)
+  plus = k(a.plus,b.plus)
+  minus = k(a.minus,b.minus)
+  SkeletonPair(plus,minus)
+end
+
+CellData.jump(a::DistributedCellField) = DistributedCellField(map_parts(jump,a.fields))
+CellData.jump(a::SkeletonPair{<:DistributedCellField}) = a.⁺ + a.⁻
+CellData.mean(a::DistributedCellField) = DistributedCellField(map_parts(mean,a.fields))
+
diff --git a/src/FESpaces.jl b/src/FESpaces.jl
index 4b1682a6..2750544e 100644
--- a/src/FESpaces.jl
+++ b/src/FESpaces.jl
@@ -257,6 +257,11 @@ end
 
 function FESpaces.interpolate(u,f::DistributedSingleFieldFESpace)
   free_values = zero_free_values(f)
+  interpolate!(u,free_values,f)
+end
+
+function FESpaces.interpolate!(
+  u,free_values::AbstractVector,f::DistributedSingleFieldFESpace)
   map_parts(f.spaces,local_views(free_values)) do V,vec
     interpolate!(u,vec,V)
   end
diff --git a/src/Geometry.jl b/src/Geometry.jl
index 35d3ff82..54d73d41 100644
--- a/src/Geometry.jl
+++ b/src/Geometry.jl
@@ -244,8 +244,16 @@ function remove_ghost_cells(glue::FaceToFaceGlue,trian,gids)
   view(trian, findall(tcell_to_mask))
 end
 
-function remove_ghost_cells(glue::SkeletonPair,trian,gids)
-  plus = remove_ghost_cells(glue.plus,trian,gids)
-  minus = remove_ghost_cells(glue.minus,trian,gids)
+function remove_ghost_cells(glue::SkeletonPair,trian::SkeletonTriangulation,gids)
+  owner_glue = glue.plus
+  plus = remove_ghost_cells(owner_glue,trian.plus,gids)
+  minus = remove_ghost_cells(owner_glue,trian.minus,gids)
   SkeletonTriangulation(plus,minus)
 end
+
+function remove_ghost_cells(glue::SkeletonPair,trian,gids)
+  owner_glue = glue.plus
+  remove_ghost_cells(owner_glue,trian,gids)
+end
+
+
diff --git a/src/MultiField.jl b/src/MultiField.jl
index 57922728..93cd8537 100644
--- a/src/MultiField.jl
+++ b/src/MultiField.jl
@@ -97,6 +97,26 @@ function FESpaces.EvaluationFunction(
   DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values)
 end
 
+function FESpaces.interpolate(objects,fe::DistributedMultiFieldFESpace)
+  free_values = zero_free_values(fe)
+  interpolate!(objects,free_values,fe)
+end
+
+function FESpaces.interpolate!(objects,free_values::AbstractVector,fe::DistributedMultiFieldFESpace)
+  local_vals = consistent_local_views(free_values,fe.gids,true)
+  part_fe_fun = map_parts(local_vals,local_views(fe)) do x,f
+    interpolate!(objects,x,f)
+  end
+  field_fe_fun = DistributedSingleFieldFEFunction[]
+  for i in 1:num_fields(fe)
+    free_values_i = restrict_to_field(fe,free_values,i)
+    fe_space_i = fe.field_fe_space[i]
+    fe_fun_i = FEFunction(fe_space_i,free_values_i)
+    push!(field_fe_fun,fe_fun_i)
+  end
+  DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values)
+end
+
 struct DistributedMultiFieldFEBasis{A,B} <: GridapType
   field_fe_basis::A
   part_fe_basis::B
@@ -137,8 +157,8 @@ function FESpaces.get_trial_fe_basis(f::DistributedMultiFieldFESpace)
   DistributedMultiFieldFEBasis(field_fe_basis,part_mbasis)
 end
 
-function FESpaces.interpolate(u,f::DistributedMultiFieldFESpace)
-  @notimplemented
+function FESpaces.TrialFESpace(objects,a::DistributedMultiFieldFESpace)
+  TrialFESpace(a,objects)
 end
 
 function FESpaces.TrialFESpace(a::DistributedMultiFieldFESpace,objects)
@@ -294,3 +314,4 @@ function propagate_to_ghost_multifield!(
   end
 end
 
+
diff --git a/test/MultiFieldTests.jl b/test/MultiFieldTests.jl
index d37fd787..b92b0679 100644
--- a/test/MultiFieldTests.jl
+++ b/test/MultiFieldTests.jl
@@ -31,9 +31,15 @@ function main(parts)
 
   VxQ = MultiFieldFESpace([V,Q])
   UxP = MultiFieldFESpace([U,P]) # This generates again the global numbering
-  UxP = TrialFESpace(VxQ,[u,p]) # This reuses the one computed
+  UxP = TrialFESpace([u,p],VxQ) # This reuses the one computed
+
+  uh, ph = interpolate([u,p],UxP)
+  eu = u - uh
+  ep = p - ph
 
   dΩ = Measure(Ω,2*k)
+  @test sqrt(sum(∫( eu⋅eu )dΩ)) < 1.0e-9
+  @test sqrt(sum(∫( eu⋅eu )dΩ)) < 1.0e-9
 
   a((u,p),(v,q)) = ∫( ∇(v)⊙∇(u) - q*(∇⋅u) - (∇⋅v)*p )*dΩ
   l((v,q)) = ∫( v⋅f - q*g )*dΩ
diff --git a/test/PoissonTests.jl b/test/PoissonTests.jl
new file mode 100644
index 00000000..6d2ad96f
--- /dev/null
+++ b/test/PoissonTests.jl
@@ -0,0 +1,81 @@
+module PoissonTests
+
+using Gridap
+using Gridap.Algebra
+using Gridap.FESpaces
+using GridapDistributed
+using PartitionedArrays
+using Test
+
+function main(parts)
+
+  output = mkpath(joinpath(@__DIR__,"output"))
+
+  domain = (0,4,0,4)
+  cells = (4,4)
+  model = CartesianDiscreteModel(parts,domain,cells)
+
+  labels = get_face_labeling(model)
+  add_tag_from_tags!(labels,"dirichlet",[1,2,3,5,7])
+  add_tag_from_tags!(labels,"neumann",[4,6,8])
+
+  Ω = Triangulation(model)
+  Γn = Boundary(model,tags="neumann")
+  n_Γn = get_normal_vector(Γn)
+
+  k = 2
+  u((x,y)) = (x+y)^k
+  f(x) = -Δ(u,x)
+  g = n_Γn⋅∇(u)
+
+  reffe = ReferenceFE(lagrangian,Float64,k)
+  V = TestFESpace(model,reffe,dirichlet_tags="dirichlet")
+  U = TrialFESpace(u,V)
+
+  dΩ = Measure(Ω,2*k)
+  dΓn = Measure(Γn,2*k)
+
+  a(u,v) = ∫( ∇(v)⋅∇(u) )dΩ
+  l(v) = ∫( v*f )dΩ + ∫( v*g )dΓn
+  op = AffineFEOperator(a,l,U,V)
+
+  solver = LinearFESolver(BackslashSolver())
+  uh = solve(solver,op)
+  eh = u - uh
+  @test sqrt(sum( ∫(abs2(eh))dΩ )) < 1.0e-9
+
+  # Now with DG
+  h = 2
+  γ = 10
+
+  V_dg = FESpace(model,reffe,conformity=:L2)
+
+  Λ = Skeleton(model)
+  Γd = Boundary(model,tags="dirichlet")
+
+  dΛ = Measure(Λ,2*k)
+  dΓd = Measure(Γd,2*k)
+
+  n_Γd = get_normal_vector(Γd)
+  n_Λ = get_normal_vector(Λ)
+
+  a_dg(u,v) =
+    ∫( ∇(v)⋅∇(u) )*dΩ +
+    ∫( (γ/h)*v*u  - v*(n_Γd⋅∇(u)) - (n_Γd⋅∇(v))*u )*dΓd +
+    ∫( (γ/h)*jump(v*n_Λ)⋅jump(u*n_Λ) -
+       jump(v*n_Λ)⋅mean(∇(u)) -
+       mean(∇(v))⋅jump(u*n_Λ) )*dΛ
+  
+  l_dg(v) =
+    ∫( v*f )*dΩ +
+    ∫( v*g )dΓn +
+    ∫( (γ/h)*v*u - (n_Γd⋅∇(v))*u )*dΓd
+
+  op = AffineFEOperator(a_dg,l_dg,V_dg,V_dg)
+  uh = solve(solver,op)
+  eh = u - uh
+  @test sqrt(sum( ∫(abs2(eh))dΩ )) < 1.0e-9
+
+end
+
+end # module
diff --git a/test/sequential/PoissonTests.jl b/test/sequential/PoissonTests.jl
new file mode 100644
index 00000000..2ada4780
--- /dev/null
+++ b/test/sequential/PoissonTests.jl
@@ -0,0 +1,12 @@
+module PoissonTestsSeq
+
+using PartitionedArrays
+
+include("../PoissonTests.jl")
+
+parts = get_part_ids(sequential,(2,2))
+PoissonTests.main(parts)
+
+end # module
+
+
diff --git a/test/sequential/runtests.jl b/test/sequential/runtests.jl
index 34e1ba92..884430d7 100644
--- a/test/sequential/runtests.jl
+++ b/test/sequential/runtests.jl
@@ -10,4 +10,6 @@ using Test
 
 @time @testset "MultiField" begin include("MultiFieldTests.jl") end
 
+@time @testset "Poisson" begin include("PoissonTests.jl") end
+
 end # module