diff --git a/Project.toml b/Project.toml
index 416292f..b59f2ff 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,9 +16,9 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 ArgParse = "1"
-FillArrays = "0.8.4, 0.9, 0.10, 0.11, 0.12"
-Gridap = "0.17.18"
-GridapDistributed = "0.3"
+FillArrays = "0.8.4, 0.9, 0.10, 0.11, 0.12, 1"
+Gridap = "0.17.20"
+GridapDistributed = "0.3.1"
 MPI = "0.20"
 P4est_wrapper = "0.2.0"
 PartitionedArrays = "0.3.3"
diff --git a/README.md b/README.md
index ecedf9a..dab97de 100644
--- a/README.md
+++ b/README.md
@@ -10,7 +10,7 @@
 
 ## Purpose 
 
-The purpose of this package is to provide a `DistributedDiscreteModel` implementation (a parallel mesh data structure, see `GridapDistributed.jl` for more details) able to handle forests of quadtrees/octrees of the computational domain. To this end, it leverages the [`p4est` software library](https://p4est.github.io/) meshing engine under the hood (click [here](https://github.com/gridap/GridapP4est.jl/blob/main/test/UniformlyRefinedForestOfOctreesDiscreteModelsTests.jl) for an example).
+The purpose of this package is to provide a `DistributedDiscreteModel` implementation (a parallel mesh data structure, see `GridapDistributed.jl` for more details) able to handle forests of quadtrees/octrees of the computational domain. To this end, it leverages the [`p4est` software library](https://p4est.github.io/) meshing engine under the hood (click [here](https://github.com/gridap/GridapP4est.jl/blob/main/test/PoissonUniformlyRefinedOctreeModelsTests.jl) for an example).
 
 ## Build 
 
diff --git a/src/FESpaces.jl b/src/FESpaces.jl
index d8da386..f605120 100644
--- a/src/FESpaces.jl
+++ b/src/FESpaces.jl
@@ -14,6 +14,121 @@ function _build_constraint_coefficients_matrix_in_ref_space(Dc, reffe::Tuple{<:L
     ref_constraints = evaluate(dof_basis_h_refined, coarse_shape_funs)
 end
 
+function _generate_face_subface_ldof_to_cell_ldof(Df,Dc,reffe::Tuple{<:Lagrangian,Any,Any})
+    cell_polytope = (Dc == 2) ? QUAD : HEX
+    rr = Gridap.Adaptivity.RedRefinementRule(cell_polytope)
+    order=reffe[2][2]
+    Gridap.Adaptivity.get_face_subface_ldof_to_cell_ldof(rr, Tuple(fill(order,Dc)), Df)
+end 
+
+const _coarse_faces_to_child_ids_2D=[1 2; 3 4; 1 3; 2 4]
+const _coarse_faces_to_child_ids_3D=[1 2 3 4; 5 6 7 8; 1 2 5 6; 3 4 7 8; 1 3 5 7; 2 4 6 8; ]
+
+function _generate_unit_hypercube_model(Dc)
+    @assert Dc==2 || Dc==3
+    if (Dc==2)
+      modelH=CartesianDiscreteModel((0,1,0,1),(1,1))
+    else
+      modelH=CartesianDiscreteModel((0,1,0,1,0,1),(1,1,1))   
+    end
+    modelH
+end 
+
+function _allocate_face_subface_ldof_to_cell_ldof(num_faces, num_subfaces, num_dofs_x_face)
+    face_subface_ldof_to_cell_ldof = Vector{Vector{Vector{Int}}}(undef,num_faces)
+    for face=1:num_faces
+        face_subface_ldof_to_cell_ldof[face]=Vector{Vector{Int}}(undef,num_subfaces)
+        for subface=1:num_subfaces
+            face_subface_ldof_to_cell_ldof[face][subface]=Vector{Int}(undef,num_dofs_x_face)
+        end
+    end 
+    face_subface_ldof_to_cell_ldof
+end 
+
+function _generate_face_subface_ldof_to_cell_ldof(Df,Dc,reffe::Tuple{<:RaviartThomas,Any,Any})
+    cell_polytope = (Dc == 2) ? QUAD : HEX
+    coarse_faces_to_child_ids = (Dc == 2) ? _coarse_faces_to_child_ids_2D : 
+                                            _coarse_faces_to_child_ids_3D
+    
+    if (Df==Dc-1) # Facets    
+        basis, reffe_args, reffe_kwargs = reffe
+        cell_reffe = ReferenceFE(cell_polytope, basis, reffe_args...; reffe_kwargs...)
+
+        # TO-DO: How can modelH be created such that it is tailored to cell_polytope?
+        modelH= _generate_unit_hypercube_model(Dc)
+        modelh=refine(modelH,2)
+        RTh=TestFESpace(modelh,reffe)
+
+        num_faces = 2*Dc
+        num_subfaces = 2^(Dc-1)
+        first_face = get_offset(get_polytope(cell_reffe),Df)
+
+        face_own_dofs=get_face_own_dofs(cell_reffe)
+        num_dofs_x_face = length(face_own_dofs[first_face+1])
+        
+        face_subface_ldof_to_cell_ldof=_allocate_face_subface_ldof_to_cell_ldof(num_faces,num_subfaces,num_dofs_x_face)
+
+        RTh_cell_dof_ids=get_cell_dof_ids(RTh)
+        for coarse_face_id=1:num_faces 
+            for (subface,child_id) in enumerate(coarse_faces_to_child_ids[coarse_face_id,:])
+                @debug "coarse_face_id: $(coarse_face_id), subface: $(subface), child_id: $(child_id)"
+                cell_dof_ids=RTh_cell_dof_ids[child_id]
+                for (i,dof) in enumerate(cell_dof_ids[face_own_dofs[first_face+coarse_face_id]])
+                    @debug "i: $(i), dof: $(dof)"
+                    face_subface_ldof_to_cell_ldof[coarse_face_id][subface][i]=dof 
+                end
+            end 
+        end
+    else
+        @assert Df==1 && Dc==3 # Edges in 3D 
+        num_edges=12
+        num_subedges=2 
+        num_dofs_x_edge=0
+        face_subface_ldof_to_cell_ldof=_allocate_face_subface_ldof_to_cell_ldof(num_edges,
+                                                                                num_subedges,
+                                                                                num_dofs_x_edge)        
+    end 
+    face_subface_ldof_to_cell_ldof
+end 
+
+function _build_constraint_coefficients_matrix_in_ref_space(Dc, reffe::Tuple{<:RaviartThomas,Any,Any})
+    cell_polytope = Dc == 2 ? QUAD : HEX
+    basis, reffe_args, reffe_kwargs = reffe
+    cell_reffe = ReferenceFE(cell_polytope, basis, reffe_args...; reffe_kwargs...)
+    
+    # TO-DO: How can modelH be created such that it is tailored to cell_polytope?
+    modelH= _generate_unit_hypercube_model(Dc)
+    modelh=refine(modelH,2)
+  
+    VH=TestFESpace(modelH,cell_reffe)
+    Vh=TestFESpace(modelh,cell_reffe)
+
+    uH=get_fe_basis(VH)
+    uHh=change_domain(uH,get_triangulation(modelh),ReferenceDomain())
+    σRTh=Gridap.FESpaces.get_fe_dof_basis(Vh)
+    ref_constraints_contribs=σRTh(uHh) 
+
+    ref_constraints = Matrix{Float64}(undef,num_free_dofs(Vh),num_free_dofs(VH))
+    cell_dof_ids = get_cell_dof_ids(Vh)
+    cache_cell_dof_ids = array_cache(cell_dof_ids)
+    cache_ref_constraints_contribs = array_cache(ref_constraints_contribs) 
+    for cell=1:length(cell_dof_ids)
+       current_cell_dof_ids=getindex!(cache_cell_dof_ids,cell_dof_ids,cell)
+       current_ref_constraints_contribs=getindex!(cache_ref_constraints_contribs,ref_constraints_contribs,cell)
+       ref_constraints[current_cell_dof_ids,:]=current_ref_constraints_contribs
+    end
+    # We change the sign of the constraints coefficients here so that 
+    # the unit normal to hanging faces matches the unit normal 
+    # of the owner face. This way, we properly glue the global DoFs at 
+    # both sides of the interface of cells at different refinement level.
+    # We could have done this instead by adjusting the sign_flips for hanging
+    # faces in the FESpace constructor, although we decide it to do it here 
+    # because it is way simpler.
+    @. ref_constraints = -ref_constraints
+    ref_constraints
+end
+
+
 # To-think: might this info go to the glue? 
 # If it is required in different scenarios, I would say it may make sense
 function _generate_hanging_faces_to_cell_and_lface(num_regular_faces,
@@ -108,6 +223,34 @@ function face_lid_within_dim(::Type{Val{Dc}}, face_lid) where {Dc}
     end
 end
 
+function _restrict_face_dofs_to_face_dim(cell_reffe,Df)
+    polytope = get_polytope(cell_reffe)
+    first_face = Gridap.ReferenceFEs.get_offset(polytope,Df)+1
+    face_own_dofs=get_face_own_dofs(cell_reffe)
+    first_face_faces = Gridap.ReferenceFEs.get_faces(polytope)[first_face]
+    face_dofs_to_face_dim = face_own_dofs[first_face_faces]
+    touched=Dict{Int,Int}()
+    current=1
+    for (i,face_dofs) in enumerate(face_dofs_to_face_dim)
+        for (j,dof) in enumerate(face_dofs)
+            if haskey(touched,dof)
+                face_dofs[j]=touched[dof]
+            else
+                touched[dof]=current
+                face_dofs[j]=touched[dof]
+                current=current+1
+            end
+        end
+    end
+    face_dofs_to_face_dim
+end 
+function _num_face_own_dofs(cell_reffe,Df)
+  polytope = get_polytope(cell_reffe)
+  first_face = Gridap.ReferenceFEs.get_offset(polytope,Df)+1
+  face_own_dofs=get_face_own_dofs(cell_reffe)
+  length(face_own_dofs[first_face])
+end 
+
 function _generate_constraints!(Df,
     Dc,
     cell_faces,
@@ -510,40 +653,68 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},
         end
 
         subface_own_dofs = Vector{Vector{Vector{Int}}}(undef, Dc - 1)
-        subface_own_dofs[Dc-1] = Gridap.ReferenceFEs.get_face_own_dofs(face_reffe)
+        subface_own_dofs[Dc-1] = _restrict_face_dofs_to_face_dim(cell_reffe,Dc-1)
+        #Gridap.ReferenceFEs.get_face_own_dofs(face_reffe)
         if (Dc == 3)
-            subface_own_dofs[1] = Gridap.ReferenceFEs.get_face_own_dofs(edge_reffe)
+            subface_own_dofs[1] = _restrict_face_dofs_to_face_dim(cell_reffe,1)
+            #Gridap.ReferenceFEs.get_face_own_dofs(edge_reffe)
+        end
+        if (_num_face_own_dofs(cell_reffe,0)>0)
+            _generate_constraints!(0,
+                Dc,
+                [gridap_cell_faces[i] for i = 1:Dc],
+                num_hanging_faces[1],
+                hanging_faces_to_cell[1],
+                hanging_faces_to_lface[1],
+                hanging_faces_owner_face_dofs[1],
+                hanging_faces_glue[1],
+                face_subface_ldof_to_cell_ldof,
+                face_dofs,
+                face_own_dofs,
+                subface_own_dofs,
+                cell_dof_ids,
+                node_permutations,
+                owner_faces_pindex,
+                owner_faces_lids,
+                ref_constraints,
+                sDOF_to_dof,
+                sDOF_to_dofs,
+                sDOF_to_coeffs)
         end
-        _generate_constraints!(0,
-            Dc,
-            [gridap_cell_faces[i] for i = 1:Dc],
-            num_hanging_faces[1],
-            hanging_faces_to_cell[1],
-            hanging_faces_to_lface[1],
-            hanging_faces_owner_face_dofs[1],
-            hanging_faces_glue[1],
-            face_subface_ldof_to_cell_ldof,
-            face_dofs,
-            face_own_dofs,
-            subface_own_dofs,
-            cell_dof_ids,
-            node_permutations,
-            owner_faces_pindex,
-            owner_faces_lids,
-            ref_constraints,
-            sDOF_to_dof,
-            sDOF_to_dofs,
-            sDOF_to_coeffs)
 
         if (Dc == 3)
-            _generate_constraints!(1,
+            if (_num_face_own_dofs(cell_reffe,1)>0)
+                _generate_constraints!(1,
+                    Dc,
+                    [gridap_cell_faces[i] for i = 1:Dc],
+                    num_hanging_faces[2],
+                    hanging_faces_to_cell[2],
+                    hanging_faces_to_lface[2],
+                    hanging_faces_owner_face_dofs[2],
+                    hanging_faces_glue[2],
+                    face_subface_ldof_to_cell_ldof,
+                    face_dofs,
+                    face_own_dofs,
+                    subface_own_dofs,
+                    cell_dof_ids,
+                    node_permutations,
+                    owner_faces_pindex,
+                    owner_faces_lids,
+                    ref_constraints,
+                    sDOF_to_dof,
+                    sDOF_to_dofs,
+                    sDOF_to_coeffs)
+            end
+        end
+        if (_num_face_own_dofs(cell_reffe,Dc-1)>0)
+            _generate_constraints!(Dc - 1,
                 Dc,
                 [gridap_cell_faces[i] for i = 1:Dc],
-                num_hanging_faces[2],
-                hanging_faces_to_cell[2],
-                hanging_faces_to_lface[2],
-                hanging_faces_owner_face_dofs[2],
-                hanging_faces_glue[2],
+                num_hanging_faces[Dc],
+                hanging_faces_to_cell[Dc],
+                hanging_faces_to_lface[Dc],
+                hanging_faces_owner_face_dofs[Dc],
+                hanging_faces_glue[Dc],
                 face_subface_ldof_to_cell_ldof,
                 face_dofs,
                 face_own_dofs,
@@ -557,26 +728,6 @@ function generate_constraints(dmodel::OctreeDistributedDiscreteModel{Dc},
                 sDOF_to_dofs,
                 sDOF_to_coeffs)
         end
-        _generate_constraints!(Dc - 1,
-            Dc,
-            [gridap_cell_faces[i] for i = 1:Dc],
-            num_hanging_faces[Dc],
-            hanging_faces_to_cell[Dc],
-            hanging_faces_to_lface[Dc],
-            hanging_faces_owner_face_dofs[Dc],
-            hanging_faces_glue[Dc],
-            face_subface_ldof_to_cell_ldof,
-            face_dofs,
-            face_own_dofs,
-            subface_own_dofs,
-            cell_dof_ids,
-            node_permutations,
-            owner_faces_pindex,
-            owner_faces_lids,
-            ref_constraints,
-            sDOF_to_dof,
-            sDOF_to_dofs,
-            sDOF_to_coeffs)
         sDOF_to_dof, Gridap.Arrays.Table(sDOF_to_dofs), Gridap.Arrays.Table(sDOF_to_coeffs)
     end |> tuple_of_arrays
 end
@@ -619,54 +770,74 @@ function fe_space_with_linear_constraints_cell_dof_ids(Uc::FESpaceWithLinearCons
     Gridap.Arrays.Table(Uc_cell_dof_ids_data, U_cell_dof_ids.ptrs)
 end
 
-# Generates a new DistributedSingleFieldFESpace composed 
-# by local FE spaces with linear multipoint constraints added
-function Gridap.FESpaces.FESpace(model::OctreeDistributedDiscreteModel{Dc}, 
-                                 reffe::Tuple{Gridap.ReferenceFEs.ReferenceFEName,Any,Any}; 
-                                 kwargs...) where {Dc}
-    order = reffe[2][2]
-    spaces_wo_constraints = map(local_views(model)) do m
-        FESpace(m, reffe; kwargs...)
-    end
-    ref_constraints = _build_constraint_coefficients_matrix_in_ref_space(Dc, reffe)
-    cell_polytope = (Dc == 2) ? QUAD : HEX
-    rr = Gridap.Adaptivity.RedRefinementRule(cell_polytope)
-    face_subface_ldof_to_cell_ldof = Vector{Vector{Vector{Vector{Int32}}}}(undef, Dc-1)
-    face_subface_ldof_to_cell_ldof[Dc-1] =
-        Gridap.Adaptivity.get_face_subface_ldof_to_cell_ldof(rr, Tuple(fill(order,Dc)), Dc-1)
-    if (Dc == 3)
-        face_subface_ldof_to_cell_ldof[1] =
-            Gridap.Adaptivity.get_face_subface_ldof_to_cell_ldof(rr, Tuple(fill(order,Dc)), 1)
-    end
-    sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs =
-        generate_constraints(model, spaces_wo_constraints, reffe, ref_constraints, face_subface_ldof_to_cell_ldof)
-
-    spaces_w_constraints = map(spaces_wo_constraints,
-        sDOF_to_dof,
-        sDOF_to_dofs,
-        sDOF_to_coeffs) do V, sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs
-        @debug "fe_space_wo_constraints_cell_dof_ids=$(get_cell_dof_ids(V))"
-        Vc = FESpaceWithLinearConstraints(sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, V)
-    end
-
-    local_cell_dof_ids = map(spaces_w_constraints,sDOF_to_dof) do Vc,sDOF_to_dof
+function _add_constraints(model::GridapDistributed.DistributedDiscreteModel{Dc},
+                          reffe,
+                          spaces_wo_constraints;
+                          conformity=nothing,
+                          kwargs...) where {Dc}
+    if (conformity!=nothing && conformity!=:L2)
+      ref_constraints = _build_constraint_coefficients_matrix_in_ref_space(Dc, reffe)
+
+      face_subface_ldof_to_cell_ldof = Vector{Vector{Vector{Vector{Int32}}}}(undef, Dc-1)
+      face_subface_ldof_to_cell_ldof[Dc-1] = _generate_face_subface_ldof_to_cell_ldof(Dc-1, Dc, reffe)
+      if (Dc == 3)
+       face_subface_ldof_to_cell_ldof[1] =
+         _generate_face_subface_ldof_to_cell_ldof(1, Dc, reffe)
+      end
+      sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs =
+          generate_constraints(model, spaces_wo_constraints, reffe, ref_constraints, face_subface_ldof_to_cell_ldof)
+      spaces_w_constraints = map(spaces_wo_constraints,
+         sDOF_to_dof,
+         sDOF_to_dofs,
+         sDOF_to_coeffs) do V, sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs
+         @debug "[$(MPI.Comm_rank(MPI.COMM_WORLD))]: fe_space_wo_constraints_cell_dof_ids=$(get_cell_dof_ids(V))"
+         Vc = FESpaceWithLinearConstraints(sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, V)
+      end
+      local_cell_dof_ids = map(spaces_w_constraints,sDOF_to_dof) do Vc,sDOF_to_dof
         result = fe_space_with_linear_constraints_cell_dof_ids(Vc,sDOF_to_dof)
-        @debug "fe_space_with_linear_constraints_cell_dof_ids=$(result)"
+        @debug "[$(MPI.Comm_rank(MPI.COMM_WORLD))]: fe_space_with_linear_constraints_cell_dof_ids=$(result)"
         result
+      end
+    else 
+      spaces_w_constraints=spaces_wo_constraints
+      local_cell_dof_ids=map(get_cell_dof_ids,spaces_w_constraints)
     end
     nldofs = map(num_free_dofs,spaces_w_constraints)
     cell_gids = get_cell_gids(model)
     gids = GridapDistributed.generate_gids(cell_gids,local_cell_dof_ids,nldofs)
     map(partition(gids)) do indices 
         @debug "[$(part_id(indices))]: l2g_cell_gids=$(local_to_global(indices))"
-        @debug "[$(part_id(indices))]: l2o_cell_gids=$(local_to_owner(indices))"
+        @debug "[$(part_id(indices))]: l2o_owner=$(local_to_owner(indices))"
     end 
     vector_type = GridapDistributed._find_vector_type(spaces_w_constraints,gids)
     GridapDistributed.DistributedSingleFieldFESpace(spaces_w_constraints,gids,vector_type)
 end
 
-function Gridap.FESpaces.FESpace(model::OctreeDistributedDiscreteModel, 
+
+# Generates a new DistributedSingleFieldFESpace composed 
+# by local FE spaces with linear multipoint constraints added
+function Gridap.FESpaces.FESpace(model::OctreeDistributedDiscreteModel{Dc}, 
+                                 reffe::Tuple{Gridap.ReferenceFEs.Lagrangian,Any,Any}; 
+                                 kwargs...) where {Dc}
+    spaces_wo_constraints = map(local_views(model)) do m
+        FESpace(m, reffe; kwargs...)
+    end
+    _add_constraints(model,reffe,spaces_wo_constraints;kwargs...)
+end
+
+function Gridap.FESpaces.FESpace(model::OctreeDistributedDiscreteModel{Dc}, 
                                  reffe::Tuple{Gridap.ReferenceFEs.RaviartThomas,Any,Any}; 
-                                 kwargs...)
-    FESpace(model.dmodel,reffe; kwargs...)
+                                 conformity=nothing,kwargs...) where {Dc}
+
+    cell_reffes = map(local_views(model.dmodel)) do m
+        basis,reffe_args,reffe_kwargs = reffe
+        cell_reffe = ReferenceFE(m,basis,reffe_args...;reffe_kwargs...)
+    end
+    sign_flips=GridapDistributed._generate_sign_flips(model.dmodel,cell_reffes)
+    spaces_wo_constraints = map(local_views(model.dmodel),sign_flips,cell_reffes) do m,sign_flip,cell_reffe
+       conf = Conformity(Gridap.Fields.testitem(cell_reffe),conformity)
+       cell_fe = CellFE(m,cell_reffe,conf,sign_flip)
+       FESpace(m, cell_fe; kwargs...)
+    end
+    _add_constraints(model,reffe,spaces_wo_constraints;conformity=conformity,kwargs...)
 end
diff --git a/src/GridapFixes.jl b/src/GridapFixes.jl
index 10d49b8..41dc1fb 100644
--- a/src/GridapFixes.jl
+++ b/src/GridapFixes.jl
@@ -4,4 +4,22 @@ function Base.map(::typeof(Gridap.Arrays.testitem),
   a1=Vector{eltype(eltype(a[1]))}(undef,size(a2,1))
   a1.=zero(Gridap.Arrays.testitem(a1))
   (a1,a2)
+end
+
+# Required to transfer fine-grid VECTOR-VALUED fields into coarse-grid 
+function Gridap.Adaptivity.FineToCoarseField(fine_fields::AbstractArray{<:Gridap.Fields.Field},
+                                              rrule::Gridap.Adaptivity.RefinementRule,
+                                              child_ids::AbstractArray{<:Integer})
+  
+  grid=Gridap.Adaptivity.get_ref_grid(rrule)
+  D=num_cell_dims(grid)
+  x=zero(Point{D,Float64})
+  ffx=lazy_map(evaluate,fine_fields,Fill([x],length(fine_fields)))
+  ffx=ffx[1]
+  fields = Vector{Gridap.Fields.Field}(undef,Gridap.Adaptivity.num_subcells(rrule))
+  fields = fill!(fields,Gridap.Fields.ConstantField(zero(eltype(ffx))))
+  for (k,id) in enumerate(child_ids)
+    fields[id] = fine_fields[k]
+  end
+  return Gridap.Adaptivity.FineToCoarseField(fields,rrule)
 end
\ No newline at end of file
diff --git a/test/CoarseDiscreteModelsTools.jl b/test/CoarseDiscreteModelsTools.jl
new file mode 100644
index 0000000..c5a636f
--- /dev/null
+++ b/test/CoarseDiscreteModelsTools.jl
@@ -0,0 +1,120 @@
+  function setup_model(::Type{Val{3}}, perm)
+  #               5 +--------+ 7 
+  #                /        /|
+  #               /        / |
+  #            6 +--------+  |
+  #              |        |  |
+  #              |  1     |  + 3 
+  #              |        | /
+  #              |        |/
+  #            2 +--------+ 4
+
+  #     6  +--------+ 8 
+  #       /        /|
+  #      /        / |
+  #  11 +--------+  |
+  #     |        |  |
+  #     |  2     |  + 4
+  #     |        | /
+  #     |        |/
+  #   9 +--------+ 10
+    ptr  = [ 1, 9, 17 ]
+    if (perm==1)
+      data = [ 1,2,3,4,5,6,7,8, 2,9,4,10,6,11,8,12 ]
+    elseif (perm==2)
+      data = [ 1,2,3,4,5,6,7,8, 10,12,4,8,9,11,2,6 ]
+    elseif (perm==3)
+      data = [ 1,2,3,4,5,6,7,8, 12,11,8,6,10,9,4,2 ]
+    elseif (perm==4) 
+      data = [ 1,2,3,4,5,6,7,8, 11,9,6,2,12,10,8,4 ]
+    end  
+    cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
+    node_coordinates = Vector{Point{3,Float64}}(undef,12)
+    node_coordinates[1]=Point{3,Float64}(0.0,0.0,0.0)
+    node_coordinates[2]=Point{3,Float64}(1.0,0.0,0.0)
+    node_coordinates[3]=Point{3,Float64}(0.0,1.0,0.0)
+    node_coordinates[4]=Point{3,Float64}(1.0,1.0,0.0)
+    node_coordinates[5]=Point{3,Float64}(0.0,0.0,1.0)
+    node_coordinates[6]=Point{3,Float64}(1.0,0.0,1.0)
+    node_coordinates[7]=Point{3,Float64}(0.0,1.0,1.0)
+    node_coordinates[8]=Point{3,Float64}(1.0,1.0,1.0)
+    node_coordinates[9]=Point{3,Float64}(2.0,0.0,0.0)
+    node_coordinates[10]=Point{3,Float64}(2.0,1.0,0.0)
+    node_coordinates[11]=Point{3,Float64}(2.0,0.0,1.0)
+    node_coordinates[12]=Point{3,Float64}(2.0,1.0,1.0)
+
+    polytope=HEX
+    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())
+    m=Gridap.Geometry.UnstructuredDiscreteModel(grid)
+    labels = get_face_labeling(m)
+    labels.d_to_dface_to_entity[1].=2
+    if (perm==1 || perm==2)
+      labels.d_to_dface_to_entity[2].=2
+      labels.d_to_dface_to_entity[3].=[2,2,2,2,2,1,2,2,2,2,2]
+    elseif (perm==3 || perm==4)
+      labels.d_to_dface_to_entity[2].=2
+      labels.d_to_dface_to_entity[3].=[2,2,2,2,2,1,2,2,2,2,2]
+    end
+    labels.d_to_dface_to_entity[4].=1  
+    add_tag!(labels,"boundary",[2])
+    add_tag!(labels,"interior",[1])
+    m
+  end
+
+  function setup_model(::Type{Val{2}}, perm)
+    @assert perm ∈ (1,2,3,4)
+    #
+    #  3-------4-------6
+    #  |       |       |
+    #  |       |       |
+    #  |       |       |
+    #  1-------2-------5
+    #
+        ptr  = [ 1, 5, 9 ]
+        if (perm==1)
+          data = [ 1,2,3,4, 2,5,4,6 ]
+        elseif (perm==2)
+          data = [ 1,2,3,4, 6,4,5,2 ]
+        elseif (perm==3)
+          data = [ 4,3,2,1, 2,5,4,6 ]
+        elseif (perm==4) 
+          data = [ 4,3,2,1, 6,4,5,2 ]
+        end  
+        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}(1.0,0.0)
+        node_coordinates[3]=Point{2,Float64}(0.0,1.0)
+        node_coordinates[4]=Point{2,Float64}(1.0,1.0)
+        node_coordinates[5]=Point{2,Float64}(2.0,0.0)
+        node_coordinates[6]=Point{2,Float64}(2.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())
+        m=Gridap.Geometry.UnstructuredDiscreteModel(grid)
+        labels = get_face_labeling(m)
+        labels.d_to_dface_to_entity[1].=2
+        if (perm==1 || perm==2)
+          labels.d_to_dface_to_entity[2].=[2,2,2,1,2,2,2]
+        elseif (perm==3 || perm==4)
+          labels.d_to_dface_to_entity[2].=[2,2,1,2,2,2,2] 
+        end
+        labels.d_to_dface_to_entity[3].=1
+        add_tag!(labels,"boundary",[2])
+        add_tag!(labels,"interior",[1])
+        m
+  end
\ No newline at end of file
diff --git a/test/DarcyNonConformingOctreeModelsTests.jl b/test/DarcyNonConformingOctreeModelsTests.jl
new file mode 100644
index 0000000..bfe0c94
--- /dev/null
+++ b/test/DarcyNonConformingOctreeModelsTests.jl
@@ -0,0 +1,282 @@
+module DarcyNonConformingOctreeModelsTests
+  using P4est_wrapper
+  using GridapP4est
+  using Gridap
+  using PartitionedArrays
+  using GridapDistributed
+  using MPI
+  using Gridap.FESpaces
+  using FillArrays
+  using Logging
+  using Test
+
+  function test_transfer_ops_and_redistribute(ranks,
+                                              dmodel::GridapDistributed.DistributedDiscreteModel{Dc},
+                                              order) where Dc
+    ref_coarse_flags=map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices
+        flags=zeros(Cint,length(indices))
+        flags.=nothing_flag
+        
+        flags[1]=refine_flag
+        flags[own_length(indices)]=refine_flag
+
+        # To create some unbalance
+        if (rank%2==0 && own_length(indices)>1)
+            flags[div(own_length(indices),2)]=refine_flag
+        end
+        flags
+    end 
+    fmodel,glue=adapt(dmodel,ref_coarse_flags);
+    
+    # Solve coarse
+    xH,XH=solve_darcy(dmodel,order)
+    check_error_darcy(dmodel,order,xH)
+    uH,pH=xH
+    UH,PH=XH
+
+    # Solve fine
+    xh,Xh=solve_darcy(fmodel,order)
+    check_error_darcy(fmodel,order,xh)
+    uh,ph=xh
+    Uh,Ph=Xh
+
+    Ωh = Triangulation(fmodel)
+    degree = 2*(order+1)
+    dΩh = Measure(Ωh,degree)
+
+    # prolongation via interpolation
+    uHh=interpolate(uH,Uh)   
+    e = uh - uHh
+    el2 = sqrt(sum( ∫( e⋅e )*dΩh ))
+    tol=1e-6
+    println("[INTERPOLATION] el2 < tol: $(el2) < $(tol)")
+    @assert el2 < tol
+
+    # prolongation via L2-projection 
+    # Coarse FEFunction -> Fine FEFunction, by projection
+    ahp(u,v)  = ∫(v⋅u)*dΩh
+    lhp(v)    = ∫(v⋅uH)*dΩh
+    oph      = AffineFEOperator(ahp,lhp,Uh,Uh)
+    uHh      = solve(oph)
+    e = uh - uHh
+    el2 = sqrt(sum( ∫( e⋅e )*dΩh ))
+    println("[L2 PROJECTION] el2 < tol: $(el2) < $(tol)")
+    @assert el2 < tol
+
+    # restriction via interpolation
+    uhH=interpolate(uh,UH) 
+    e = uH - uhH
+    el2 = sqrt(sum( ∫( e⋅e )*dΩh ))
+    println("[INTERPOLATION] el2 < tol: $(el2) < $(tol)")
+    @assert el2 < tol
+
+    # restriction via L2-projection
+    ΩH = Triangulation(dmodel)
+    degree = 2*(order+1)
+    dΩH = Measure(ΩH,degree)
+    
+    dΩhH = Measure(ΩH,Ωh,2*order)
+    aHp(u,v) = ∫(v⋅u)*dΩH
+    lHp(v)   = ∫(v⋅uh)*dΩhH
+    oph     = AffineFEOperator(aHp,lHp,UH,UH)
+    uhH     = solve(oph)
+    e       = uH - uhH
+    el2     = sqrt(sum( ∫( e⋅e )*dΩH ))
+
+    fmodel_red, red_glue=GridapDistributed.redistribute(fmodel);
+    xh_red,Xh_red=solve_darcy(fmodel_red,order)
+    check_error_darcy(fmodel_red,order,xh_red)
+    uhred,phred=xh_red
+    Uh_red,Ph_red=Xh_red
+
+
+    trian = Triangulation(fmodel_red)
+    degree = 2*(order+1)
+    dΩhred = Measure(trian,degree)
+
+    u_ex, p_ex, f_ex=get_analytical_functions(Dc)
+
+    uhred2 = GridapDistributed.redistribute_fe_function(uh,Uh_red,fmodel_red,red_glue)
+    
+    e = u_ex - uhred2
+    el2 = sqrt(sum( ∫( e⋅e )*dΩhred ))
+    println("[REDISTRIBUTE SOLUTION] el2 < tol: $(el2) < $(tol)")
+    @assert el2 < tol
+
+    fmodel_red
+  end
+
+  function test_refine_and_coarsen_at_once(ranks,
+            dmodel::OctreeDistributedDiscreteModel{Dc},
+            order) where Dc
+    degree = 2*order+1
+    ref_coarse_flags=map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices
+        flags=zeros(Cint,length(indices))
+        flags.=nothing_flag        
+        if (rank==1)
+           flags[1:2^Dc].=coarsen_flag
+        else 
+           flags[own_length(indices)]=refine_flag
+        end 
+        flags
+    end
+    fmodel,glue=adapt(dmodel,ref_coarse_flags);
+
+    # Solve coarse
+    xH,XH=solve_darcy(dmodel,order)
+    check_error_darcy(dmodel,order,xH)
+    uH,pH=xH
+
+    # Solve fine
+    xh,Xh=solve_darcy(fmodel,order)
+    check_error_darcy(fmodel,order,xh)
+    uh,ph=xh
+
+    # prolongation via interpolation
+    Uh,Ph=Xh
+    uHh=interpolate(uH,Uh)
+    e = uh - uHh
+
+    trian = Triangulation(fmodel)
+    degree = 2*(order+1)
+    dΩh = Measure(trian,degree)
+
+    el2 = sqrt(sum( ∫( e⋅e )*dΩh ))
+    tol=1e-7
+    @assert el2 < tol
+  end
+
+  function test_2d(ranks,order)
+    coarse_model=CartesianDiscreteModel((0,1,0,1),(1,1))
+    dmodel=OctreeDistributedDiscreteModel(ranks,coarse_model,2)
+    test_refine_and_coarsen_at_once(ranks,dmodel,order)
+    rdmodel=dmodel
+    for i=1:5
+     rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order)
+    end
+  end 
+
+  function test_3d(ranks,order)
+    coarse_model=CartesianDiscreteModel((0,1,0,1,0,1),(1,1,1))
+    dmodel=OctreeDistributedDiscreteModel(ranks,coarse_model,2)
+    test_refine_and_coarsen_at_once(ranks,dmodel,order)
+    rdmodel=dmodel
+    for i=1:5
+      rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order)
+    end
+  end
+
+  u_ex_2D(x) = VectorValue(2*x[1],x[1]+x[2])
+  p_ex_2D(x) = x[1]-x[2]
+  f_ex_2D(x) = u_ex_2D(x) + ∇(p_ex_2D)(x)
+  u_ex_3D(x) = VectorValue(x[1],x[2],x[3])
+  p_ex_3D(x) = x[2]
+  f_ex_3D(x) = u_ex_3D(x) + ∇(p_ex_3D)(x)
+  
+  function get_analytical_functions(Dc)
+    if Dc==2
+      return u_ex_2D, p_ex_2D, f_ex_2D
+    else
+      @assert Dc==3
+      return u_ex_3D, p_ex_3D, f_ex_3D
+    end
+  end
+
+  include("CoarseDiscreteModelsTools.jl")
+
+  function solve_darcy(model::GridapDistributed.DistributedDiscreteModel{Dc},order) where {Dc}
+    if (Dc==2)
+      dirichlet_tags=[5,6]
+      neumanntags = [7,8]
+    else
+      @assert Dc==3
+      dirichlet_tags=[21,22,23]
+      neumanntags=[24,25,26]
+    end 
+
+    u_ex, p_ex, f_ex=get_analytical_functions(Dc)
+
+    V = FESpace(model,
+                ReferenceFE(raviart_thomas,Float64,order),
+                conformity=:Hdiv,
+                dirichlet_tags=dirichlet_tags)
+    
+    Q = FESpace(model,
+                ReferenceFE(lagrangian,Float64,order); 
+                conformity=:L2)
+    
+    U = TrialFESpace(V,u_ex)
+    P = TrialFESpace(Q)
+    
+    Y = MultiFieldFESpace([V, Q])
+    X = MultiFieldFESpace([U, P])
+    
+    trian = Triangulation(model)
+    degree = 2*(order+1)
+    dΩ = Measure(trian,degree)
+    
+    btrian = Boundary(model,tags=neumanntags)
+    degree = 2*(order+1)
+    dΓ = Measure(btrian,degree)
+    nb = get_normal_vector(btrian)
+    
+    a((u, p),(v, q)) = ∫(u⋅v)dΩ +∫(q*(∇⋅u))dΩ-∫((∇⋅v)*p)dΩ
+    b(( v, q)) = ∫( v⋅f_ex + q*(∇⋅u_ex))dΩ - ∫((v⋅nb)*p_ex )dΓ
+
+    op = AffineFEOperator(a,b,X,Y)
+    xh = solve(op)
+    xh,X
+  end 
+
+  function check_error_darcy(model::GridapDistributed.DistributedDiscreteModel{Dc},order,xh) where {Dc}
+    trian = Triangulation(model)
+    degree = 2*(order+1)
+    dΩ = Measure(trian,degree)
+
+    uh, ph = xh
+
+    u_ex, p_ex, f_ex = get_analytical_functions(Dc)
+    
+    eu = u_ex - uh
+    ep = p_ex - ph
+
+    l2(v) = sqrt(sum(∫(v⋅v)*dΩ))
+    h1(v) = sqrt(sum(∫(v*v + ∇(v)⋅∇(v))*dΩ))
+    
+    eu_l2 = l2(eu)
+    ep_l2 = l2(ep)
+    ep_h1 = h1(ep)
+    
+    tol = 1.0e-9
+    @test eu_l2 < tol
+    @test ep_l2 < tol
+    @test ep_h1 < tol
+  end 
+
+  function driver(ranks,coarse_model,order)
+    model=OctreeDistributedDiscreteModel(ranks,coarse_model,1)
+    ref_coarse_flags=map(ranks,partition(get_cell_gids(model.dmodel))) do rank,indices
+        flags=zeros(Cint,length(indices))
+        flags.=nothing_flag        
+        #flags[1:2^2].=coarsen_flag
+        flags[own_length(indices)]=refine_flag 
+        flags
+    end
+    fmodel,glue=adapt(model,ref_coarse_flags)
+    xh,Xh = solve_darcy(fmodel,order)
+    check_error_darcy(fmodel,order,xh)
+  end 
+
+  function run(distribute)
+    # debug_logger = ConsoleLogger(stderr, Logging.Debug)
+    # global_logger(debug_logger); # Enable the debug logger globally
+    np    = MPI.Comm_size(MPI.COMM_WORLD)
+    ranks = distribute(LinearIndices((np,)))
+
+    order=1
+    test_2d(ranks,order)
+
+    order=1
+    test_3d(ranks,order)
+  end
+end
\ No newline at end of file
diff --git a/test/NonConformingOctreeDistributedDiscreteModelsTests.jl b/test/PoissonNonConformingOctreeModelsTests.jl
similarity index 60%
rename from test/NonConformingOctreeDistributedDiscreteModelsTests.jl
rename to test/PoissonNonConformingOctreeModelsTests.jl
index db7dd1f..9c7a9eb 100644
--- a/test/NonConformingOctreeDistributedDiscreteModelsTests.jl
+++ b/test/PoissonNonConformingOctreeModelsTests.jl
@@ -1,4 +1,4 @@
-module NonConformingOctreeDistributedDiscreteModelsTests
+module PoissonNonConformingOctreeModelsTests
   using P4est_wrapper
   using GridapP4est
   using Gridap
@@ -9,127 +9,7 @@ module NonConformingOctreeDistributedDiscreteModelsTests
   using FillArrays
   using Logging
 
-  function setup_model(::Type{Val{3}}, perm)
-  #               5 +--------+ 7 
-  #                /        /|
-  #               /        / |
-  #            6 +--------+  |
-  #              |        |  |
-  #              |  1     |  + 3 
-  #              |        | /
-  #              |        |/
-  #            2 +--------+ 4
-
-  #     6  +--------+ 8 
-  #       /        /|
-  #      /        / |
-  #  11 +--------+  |
-  #     |        |  |
-  #     |  2     |  + 4
-  #     |        | /
-  #     |        |/
-  #   9 +--------+ 10
-    ptr  = [ 1, 9, 17 ]
-    if (perm==1)
-      data = [ 1,2,3,4,5,6,7,8, 2,9,4,10,6,11,8,12 ]
-    elseif (perm==2)
-      data = [ 1,2,3,4,5,6,7,8, 10,12,4,8,9,11,2,6 ]
-    elseif (perm==3)
-      data = [ 1,2,3,4,5,6,7,8, 12,11,8,6,10,9,4,2 ]
-    elseif (perm==4) 
-      data = [ 1,2,3,4,5,6,7,8, 11,9,6,2,12,10,8,4 ]
-    end  
-    cell_vertex_lids = Gridap.Arrays.Table(data,ptr)
-    node_coordinates = Vector{Point{3,Float64}}(undef,12)
-    node_coordinates[1]=Point{3,Float64}(0.0,0.0,0.0)
-    node_coordinates[2]=Point{3,Float64}(1.0,0.0,0.0)
-    node_coordinates[3]=Point{3,Float64}(0.0,1.0,0.0)
-    node_coordinates[4]=Point{3,Float64}(1.0,1.0,0.0)
-    node_coordinates[5]=Point{3,Float64}(0.0,0.0,1.0)
-    node_coordinates[6]=Point{3,Float64}(1.0,0.0,1.0)
-    node_coordinates[7]=Point{3,Float64}(0.0,1.0,1.0)
-    node_coordinates[8]=Point{3,Float64}(1.0,1.0,1.0)
-    node_coordinates[9]=Point{3,Float64}(2.0,0.0,0.0)
-    node_coordinates[10]=Point{3,Float64}(2.0,1.0,0.0)
-    node_coordinates[11]=Point{3,Float64}(2.0,0.0,1.0)
-    node_coordinates[12]=Point{3,Float64}(2.0,1.0,1.0)
-
-    polytope=HEX
-    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())
-    m=Gridap.Geometry.UnstructuredDiscreteModel(grid)
-    labels = get_face_labeling(m)
-    labels.d_to_dface_to_entity[1].=2
-    if (perm==1 || perm==2)
-      labels.d_to_dface_to_entity[2].=2
-      labels.d_to_dface_to_entity[3].=[2,2,2,2,2,1,2,2,2,2,2]
-    elseif (perm==3 || perm==4)
-      labels.d_to_dface_to_entity[2].=2
-      labels.d_to_dface_to_entity[3].=[2,2,2,2,2,1,2,2,2,2,2]
-    end
-    labels.d_to_dface_to_entity[4].=1  
-    add_tag!(labels,"boundary",[2])
-    add_tag!(labels,"interior",[1])
-    m
-  end
-
-  function setup_model(::Type{Val{2}}, perm)
-    @assert perm ∈ (1,2,3,4)
-    #
-    #  3-------4-------6
-    #  |       |       |
-    #  |       |       |
-    #  |       |       |
-    #  1-------2-------5
-    #
-        ptr  = [ 1, 5, 9 ]
-        if (perm==1)
-          data = [ 1,2,3,4, 2,5,4,6 ]
-        elseif (perm==2)
-          data = [ 1,2,3,4, 6,4,5,2 ]
-        elseif (perm==3)
-          data = [ 4,3,2,1, 2,5,4,6 ]
-        elseif (perm==4) 
-          data = [ 4,3,2,1, 6,4,5,2 ]
-        end  
-        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}(1.0,0.0)
-        node_coordinates[3]=Point{2,Float64}(0.0,1.0)
-        node_coordinates[4]=Point{2,Float64}(1.0,1.0)
-        node_coordinates[5]=Point{2,Float64}(2.0,0.0)
-        node_coordinates[6]=Point{2,Float64}(2.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())
-        m=Gridap.Geometry.UnstructuredDiscreteModel(grid)
-        labels = get_face_labeling(m)
-        labels.d_to_dface_to_entity[1].=2
-        if (perm==1 || perm==2)
-          labels.d_to_dface_to_entity[2].=[2,2,2,1,2,2,2]
-        elseif (perm==3 || perm==4)
-          labels.d_to_dface_to_entity[2].=[2,2,1,2,2,2,2] 
-        end
-        labels.d_to_dface_to_entity[3].=1
-        add_tag!(labels,"boundary",[2])
-        add_tag!(labels,"interior",[1])
-        m
-  end
-
+  include("CoarseDiscreteModelsTools.jl")
 
   function test_transfer_ops_and_redistribute(ranks,dmodel,order)
     # Define manufactured functions
diff --git a/test/UniformlyRefinedForestOfOctreesDiscreteModelsTests.jl b/test/PoissonUniformlyRefinedOctreeModelsTests.jl
similarity index 97%
rename from test/UniformlyRefinedForestOfOctreesDiscreteModelsTests.jl
rename to test/PoissonUniformlyRefinedOctreeModelsTests.jl
index 3bb4099..6010449 100644
--- a/test/UniformlyRefinedForestOfOctreesDiscreteModelsTests.jl
+++ b/test/PoissonUniformlyRefinedOctreeModelsTests.jl
@@ -1,4 +1,4 @@
-module UniformlyRefinedForestOfOctreesDiscreteModelsTests
+module PoissonUniformlyRefinedOctreeModelsTests
   using MPI
   using Gridap
   using Gridap.Algebra
diff --git a/test/mpi/NonConformingOctreeDistributedDiscreteModelsTests.jl b/test/mpi/DarcyNonConformingOctreeModelsTests.jl
similarity index 52%
rename from test/mpi/NonConformingOctreeDistributedDiscreteModelsTests.jl
rename to test/mpi/DarcyNonConformingOctreeModelsTests.jl
index f420dbf..5df77c2 100644
--- a/test/mpi/NonConformingOctreeDistributedDiscreteModelsTests.jl
+++ b/test/mpi/DarcyNonConformingOctreeModelsTests.jl
@@ -2,8 +2,8 @@
 using MPI
 using PartitionedArrays
 
-include("../NonConformingOctreeDistributedDiscreteModelsTests.jl")
-import .NonConformingOctreeDistributedDiscreteModelsTests as TestModule
+include("../DarcyNonConformingOctreeModelsTests.jl")
+import .DarcyNonConformingOctreeModelsTests as TestModule
 
 if !MPI.Initialized()
   MPI.Init()
diff --git a/test/mpi/PoissonNonConformingOctreeModelsTests.jl b/test/mpi/PoissonNonConformingOctreeModelsTests.jl
new file mode 100644
index 0000000..7332a89
--- /dev/null
+++ b/test/mpi/PoissonNonConformingOctreeModelsTests.jl
@@ -0,0 +1,17 @@
+
+using MPI
+using PartitionedArrays
+
+include("../PoissonNonConformingOctreeModelsTests.jl")
+import .PoissonNonConformingOctreeModelsTests as TestModule
+
+if !MPI.Initialized()
+  MPI.Init()
+end
+
+with_mpi() do distribute 
+  TestModule.run(distribute)
+end 
+
+MPI.Finalize()
+
diff --git a/test/mpi/UniformlyRefinedForestOfOctreesDiscreteModelsTests.jl b/test/mpi/PoissonUniformlyRefinedOctreeModelsTests.jl
similarity index 71%
rename from test/mpi/UniformlyRefinedForestOfOctreesDiscreteModelsTests.jl
rename to test/mpi/PoissonUniformlyRefinedOctreeModelsTests.jl
index 17bab82..74b1cd5 100644
--- a/test/mpi/UniformlyRefinedForestOfOctreesDiscreteModelsTests.jl
+++ b/test/mpi/PoissonUniformlyRefinedOctreeModelsTests.jl
@@ -2,8 +2,8 @@
 using MPI
 using PartitionedArrays
 
-include("../UniformlyRefinedForestOfOctreesDiscreteModelsTests.jl")
-import .UniformlyRefinedForestOfOctreesDiscreteModelsTests as TestModule
+include("../PoissonUniformlyRefinedOctreeModelsTests.jl")
+import .PoissonUniformlyRefinedOctreeModelsTests as TestModule
 
 if !MPI.Initialized()
   MPI.Init()
diff --git a/test/runtests.jl b/test/runtests.jl
index 31f56ab..31acee4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -33,15 +33,16 @@ function run_tests(testdir)
     testfiles = sort(filter(istest, readdir(testdir)))
     @time @testset "$f" for f in testfiles
       MPI.mpiexec() do cmd
-        if f in ["UniformlyRefinedForestOfOctreesDiscreteModelsTests.jl"]
+        if f in ["PoissonUniformlyRefinedOctreeModelsTests.jl"]
           np = [4]
           extra_args = "-s 2 2 -r 2"
         elseif f in ["OctreeDistributedDiscreteModelsTests.jl",
                      "OctreeDistributedDiscreteModelsNoEnvTests.jl"]
           np = [4]
           extra_args = ""
-        elseif f in ["NonConformingOctreeDistributedDiscreteModelsTests.jl",
-                     "AdaptivityFlagsMarkingStrategiesTests.jl"]
+        elseif f in ["PoissonNonConformingOctreeModelsTests.jl",
+                     "AdaptivityFlagsMarkingStrategiesTests.jl",
+                     "DarcyNonConformingOctreeModelsTests.jl"]
           np = [1,2,4]
           extra_args = ""
         else