diff --git a/NEWS.md b/NEWS.md
index 8f87937..2136fc0 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -7,12 +7,20 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.4.14] - 2026-03-20
+
+### Changed 
+
+- Added support for Gridap 0.20.0. Since PR[#200](https://github.com/gridap/GridapDistributed.jl/pull/200).
+
 ## [0.4.13] - 2026-03-16
 
 ### Fixed
+
 - BUG fix in `_generate_sign_flips(...)` private helper function for Nedelec elements. Since PR[#199](https://github.com/gridap/GridapDistributed.jl/pull/199).
 
 ### Added
+
 - An Hcurl projection test case. Since PR[#199](https://github.com/gridap/GridapDistributed.jl/pull/199).
 
 ## [0.4.12] - 2026-03-14
diff --git a/Project.toml b/Project.toml
index 0f6d581..a70000b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "GridapDistributed"
 uuid = "f9701e48-63b3-45aa-9a63-9bc6c271f355"
 authors = ["S. Badia <santiago.badia@monash.edu>", "A. F. Martin <alberto.f.martin@anu.edu.au>", "F. Verdugo <f.verdugo.rojano@vu.nl>"]
-version = "0.4.13"
+version = "0.4.14"
 
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
@@ -21,7 +21,7 @@ BlockArrays = "1"
 CircularArrays = "1.4.0"
 FillArrays = "1"
 ForwardDiff = "0.10, 1"
-Gridap = "0.19.9"
+Gridap = "0.20.0"
 LinearAlgebra = "1"
 MPI = "0.16, 0.17, 0.18, 0.19, 0.20"
 PartitionedArrays = "0.3.3"
diff --git a/src/DivAndCurlConformingFESpaces.jl b/src/DivAndCurlConformingFESpaces.jl
deleted file mode 100644
index 2217c07..0000000
--- a/src/DivAndCurlConformingFESpaces.jl
+++ /dev/null
@@ -1,177 +0,0 @@
-const RTorNedelec = Union{RaviartThomas,Nedelec}
-
-function FESpaces.FESpace(model::DistributedDiscreteModel,
-                          reffe::Tuple{<:RTorNedelec,Any,Any};
-                          conformity=nothing,
-                          split_own_and_ghost=false,
-                          constraint=nothing,
-                          kwargs...)
-
-  cell_reffes = map(local_views(model)) do m
-    basis,reffe_args,reffe_kwargs = reffe
-    cell_reffe = ReferenceFE(m,basis,reffe_args...;reffe_kwargs...)
-  end
-  _common_fe_space_constructor(model,cell_reffes;conformity,split_own_and_ghost,kwargs...)
-end
-
-function FESpaces.FESpace(model::DistributedDiscreteModel,
-                 reffe::GenericRefFE{<:RTorNedelec};
-                 conformity=nothing,
-                 split_own_and_ghost=false,
-                 constraint=nothing,
-                 kwargs...)
-  cell_reffes = map(local_views(model)) do m
-    Fill(reffe,num_cells(m))
-  end
-  _common_fe_space_constructor(model,cell_reffes;conformity,split_own_and_ghost,kwargs...)
-end
-
-function _setup_dmodel_and_dtrian(_trian::DistributedTriangulation)
-  trian = add_ghost_cells(_trian)
-  models = map(local_views(trian)) do t
-    get_active_model(t)
-  end
-  GenericDistributedDiscreteModel(models, generate_cell_gids(trian)), trian
-end  
-
-function FESpaces.FESpace(
-  _trian::DistributedTriangulation,
-  reffe::Tuple{<:RTorNedelec,Any,Any};
-  conformity=nothing,
-  split_own_and_ghost=false,
-  constraint=nothing,kwargs...
-)
-  dmodel, dtrian = _setup_dmodel_and_dtrian(_trian)
-  cell_reffes = map(local_views(dmodel)) do m
-    basis,reffe_args,reffe_kwargs = reffe
-    cell_reffe = ReferenceFE(m,basis,reffe_args...;reffe_kwargs...)
-  end
-  _common_fe_space_constructor(dmodel,cell_reffes,dtrian;conformity,split_own_and_ghost,kwargs...)
-end
-
-function FESpaces.FESpace(_trian::DistributedTriangulation,
-                 reffe::GenericRefFE{<:RTorNedelec};
-                 conformity=nothing,
-                 split_own_and_ghost=false,
-                 constraint=nothing,
-                 kwargs...)
-  dmodel, dtrian = _setup_dmodel_and_dtrian(_trian)
-  cell_reffes = map(local_views(dmodel)) do m
-     Fill(reffe,num_cells(m))
-  end
-  _common_fe_space_constructor(dmodel,cell_reffes,dtrian;conformity,split_own_and_ghost,kwargs...)
-end
-
-function _common_fe_space_constructor(model,cell_reffes,trian;conformity,split_own_and_ghost,kwargs...)
-  sign_flips=_generate_sign_flips(model,cell_reffes)
-  spaces = map(local_views(model),sign_flips,cell_reffes,local_views(trian)) do m,sign_flip,cell_reffe,trian
-     conf = Conformity(testitem(cell_reffe),conformity)
-     cell_fe = CellFE(m,cell_reffe,conf,sign_flip)
-     FESpace(m, cell_fe; trian=trian, kwargs...)
-  end
-  gids = generate_gids(model,spaces)
-  vector_type = _find_vector_type(spaces,gids;split_own_and_ghost=split_own_and_ghost)
-  DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
-end
-
-function _common_fe_space_constructor(model,cell_reffes;conformity,split_own_and_ghost,kwargs...)
-  sign_flips=_generate_sign_flips(model,cell_reffes)
-  spaces = map(local_views(model),sign_flips,cell_reffes) do m,sign_flip,cell_reffe
-     conf = Conformity(testitem(cell_reffe),conformity)
-     cell_fe = CellFE(m,cell_reffe,conf,sign_flip)
-     FESpace(m, cell_fe; kwargs...)
-  end
-  gids = generate_gids(model,spaces)
-  trian = DistributedTriangulation(map(get_triangulation,spaces),model)
-  vector_type = _find_vector_type(spaces,gids;split_own_and_ghost=split_own_and_ghost)
-  DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
-end
-
-function _generate_sign_flips(model, cell_gids, cell_reffes)
-sign_flips = map(local_views(model),partition(cell_gids),cell_reffes) do m, p, cell_reffe
-    D = num_cell_dims(model)
-
-    gtopo = get_grid_topology(m)
-
-    # Extract composition among cells and facets
-    cell_wise_facets_ids = get_faces(gtopo, D, D - 1)
-    cache_cell_wise_facets_ids = array_cache(cell_wise_facets_ids)
-
-    # Extract cells around facets
-    cells_around_facets = get_faces(gtopo, D - 1, D)
-    cache_cells_around_facets = array_cache(cells_around_facets)
-
-    ncells = num_cells(m)
-    ptrs = Vector{Int32}(undef,ncells+1)
-    for cell in 1:ncells
-      reffe=cell_reffe[cell]
-      ptrs[cell+1] = num_dofs(reffe)
-    end
-    PArrays.length_to_ptrs!(ptrs)
-    ndata = ptrs[end]-1
-    data  = Vector{Bool}(undef,ndata)
-    data .= false
-
-    loc_to_glo=local_to_global(p)
-    for cell in own_to_local(p)
-      sign_flip = view(data,ptrs[cell]:ptrs[cell+1]-1)
-      reffe = cell_reffe[cell]
-      D = num_dims(reffe)
-      face_own_dofs   = get_face_own_dofs(reffe)
-      facet_lid       = get_offsets(get_polytope(reffe))[D] + 1
-      cell_facets_ids = getindex!(cache_cell_wise_facets_ids,
-                                  cell_wise_facets_ids,
-                                  cell)
-      for facet_gid in cell_facets_ids
-          facet_cells_around = getindex!(cache_cells_around_facets,
-                cells_around_facets,
-                facet_gid)
-          is_slave=false
-          if length(facet_cells_around) > 1
-            mx=maximum(loc_to_glo[facet_cells_around])
-            is_slave = (loc_to_glo[cell] == mx)
-          end
-          if is_slave
-              for dof in face_own_dofs[facet_lid]
-                  sign_flip[dof] = true
-              end
-          end
-          facet_lid = facet_lid + 1
-      end
-    end
-    JaggedArray(data,ptrs)
-  end
-  cache = fetch_vector_ghost_values_cache(sign_flips,partition(cell_gids))
-  fetch_vector_ghost_values!(sign_flips,cache) |> wait
-  sign_flips
-end 
-
-
-function _generate_sign_flips(model::DistributedDiscreteModel,
-                             cell_reffes::AbstractVector{<:AbstractVector{<:GenericRefFE{<:RaviartThomas}}})
-  cell_gids  = get_cell_gids(model)
-  _generate_sign_flips(model,cell_gids,cell_reffes)
-end
-
-function _generate_sign_flips(model::DistributedDiscreteModel,
-                              cell_reffes::AbstractVector{<:AbstractVector{<:GenericRefFE{<:Nedelec}}})
-  Dc=num_cell_dims(model)
-  is_simp = false 
-  map(cell_reffes) do cell_reffes
-    cell_reffe = cell_reffes[1]
-    polytope = get_polytope(cell_reffe)
-    if (Gridap.Geometry.is_simplex(polytope))
-        is_simp=true
-    end 
-  end 
-  if (Dc==2 || is_simp)
-    map(local_views(model), cell_reffes) do model, cell_reffes
-        Gridap.FESpaces._no_sign_flip(model, cell_reffes)
-    end
-  else
-    @assert Dc==3
-    cell_gids  = get_cell_gids(model)
-    _generate_sign_flips(model,cell_gids,cell_reffes)
-  end
-end
-
diff --git a/src/FESpaces.jl b/src/FESpaces.jl
index b335447..74c16b2 100644
--- a/src/FESpaces.jl
+++ b/src/FESpaces.jl
@@ -437,22 +437,26 @@ function FESpaces.HomogeneousTrialFESpace(f::DistributedSingleFieldFESpace)
 end
 
 function generate_gids(
-  model::DistributedDiscreteModel{Dc},
+  model::DistributedDiscreteModel,
   spaces::AbstractArray{<:SingleFieldFESpace}
-) where Dc
-  cell_to_ldofs = map(get_cell_dof_ids,spaces)
-  nldofs = map(num_free_dofs,spaces)
+)
   cell_gids = get_cell_gids(model)
-  generate_gids(cell_gids,cell_to_ldofs,nldofs)
+  generate_gids(cell_gids,spaces)
 end
 
 function generate_gids(
-  trian::DistributedTriangulation{Dc},
+  trian::DistributedTriangulation,
   spaces::AbstractArray{<:SingleFieldFESpace}
-) where Dc
+)
+  cell_gids = generate_cell_gids(trian)
+  generate_gids(cell_gids,spaces)
+end
+
+function generate_gids(
+  cell_gids::PRange, spaces::AbstractArray{<:SingleFieldFESpace}
+)
   cell_to_ldofs = map(get_cell_dof_ids,spaces)
   nldofs = map(num_free_dofs,spaces)
-  cell_gids = generate_cell_gids(trian)
   generate_gids(cell_gids,cell_to_ldofs,nldofs)
 end
 
@@ -537,30 +541,76 @@ end
 
 # Factories
 
+# function FESpaces.FESpace(
+#   model::DistributedDiscreteModel,reffe;split_own_and_ghost=false,constraint=nothing,kwargs...
+# )
+#   spaces = map(local_views(model)) do m
+#     FESpace(m,reffe;kwargs...)
+#   end
+#   gids =  generate_gids(model,spaces)
+#   trian = DistributedTriangulation(map(get_triangulation,spaces),model)
+#   vector_type = _find_vector_type(spaces,gids;split_own_and_ghost=split_own_and_ghost)
+#   space = DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
+#   return _add_distributed_constraint(space,reffe,constraint)
+# end
+# 
+# function FESpaces.FESpace(
+#   _trian::DistributedTriangulation,reffe;split_own_and_ghost=false,constraint=nothing,kwargs...
+# )
+#   trian = add_ghost_cells(_trian)
+#   spaces = map(local_views(trian)) do t
+#     FESpace(t,reffe;kwargs...)
+#   end
+#   gids = generate_gids(trian,spaces)
+#   vector_type = _find_vector_type(spaces,gids;split_own_and_ghost=split_own_and_ghost)
+#   space = DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
+#   return _add_distributed_constraint(space,reffe,constraint)
+# end
+
 function FESpaces.FESpace(
-  model::DistributedDiscreteModel,reffe;split_own_and_ghost=false,constraint=nothing,kwargs...
+  model::DistributedDiscreteModel,args...;kwargs...
 )
-  spaces = map(local_views(model)) do m
-    FESpace(m,reffe;kwargs...)
-  end
-  gids =  generate_gids(model,spaces)
-  trian = DistributedTriangulation(map(get_triangulation,spaces),model)
-  vector_type = _find_vector_type(spaces,gids;split_own_and_ghost=split_own_and_ghost)
-  space = DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
-  return _add_distributed_constraint(space,reffe,constraint)
+  trian = Triangulation(with_ghost,model)
+  cell_gids = get_cell_gids(model)
+  DistributedSingleFieldFESpace(model,trian,cell_gids,args...;kwargs...)
 end
 
 function FESpaces.FESpace(
-  _trian::DistributedTriangulation,reffe;split_own_and_ghost=false,constraint=nothing,kwargs...
+  _trian::DistributedTriangulation,args...;kwargs...
 )
   trian = add_ghost_cells(_trian)
-  spaces = map(local_views(trian)) do t
-    FESpace(t,reffe;kwargs...)
+  cell_gids = generate_cell_gids(trian)
+  model = DistributedDiscreteModel(map(get_active_model,local_views(trian)), cell_gids)
+  DistributedSingleFieldFESpace(model,trian,cell_gids,args...;kwargs...)
+end
+
+function DistributedSingleFieldFESpace(
+  model::DistributedDiscreteModel,
+  trian::DistributedTriangulation,
+  cell_gids::PRange, reffe; kwargs...
+)
+  cell_reffe = map(local_views(model)) do model
+    ReferenceFE(model,reffe)
   end
-  gids = generate_gids(trian,spaces)
-  vector_type = _find_vector_type(spaces,gids;split_own_and_ghost=split_own_and_ghost)
+  DistributedSingleFieldFESpace(model,trian,cell_gids,cell_reffe;kwargs...)
+end
+
+function DistributedSingleFieldFESpace(
+  model::DistributedDiscreteModel, # Active model, not bg model
+  trian::DistributedTriangulation,
+  cell_gids::PRange, 
+  cell_reffe::AbstractArray; 
+  split_own_and_ghost=false, 
+  constraint=nothing,
+  kwargs...
+)
+  spaces = map(local_views(model),local_views(trian),cell_reffe) do model, trian, cell_reffe
+    FESpace(model,cell_reffe;trian,kwargs...)
+  end
+  gids = generate_gids(cell_gids,spaces)
+  vector_type = _find_vector_type(spaces,gids;split_own_and_ghost)
   space = DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
-  return _add_distributed_constraint(space,reffe,constraint)
+  return _add_distributed_constraint(space,cell_reffe,constraint)
 end
 
 function _find_vector_type(spaces,gids;split_own_and_ghost=false)
@@ -595,6 +645,16 @@ function _add_distributed_constraint(
   _add_distributed_constraint(F,order,constraint)
 end
 
+function _add_distributed_constraint(
+  F::DistributedFESpace,cell_reffe::AbstractArray,constraint
+)
+  reffe = map(cell_reffe) do cell_reffe
+    reffes, ctypes = compress_cell_data(cell_reffe)
+    return only(reffes)
+  end |> getany
+  _add_distributed_constraint(F,reffe,constraint)
+end
+
 function _add_distributed_constraint(F::DistributedFESpace,order::Integer,constraint)
   if isnothing(constraint)
     V = F
diff --git a/src/GridapDistributed.jl b/src/GridapDistributed.jl
index cfe4135..da79999 100644
--- a/src/GridapDistributed.jl
+++ b/src/GridapDistributed.jl
@@ -58,7 +58,7 @@ include("Visualization.jl")
 
 include("FESpaces.jl")
 
-include("DivAndCurlConformingFESpaces.jl")
+include("Pullbacks.jl")
 
 include("MultiField.jl")
 
diff --git a/src/Pullbacks.jl b/src/Pullbacks.jl
new file mode 100644
index 0000000..d8d0bba
--- /dev/null
+++ b/src/Pullbacks.jl
@@ -0,0 +1,137 @@
+
+# These should be the reffes that require a 
+# globally-computed pushforward, and thus cannot 
+# simply go through the local FESpace constructor.
+const PullbackReffes = Union{
+  GenericRefFE{RaviartThomas},
+  GenericRefFE{Nedelec},
+}
+
+function DistributedSingleFieldFESpace(
+  model::DistributedDiscreteModel, # Active model, not bg model
+  trian::DistributedTriangulation,
+  cell_gids::PRange, 
+  cell_reffe::AbstractArray{<:AbstractArray{T}};
+  labels = get_face_labeling(model), 
+  split_own_and_ghost=false, 
+  constraint=nothing,
+  conformity=nothing,
+  scale_dof=false,
+  global_meshsize=nothing,
+  kwargs...
+) where T <: PullbackReffes
+  # Construct a globally conforming CellFE
+  conf = map(cell_reffe) do cell_reffe
+    Conformity(testitem(cell_reffe),conformity)
+  end |> getany
+  cell_fe = FESpaces.CellFE(model, cell_reffe, conf; scale_dof, global_meshsize)
+
+  spaces = map(
+    local_views(model),local_views(trian),local_views(labels), cell_fe
+  ) do model, trian, labels, cell_fe
+    FESpace(model,cell_fe;trian,labels,kwargs...)
+  end
+
+  gids = generate_gids(cell_gids,spaces)
+  vector_type = _find_vector_type(spaces,gids;split_own_and_ghost)
+  space = DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
+  return _add_distributed_constraint(space,cell_reffe,constraint)
+end
+
+function FESpaces.CellFE(
+  model::DistributedDiscreteModel, cell_reffe::AbstractArray{<:AbstractArray{T}}, conformity::Conformity; kwargs...
+) where T <: ReferenceFE
+  cell_bases = FESpaces.get_cell_shapefuns_and_dof_basis(
+    model, cell_reffe, conformity; kwargs...
+  )
+  cell_fe = map(local_views(model), cell_reffe, cell_bases) do model, cell_reffe, cell_bases
+    cell_shapefuns, cell_dof_basis = cell_bases
+    FESpaces.CellFE(model, cell_reffe, cell_shapefuns, cell_dof_basis, conformity)
+  end
+  return cell_fe
+end
+
+function FESpaces.get_cell_shapefuns_and_dof_basis(
+  model::DistributedDiscreteModel, cell_reffe::AbstractArray{<:AbstractArray{T}}, conformity::Conformity; kwargs...
+) where T <: ReferenceFE
+  reffe_name = get_name(T)
+  pushforward = Pushforward(reffe_name, conformity)
+  cell_Jt = map(local_views(model)) do model
+    cell_map = get_cell_map(get_grid(model))
+    return lazy_map(Broadcasting(∇), cell_map)
+  end
+  # The cell changes are the ones that need to be globally consistent
+  # They contain sign changes, etc...
+  cell_changes = FESpaces.compute_cell_bases_changes(
+    reffe_name, pushforward, model, cell_reffe, cell_Jt
+  )
+  cell_bases = map(
+    local_views(model), cell_reffe, cell_changes, cell_Jt
+  ) do model, cell_reffe, cell_changes, cell_Jt
+    FESpaces.get_cell_shapefuns_and_dof_basis(
+      pushforward, model, cell_reffe, cell_changes, cell_Jt; kwargs...
+    )
+  end
+  return cell_bases
+end
+
+function FESpaces.compute_cell_bases_changes(
+  ::ReferenceFEName, ::ContraVariantPiolaMap, model::DistributedDiscreteModel, cell_reffe, cell_Jt
+)
+  change = FESpaces.get_sign_flip(model, cell_reffe) # equal to its transposed inverse
+  return map(c -> (c, c), change)
+end
+
+function FESpaces.compute_cell_bases_changes(
+  ::ReferenceFEName, ::CoVariantPiolaMap, model::DistributedDiscreteModel, cell_reffe, cell_Jt
+)
+  D = num_cell_dims(model)
+  poly = map(get_polytopes, local_views(model)) |> getany |> only
+  if (D==2) || is_simplex(poly)
+    # For these cases, we do not need to apply a sign flip
+    return nothing
+  elseif (D==3) && is_n_cube(poly)
+    change = FESpaces.get_sign_flip(model, cell_reffe)
+    return map(c -> (c, c), change)
+  end
+  @notimplemented
+end
+
+function FESpaces.get_sign_flip(model::DistributedDiscreteModel, cell_reffe)
+  facet_owners = FESpaces.compute_facet_owners(model)
+  map(local_views(model),cell_reffe,facet_owners) do model, cell_reffe, facet_owners
+    sign_map = FESpaces.NormalSignMap(model,facet_owners)
+    FESpaces.get_sign_flip(model, cell_reffe, sign_map)
+  end
+end
+
+function FESpaces.compute_facet_owners(model::DistributedDiscreteModel)
+  Dc = num_cell_dims(model)
+  cell_ids  = partition(get_cell_gids(model))
+  facet_ids = partition(get_face_gids(model, Dc-1))
+  facet_to_owner = map(FESpaces.compute_facet_owners, local_views(model))
+
+  # Map local owners to global ids
+  map(facet_to_owner, cell_ids) do facet_to_owner, cell_ids
+    lid_to_gid = local_to_global(cell_ids)
+    for f in eachindex(facet_to_owner)
+      lid = facet_to_owner[f]
+      facet_to_owner[f] = lid_to_gid[lid]
+    end
+  end
+
+  # Communicate true owners across processes
+  wait(consistent!(PVector(facet_to_owner, facet_ids)))
+
+  # Non-local owners will be set to zero, which 
+  # will trigger a sign flip (which is the correct behaviour)
+  map(facet_to_owner, cell_ids) do facet_to_owner, cell_ids
+    gid_to_lid = global_to_local(cell_ids)
+    for f in eachindex(facet_to_owner)
+      gid = facet_to_owner[f]
+      facet_to_owner[f] = gid_to_lid[gid]
+    end
+  end
+
+  return facet_to_owner
+end
diff --git a/test/DivAndCurlConformingTests.jl b/test/DivAndCurlConformingTests.jl
deleted file mode 100644
index 85c88e0..0000000
--- a/test/DivAndCurlConformingTests.jl
+++ /dev/null
@@ -1,148 +0,0 @@
-module DivAndCurlConformingTests
-using SparseMatricesCSR
-import Gridap: ∇, divergence, DIV
-using Gridap
-using Gridap.Algebra
-using Gridap.FESpaces
-using GridapDistributed
-using PartitionedArrays
-using Test
-using FillArrays
-
-function setup_p1_model()
-  ptr  = [ 1, 5, 9 ]
-  data = [ 1,2,3,4, 2,5,4,6  ]
-  cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
-  node_coordinates = Vector{Point{2,Float64}}(undef,6)
-
-  node_coordinates[1]=Point{2,Float64}(0.0,0.0)
-  node_coordinates[2]=Point{2,Float64}(0.5,0.0)
-  node_coordinates[3]=Point{2,Float64}(0.0,1.0)
-  node_coordinates[4]=Point{2,Float64}(0.5,1.0)
-  node_coordinates[5]=Point{2,Float64}(1.0,0.0)
-  node_coordinates[6]=Point{2,Float64}(1.0,1.0)
-
-  polytope=QUAD
-  scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1)
-  cell_types=collect(Fill(1,length(cell_vertex_lids)))
-  cell_reffes=[scalar_reffe]
-  grid = Gridap.Geometry.UnstructuredGrid(node_coordinates,
-                                          cell_vertex_lids,
-                                          cell_reffes,
-                                          cell_types,
-                                          Gridap.Geometry.NonOriented())
-  Gridap.Geometry.UnstructuredDiscreteModel(grid)
-end
-
-function setup_p2_model()
-  ptr  = [ 1, 5, 9 ]
-  data = [ 1,2,3,4, 5,1,6,3  ]
-  cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
-  node_coordinates = Vector{Point{2,Float64}}(undef,6)
-
-  node_coordinates[1]=Point{2,Float64}(0.5,0.0)
-  node_coordinates[2]=Point{2,Float64}(1.0,0.0)
-  node_coordinates[3]=Point{2,Float64}(0.5,1.0)
-  node_coordinates[4]=Point{2,Float64}(1.0,1.0)
-  node_coordinates[5]=Point{2,Float64}(0.0,0.0)
-  node_coordinates[6]=Point{2,Float64}(0.0,1.0)
-
-  polytope=QUAD
-  scalar_reffe=Gridap.ReferenceFEs.ReferenceFE(polytope,Gridap.ReferenceFEs.lagrangian,Float64,1)
-  cell_types=collect(Fill(1,length(cell_vertex_lids)))
-  cell_reffes=[scalar_reffe]
-  grid = Gridap.Geometry.UnstructuredGrid(node_coordinates,
-                                          cell_vertex_lids,
-                                          cell_reffes,
-                                          cell_types,
-                                          Gridap.Geometry.NonOriented())
-  Gridap.Geometry.UnstructuredDiscreteModel(grid)
-end
-
-function f(model,reffe,trian,das)
-    V = FESpace(model,reffe)
-    U = TrialFESpace(V)
-
-    degree = 2
-    dΩ = Measure(trian,degree)
-    a(u,v) = ∫( u⋅v )*dΩ
-
-    u  = get_trial_fe_basis(U)
-    v  = get_fe_basis(V)
-    dc = a(u,v)
-
-    assem = SparseMatrixAssembler(U,V,das)
-    data = collect_cell_matrix(U,V,a(u,v))
-    A = assemble_matrix(assem,data)
-    t1  = trian.trians.items[1]
-    t2  = trian.trians.items[2]
-    dc1 = dc.contribs.items[1]
-    dc2 = dc.contribs.items[2]
-    c1  = Gridap.CellData.get_contribution(dc1,t1)
-    c2  = Gridap.CellData.get_contribution(dc2,t2)
-    
-    tol = 1.0e-12
-    @test norm(c1[1]-c2[2]) < tol
-    @test norm(c1[2]-c2[1]) < tol
-end
-
-function main(distribute,nranks)
-   ranks = distribute(LinearIndices((prod(nranks),)))
-   @assert isa(ranks,DebugArray)
-   @assert length(ranks)==2
-
-   models=map(ranks) do rank
-      if (rank==1)
-        setup_p1_model()
-      else
-        setup_p2_model()
-      end
-    end
-
-    num_owners,owner_to_global,ghost_to_global,ghost_to_owner=map(ranks) do rank
-      if (rank==1)
-        1,[1],[2],Int32[2]
-      else
-        1,[2],[1],Int32[1]
-      end
-    end |> tuple_of_arrays
-
-    indices=map(ranks, owner_to_global, ghost_to_global, ghost_to_owner) do rank, 
-                                                                            owner_to_global,
-                                                                            ghost_to_global,
-                                                                            ghost_to_owner
-      owner = rank
-      ngdofs = 2
-      own_indices=OwnIndices(ngdofs,owner,owner_to_global)
-      ghost_indices=GhostIndices(ngdofs,ghost_to_global,ghost_to_owner)
-      OwnAndGhostIndices(own_indices,ghost_indices)
-    end
-    gids=PRange(indices)
-
-    model = GridapDistributed.DistributedDiscreteModel(models,gids)
-
-    das = FullyAssembledRows()
-    trian = Triangulation(das,model)
-
-    reffe=ReferenceFE(raviart_thomas,Float64,0)
-    f(model,reffe,trian,das)
-    f(trian,reffe,trian,das)
-    f(Triangulation(model),reffe,trian,das)
-    
-    reffe=ReferenceFE(QUAD, raviart_thomas, 0)
-    f(model,reffe,trian,das)
-    f(trian,reffe,trian,das)
-    f(Triangulation(model),reffe,trian,das)
-
-    reffe=ReferenceFE(nedelec,Float64,0)
-    f(model,reffe,trian,das)
-    f(trian,reffe,trian,das)
-    f(Triangulation(model),reffe,trian,das)
-    
-    reffe=ReferenceFE(QUAD, nedelec, 0)
-    f(model,reffe,trian,das)
-    f(trian,reffe,trian,das)
-    f(Triangulation(model),reffe,trian,das)
-  end
-
-end # module
diff --git a/test/HcurlProjectionTests.jl b/test/HcurlProjectionTests.jl
index 671be7c..b59dc71 100644
--- a/test/HcurlProjectionTests.jl
+++ b/test/HcurlProjectionTests.jl
@@ -20,78 +20,77 @@ function get_analytical_functions(Dc)
 end
 
 function solve_hcurl_projection(model::GridapDistributed.DistributedDiscreteModel{Dc},order) where {Dc}
-    u_ex, f_ex=get_analytical_functions(Dc)
+  u_ex, f_ex=get_analytical_functions(Dc)
 
-    V = FESpace(model,
-                ReferenceFE(nedelec,order),
-                conformity=:Hcurl,
-                dirichlet_tags="boundary")
-    
-    U = TrialFESpace(V,u_ex)
-    
-    trian = Triangulation(model)
-    degree = 2*(order+1)
-    dΩ = Measure(trian,degree)
-        
-    a(u,v) = ∫( (∇×u)⋅(∇×v) + u⋅v )dΩ
-    l(v) = ∫(f_ex⋅v)dΩ
+  V = FESpace(
+    model, ReferenceFE(nedelec,order), conformity=:Hcurl, dirichlet_tags="boundary"
+  )
+  
+  U = TrialFESpace(V,u_ex)
+  
+  trian = Triangulation(model)
+  degree = 2*(order+1)
+  dΩ = Measure(trian,degree)
+      
+  a(u,v) = ∫( (∇×u)⋅(∇×v) + u⋅v )dΩ
+  l(v) = ∫(f_ex⋅v)dΩ
 
-    op = AffineFEOperator(a,l,U,V)
-    if (num_free_dofs(U)==0)
-      # UMFPACK cannot handle empty linear systems
-      uh = zero(U)
-    else
-      uh = solve(op)
-    end
-    uh,U
-  end 
+  op = AffineFEOperator(a,l,U,V)
+  if (num_free_dofs(U)==0)
+    # UMFPACK cannot handle empty linear systems
+    uh = zero(U)
+  else
+    uh = solve(op)
+  end
+  uh,U
+end 
 
-  function check_error_hcurl_projection(model::GridapDistributed.DistributedDiscreteModel{Dc},order,uh) where {Dc}
-    trian = Triangulation(model)
-    degree = 2*(order+1)
-    dΩ = Measure(trian,degree)
+function check_error_hcurl_projection(model::GridapDistributed.DistributedDiscreteModel{Dc},order,uh) where {Dc}
+  trian = Triangulation(model)
+  degree = 2*(order+1)
+  dΩ = Measure(trian,degree)
 
-    u_ex, f_ex = get_analytical_functions(Dc)
-    
-    eu = u_ex - uh
+  u_ex, f_ex = get_analytical_functions(Dc)
+  
+  eu = u_ex - uh
 
-    l2(v) = sqrt(sum(∫(v⋅v)*dΩ))
-    hcurl(v) = sqrt(sum(∫(v⋅v + (∇×v)⋅(∇×v))*dΩ))
-    
-    eu_l2 = l2(eu)
-    eu_hcurl = hcurl(eu)
-    
-    tol = 1.0e-6
-    @test eu_l2 < tol
-    @test eu_hcurl < tol
+  l2(v) = sqrt(sum(∫(v⋅v)*dΩ))
+  hcurl(v) = sqrt(sum(∫(v⋅v + (∇×v)⋅(∇×v))*dΩ))
+  
+  eu_l2 = l2(eu)
+  eu_hcurl = hcurl(eu)
+  
+  tol = 1.0e-6
+  @test eu_l2 < tol
+  @test eu_hcurl < tol
 end
 
- function test_2d(ranks,parts,order)
-    domain = (0,1,0,1)
-    model = CartesianDiscreteModel(ranks,parts,domain,(4,4))
-    solve_hcurl_projection(model,order) |> x -> check_error_hcurl_projection(model,order,x[1])
- end 
+function test_2d(ranks,parts,order)
+  domain = (0,1,0,1)
+  model = CartesianDiscreteModel(ranks,parts,domain,(4,4))
+  solve_hcurl_projection(model,order) |> x -> check_error_hcurl_projection(model,order,x[1])
+end 
 
- function test_3d(ranks,parts,order)
-    domain = (0,1,0,1,0,1)
-    model = CartesianDiscreteModel(ranks,parts,domain,(4,4,4))
-    solve_hcurl_projection(model,order) |> x -> check_error_hcurl_projection(model,order,x[1])
- end 
+function test_3d(ranks,parts,order)
+  domain = (0,1,0,1,0,1)
+  model = CartesianDiscreteModel(ranks,parts,domain,(4,4,4))
+  solve_hcurl_projection(model,order) |> x -> check_error_hcurl_projection(model,order,x[1])
+end 
 
 function main(distribute,parts)
   ranks = distribute(LinearIndices((prod(parts),)))
   
   if length(parts)==2 
-      for order=0:2
-        test_2d(ranks,parts,order)
-      end
+    for order=0:2
+      test_2d(ranks,parts,order)
+    end
   elseif length(parts)==3
-      for order=0:2
-        test_3d(ranks,parts,order)
-      end
+    for order=0:2
+      test_3d(ranks,parts,order)
+    end
   else 
       @assert false 
   end 
 end
-  
+
 end #module
diff --git a/test/PullbackTests.jl b/test/PullbackTests.jl
new file mode 100644
index 0000000..3c12f37
--- /dev/null
+++ b/test/PullbackTests.jl
@@ -0,0 +1,143 @@
+module PullbackTests
+using SparseMatricesCSR
+
+using Gridap
+using Gridap.Geometry
+using Gridap.Algebra
+using Gridap.FESpaces
+using Gridap.ReferenceFEs
+using GridapDistributed
+using PartitionedArrays
+using Test
+using FillArrays
+
+using Gridap: ∇, divergence, DIV
+
+function setup_p1_model()
+  ptr  = [ 1, 5, 9 ]
+  data = [ 1,2,3,4, 2,5,4,6  ]
+  cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
+
+  node_coordinates = Vector{Point{2,Float64}}(undef,6)
+  node_coordinates[1]=Point{2,Float64}(0.0,0.0)
+  node_coordinates[2]=Point{2,Float64}(0.5,0.0)
+  node_coordinates[3]=Point{2,Float64}(0.0,1.0)
+  node_coordinates[4]=Point{2,Float64}(0.5,1.0)
+  node_coordinates[5]=Point{2,Float64}(1.0,0.0)
+  node_coordinates[6]=Point{2,Float64}(1.0,1.0)
+
+  reffe = ReferenceFEs.ReferenceFE(QUAD,ReferenceFEs.lagrangian,Float64,1)
+  cell_types = collect(Fill(1,length(cell_vertex_lids)))
+  grid = Geometry.UnstructuredGrid(
+    node_coordinates, cell_vertex_lids, [reffe], cell_types, Geometry.NonOriented()
+  )
+  Geometry.UnstructuredDiscreteModel(grid)
+end
+
+function setup_p2_model()
+  ptr  = [ 1, 5, 9 ]
+  data = [ 1,2,3,4, 5,1,6,3  ]
+  cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
+  
+  node_coordinates = Vector{Point{2,Float64}}(undef,6)
+  node_coordinates[1]=Point{2,Float64}(0.5,0.0)
+  node_coordinates[2]=Point{2,Float64}(1.0,0.0)
+  node_coordinates[3]=Point{2,Float64}(0.5,1.0)
+  node_coordinates[4]=Point{2,Float64}(1.0,1.0)
+  node_coordinates[5]=Point{2,Float64}(0.0,0.0)
+  node_coordinates[6]=Point{2,Float64}(0.0,1.0)
+
+  reffe = Gridap.ReferenceFEs.ReferenceFE(QUAD,ReferenceFEs.lagrangian,Float64,1)
+  cell_types = collect(Fill(1,length(cell_vertex_lids)))
+  grid = Geometry.UnstructuredGrid(
+    node_coordinates, cell_vertex_lids, [reffe], cell_types, Geometry.NonOriented()
+  )
+  return Geometry.UnstructuredDiscreteModel(grid)
+end
+
+function f(model,reffe,trian,das)
+  V = FESpace(model,reffe)
+  U = TrialFESpace(V)
+
+  degree = 2
+  dΩ = Measure(trian,degree)
+  a(u,v) = ∫( u⋅v )*dΩ
+
+  u  = get_trial_fe_basis(U)
+  v  = get_fe_basis(V)
+  dc = a(u,v)
+
+  assem = SparseMatrixAssembler(U,V,das)
+  data = collect_cell_matrix(U,V,a(u,v))
+  A = assemble_matrix(assem,data)
+  t1  = trian.trians.items[1]
+  t2  = trian.trians.items[2]
+  dc1 = dc.contribs.items[1]
+  dc2 = dc.contribs.items[2]
+  c1  = Gridap.CellData.get_contribution(dc1,t1)
+  c2  = Gridap.CellData.get_contribution(dc2,t2)
+  
+  tol = 1.0e-12
+  @test norm(c1[1]-c2[2]) < tol
+  @test norm(c1[2]-c2[1]) < tol
+end
+
+function main(distribute,nranks)
+  ranks = distribute(LinearIndices((prod(nranks),)))
+  @assert isa(ranks,DebugArray)
+  @assert length(ranks)==2
+
+  models = map(ranks) do rank
+    if (rank==1)
+      setup_p1_model()
+    else
+      setup_p2_model()
+    end
+  end
+
+  num_owners, owner_to_global, ghost_to_global, ghost_to_owner = map(ranks) do rank
+    if (rank==1)
+      1,[1],[2],Int32[2]
+    else
+      1,[2],[1],Int32[1]
+    end
+  end |> tuple_of_arrays
+
+  indices = map(
+    ranks, owner_to_global, ghost_to_global, ghost_to_owner
+    ) do rank, owner_to_global, ghost_to_global, ghost_to_owner
+    owner = rank
+    ngdofs = 2
+    own_indices=OwnIndices(ngdofs,owner,owner_to_global)
+    ghost_indices=GhostIndices(ngdofs,ghost_to_global,ghost_to_owner)
+    OwnAndGhostIndices(own_indices,ghost_indices)
+  end
+  gids = PRange(indices)
+
+  model = GridapDistributed.DistributedDiscreteModel(models,gids)
+
+  das = FullyAssembledRows()
+  trian = Triangulation(das,model)
+
+  reffe = ReferenceFE(raviart_thomas,Float64,0)
+  f(model,reffe,trian,das)
+  f(trian,reffe,trian,das)
+  f(Triangulation(model),reffe,trian,das)
+  
+  reffe = ReferenceFE(QUAD, raviart_thomas, 0)
+  f(model,reffe,trian,das)
+  f(trian,reffe,trian,das)
+  f(Triangulation(model),reffe,trian,das)
+
+  reffe = ReferenceFE(nedelec,Float64,0)
+  f(model,reffe,trian,das)
+  f(trian,reffe,trian,das)
+  f(Triangulation(model),reffe,trian,das)
+  
+  reffe = ReferenceFE(QUAD, nedelec, 0)
+  f(model,reffe,trian,das)
+  f(trian,reffe,trian,das)
+  f(Triangulation(model),reffe,trian,das)
+end
+
+end # module
diff --git a/test/sequential/DivAndCurlConformingTests.jl b/test/sequential/DivAndCurlConformingTests.jl
deleted file mode 100644
index e317b58..0000000
--- a/test/sequential/DivAndCurlConformingTests.jl
+++ /dev/null
@@ -1,9 +0,0 @@
-module DivConformingTestsSeq
-using PartitionedArrays
-include("../DivAndCurlConformingTests.jl")
-
-with_debug() do distribute
-  DivAndCurlConformingTests.main(distribute,2)
-end
-
-end # module
diff --git a/test/sequential/FESpacesTests.jl b/test/sequential/FESpacesTests.jl
index 3463c83..8753797 100644
--- a/test/sequential/FESpacesTests.jl
+++ b/test/sequential/FESpacesTests.jl
@@ -1,7 +1,9 @@
 module FESpacesTestsSeq
 using PartitionedArrays
+
 include("../FESpacesTests.jl")
 with_debug() do distribute
   FESpacesTests.main(distribute,(2,2))
 end
+
 end # module
diff --git a/test/sequential/PullbackTests.jl b/test/sequential/PullbackTests.jl
new file mode 100644
index 0000000..8f224fc
--- /dev/null
+++ b/test/sequential/PullbackTests.jl
@@ -0,0 +1,9 @@
+module PullbackTestsSeq
+using PartitionedArrays
+
+include("../PullbackTests.jl")
+with_debug() do distribute
+  PullbackTests.main(distribute,2)
+end
+
+end # module
\ No newline at end of file
diff --git a/test/sequential/runtests.jl b/test/sequential/runtests.jl
index eedc5a3..aee13a6 100644
--- a/test/sequential/runtests.jl
+++ b/test/sequential/runtests.jl
@@ -11,6 +11,7 @@ end
 
 if TESTCASE ∈ ("all", "seq-fespaces")
   @time @testset "FESpaces"         begin include("FESpacesTests.jl") end
+  @time @testset "Pullbacks"        begin include("PullbackTests.jl") end
   @time @testset "MultiField"       begin include("MultiFieldTests.jl") end
   @time @testset "issue_142"        begin include("issue_142.jl") end
   @time @testset "ZeroMeanFESpaces" begin include("ZeroMeanFESpacesTests.jl") end
@@ -20,7 +21,6 @@ end
 if TESTCASE ∈ ("all", "seq-physics")
   @time @testset "Poisson"                     begin include("PoissonTests.jl") end
   @time @testset "PLaplacian"                  begin include("PLaplacianTests.jl") end
-  @time @testset "DivAndCurlConformingTests"   begin include("DivAndCurlConformingTests.jl") end
   @time @testset "SurfaceCouplingTests"        begin include("SurfaceCouplingTests.jl") end
   @time @testset "StokesHdivDGTests"           begin include("StokesHdivDGTests.jl") end
   @time @testset "StokesOpenBoundary"          begin include("StokesOpenBoundaryTests.jl") end