diff --git a/include/entity.h b/include/entity.h
index 8c052f1fcdbbfc1ee52796baf23e0e41c4b7f1da..45447fa764dc4c5801229d9ec629864c2ab1b2e3 100644
--- a/include/entity.h
+++ b/include/entity.h
@@ -27,10 +27,18 @@ class entity
     int     parent_dim;
     int     parent_tag;
 
-    size_t                              m_dof_base;
-    size_t                              m_flux_base;
-    size_t                              m_index_base;
-    size_t                              m_entity_number;
+    /* Each entity has model offsets, which are offsets relative to the
+     * computational unit, and world offsets, which are offsets needed
+     * to recover the whole-world vectors of DOFs. */
+    size_t                              m_dof_base_model;
+    size_t                              m_flux_base_model;
+    size_t                              m_index_base_model;
+    size_t                              m_entity_number_model;
+
+    size_t                              m_dof_base_world;
+    size_t                              m_flux_base_world;
+    size_t                              m_index_base_world;
+    size_t                              m_entity_number_world;
 
     entity_ordering                     cur_elem_ordering;
 
@@ -82,11 +90,14 @@ public:
     const physical_element&     cell(size_t) const;
     const reference_element&    cell_refelem(size_t) const;
     const reference_element&    cell_refelem(const physical_element&) const;
-    size_t                      cell_global_index(size_t) const;
+    size_t                      cell_model_index(size_t) const;
+    size_t                      cell_world_index(size_t) const;
     size_t                      cell_local_index_by_gmsh(size_t) const;
-    size_t                      cell_global_index_by_gmsh(size_t) const;
+    size_t                      cell_model_index_by_gmsh(size_t) const;
+    size_t                      cell_world_index_by_gmsh(size_t) const;
     size_t                      cell_local_dof_offset(size_t) const;
-    size_t                      cell_global_dof_offset(size_t) const;
+    size_t                      cell_model_dof_offset(size_t) const;
+    size_t                      cell_world_dof_offset(size_t) const;
 
     size_t                      num_faces(void) const;
     std::vector<size_t>         face_tags() const;
@@ -96,8 +107,10 @@ public:
 
     size_t                      num_dofs(void) const;
     size_t                      num_fluxes(void) const;
-    void                        base(size_t, size_t, size_t);
+    void                        base_model(size_t, size_t, size_t);
+    void                        base_world(size_t, size_t, size_t);
     size_t                      dof_base(void) const;
+    size_t                      dof_base_world(void) const;
     size_t                      flux_base(void) const;
     size_t                      index_base(void) const;
     size_t                      number(void) const;
diff --git a/include/gmsh_io.h b/include/gmsh_io.h
index 51be9f9f9245737f9fb5e2bd4fb2865cfcd3be7b..8d63a7576f8850d0b14a041b8214f5b59bcd626f 100644
--- a/include/gmsh_io.h
+++ b/include/gmsh_io.h
@@ -82,42 +82,23 @@ struct boundary_descriptor
 #define IMPLIES(a, b) ((not(a)) or (b))
 #define IFF(a,b) (IMPLIES(a,b) and IMPLIES(b,a))
 
-#ifdef USE_MPI
-struct dof_range
-{
-    int                     tag;
-    int                     dim;
-    size_t                  base;
-    std::vector<size_t>     dof_offsets;
-
-    //dof_range(dof_range&&) = default;
-
-    dof_range(const entity& e, size_t b)
-    {
-        tag = e.gmsh_tag();
-        dim = e.gmsh_dim();
-        base = b;
-        dof_offsets.resize( e.num_dofs() );
-        auto dofs_per_cell = e.num_dofs()/e.num_cells();
-        assert( dofs_per_cell*e.num_cells() == e.num_dofs() );
-        for (size_t iT = 0; iT < e.num_cells(); iT++)
-            for (size_t iD = 0; iD < dofs_per_cell; iD++)
-                dof_offsets[iT*dofs_per_cell+iD] = e.cell_global_dof_offset(iT)+iD;
-    }
-};
-
-inline std::ostream& operator<<(std::ostream& os, const dof_range& dr)
-{
-    os << "Tag " << dr.tag << ", dim " << dr.dim << ": (";
-    os << dr.base << ", " << dr.dof_offsets.size() << ")";
-    return os;
-}
-#endif
-
 class model
 {
+    /* The objects of the class 'model' are always local to the processor
+     * and they contain only the entities pertaining to the specific
+     * processor. In the 'model' class however there are additional
+     * data structures, which are needed by Rank 0 to keep track of
+     * the whole GMSH unpartitioned model. In other words, the 'model'
+     * class reflects more a GMSH partition. This mismatch in terminology
+     * is due to the fact that MPI was added in a second time.
+     * In the code, the term 'model' will thus refer to the GMSH partition,
+     * and will be a single computational unit.
+     * The collection of the models will be referred as 'world'.
+     */
+
     int                                         geometric_order;
     int                                         approximation_order;
+
 #ifdef USE_MPI
     int                                         num_partitions;
 #endif /* USE_MPI */
@@ -141,10 +122,10 @@ class model
 
     /* Map from 3D entity tag to partition */
     std::map<int, int>                          partition_inverse_map;
-public:
-    /* RANK 0 only: keeps global (process-wide) DOF numbering and positioning */
-    std::map<int, dof_range>                    global_dof_ranges;
-private:
+
+    /* RANK 0 only: keeps track of all the existing entities */
+    std::vector<int>                            all_entities_tags;
+
     /* RANK 0 only: copy of remote entities for postpro */
     std::map<int, std::vector<entity>>          remote_entities;
 #endif /* USE_MPI */
@@ -184,20 +165,25 @@ public:
     void populate_from_gmsh(void);
 
 #ifdef USE_MPI
-    const std::map<int, dof_range>& dof_ranges(void) const
+    const std::vector<int>& world_entities_tags(void) const
     {
-        return global_dof_ranges;
+        return all_entities_tags;
     }
 
-    size_t global_num_dofs() const
+    size_t num_dofs_world() const
     {
         size_t ret = 0;
-        for (auto& dr : global_dof_ranges)
-            ret += dr.second.dof_offsets.size();
+        for (auto& e : entities)
+            ret += e.num_dofs();
+        
+        for (auto& re : remote_entities)
+            for (auto& e : re.second)
+                ret += e.num_dofs();
+
         return ret;
     }
 
-    size_t global_num_cells() const
+    size_t num_cells_world() const
     {
         size_t ret = 0;
         for (auto& e : entities)
@@ -215,8 +201,8 @@ public:
         return partition_inverse_map.at(tag)-1;
     }
 
-    entity&         entity_global_lookup(int);
-    const entity&   entity_global_lookup(int) const;
+    entity&         entity_world_lookup(int);
+    const entity&   entity_world_lookup(int) const;
 #endif /* USE_MPI */
 
     const element_connectivity<element_key>& connectivity() const
diff --git a/src/entity.cpp b/src/entity.cpp
index 8e1d58f48f4e3851664572c89e8efeb3c6d95dd4..f9055efc1de666d66ed53e3ca24b216cab2bcaf0 100644
--- a/src/entity.cpp
+++ b/src/entity.cpp
@@ -9,7 +9,9 @@
 
 entity::entity(const entity_params& ep)
     : dim(ep.dim), tag(ep.tag), elemType(ep.etype), g_order(ep.gorder),
-      a_order(ep.aorder), m_dof_base(0), m_flux_base(0), m_index_base(0),
+      a_order(ep.aorder), m_dof_base_model(0), m_flux_base_model(0),
+      m_index_base_model(0), m_entity_number_model(0), m_dof_base_world(0),
+      m_flux_base_world(0), m_index_base_world(0), m_entity_number_world(0),
       cur_elem_ordering(entity_ordering::GMSH)
 {
     /* Prepare reference elements */
@@ -98,7 +100,7 @@ entity::project(const scalar_function& function) const
 }
 
 /* Project a function on this entity and store the partial result in
- * the global vector pf at the offset corresponding to this entity. */
+ * the model-local vector pf at the offset corresponding to this entity. */
 void
 entity::project(const scalar_function& function, vecxd& pf) const
 {
@@ -124,12 +126,12 @@ entity::project(const scalar_function& function, vecxd& pf) const
             rhs += pqp.weight() * function(pqp.point()) * phi; 
         }
 
-        pf.segment(m_dof_base + iT*num_bf, num_bf) = mass.ldlt().solve(rhs);
+        pf.segment(m_dof_base_model + iT*num_bf, num_bf) = mass.ldlt().solve(rhs);
     }
 }
 
 /* Project a function on this entity and store the partial result in
- * the global vector pf at the offset corresponding to this entity. */
+ * the model-local vector pf at the offset corresponding to this entity. */
 void
 entity::project(const vector_function& function, vecxd& pfx, vecxd& pfy, vecxd& pfz) const
 {
@@ -162,9 +164,9 @@ entity::project(const vector_function& function, vecxd& pfx, vecxd& pfy, vecxd&
 
         Eigen::LDLT<matxd> mass_ldlt;
         mass_ldlt.compute(mass);
-        pfx.segment(m_dof_base + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_x);
-        pfy.segment(m_dof_base + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_y);
-        pfz.segment(m_dof_base + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_z);
+        pfx.segment(m_dof_base_model + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_x);
+        pfy.segment(m_dof_base_model + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_y);
+        pfz.segment(m_dof_base_model + iT*num_bf, num_bf) = mass_ldlt.solve(rhs_z);
     }
 }
 
@@ -206,9 +208,9 @@ entity::project(const vector_function& function, double *pfx, double *pfy, doubl
 
         for (size_t i = 0; i < num_bf; i++)
         {
-            pfx[m_dof_base + iT*num_bf + i] = p_x[i];
-            pfy[m_dof_base + iT*num_bf + i] = p_y[i];
-            pfz[m_dof_base + iT*num_bf + i] = p_z[i];
+            pfx[m_dof_base_model + iT*num_bf + i] = p_x[i];
+            pfy[m_dof_base_model + iT*num_bf + i] = p_y[i];
+            pfz[m_dof_base_model + iT*num_bf + i] = p_z[i];
         }
     }
 }
@@ -267,10 +269,17 @@ entity::cell_refelem(const physical_element& pe) const
 }
 
 size_t
-entity::cell_global_index(size_t iT) const
+entity::cell_model_index(size_t iT) const
 {
     assert(iT < physical_cells.size());
-    return m_index_base + iT;
+    return m_index_base_model + iT;
+}
+
+size_t
+entity::cell_world_index(size_t iT) const
+{
+    assert(iT < physical_cells.size());
+    return m_index_base_world + iT;
 }
 
 size_t
@@ -281,12 +290,18 @@ entity::cell_local_index_by_gmsh(size_t iT) const
 }
 
 size_t
-entity::cell_global_index_by_gmsh(size_t iT) const
+entity::cell_model_index_by_gmsh(size_t iT) const
 {
     assert(iT < physical_cells.size());
-    return m_index_base + physical_cells[iT].original_position();
+    return m_index_base_model + physical_cells[iT].original_position();
 }
 
+size_t
+entity::cell_world_index_by_gmsh(size_t iT) const
+{
+    assert(iT < physical_cells.size());
+    return m_index_base_world + physical_cells[iT].original_position();
+}
 
 size_t
 entity::cell_local_dof_offset(size_t iT) const
@@ -296,10 +311,17 @@ entity::cell_local_dof_offset(size_t iT) const
 }
 
 size_t
-entity::cell_global_dof_offset(size_t iT) const
+entity::cell_model_dof_offset(size_t iT) const
+{
+    assert(reference_cells.size() > 0);
+    return m_dof_base_model + iT * reference_cells[0].num_basis_functions();
+}
+
+size_t
+entity::cell_world_dof_offset(size_t iT) const
 {
     assert(reference_cells.size() > 0);
-    return m_dof_base + iT * reference_cells[0].num_basis_functions();
+    return m_dof_base_world + iT * reference_cells[0].num_basis_functions();
 }
 
 const reference_element&
@@ -369,41 +391,55 @@ entity::num_fluxes(void) const
 }
 
 void
-entity::base(size_t d_base, size_t f_base, size_t i_base)
+entity::base_model(size_t d_base, size_t f_base, size_t i_base)
 {
-    m_dof_base = d_base;
-    m_flux_base = f_base;
-    m_index_base = i_base;
+    m_dof_base_model = d_base;
+    m_flux_base_model = f_base;
+    m_index_base_model = i_base;
+}
+
+void
+entity::base_world(size_t d_base, size_t f_base, size_t i_base)
+{
+    m_dof_base_world = d_base;
+    m_flux_base_world = f_base;
+    m_index_base_world = i_base;
 }
 
 size_t
 entity::dof_base(void) const
 {
-    return m_dof_base;
+    return m_dof_base_model;
+}
+
+size_t
+entity::dof_base_world(void) const
+{
+    return m_dof_base_world;
 }
 
 size_t
 entity::flux_base(void) const
 {
-    return m_flux_base;
+    return m_flux_base_model;
 }
 
 size_t
 entity::index_base(void) const
 {
-    return m_index_base;
+    return m_index_base_model;
 }
 
 size_t
 entity::number(void) const
 {
-    return m_entity_number;
+    return m_entity_number_model;
 }
 
 void
 entity::number(size_t n)
 {
-    m_entity_number = n;
+    m_entity_number_model = n;
 }
 
 matxd
@@ -848,8 +884,8 @@ entity::populate_entity_data(entity_data_cpu& ed, const model& mod) const
     ed.num_fluxes = reference_faces[0].num_basis_functions();
     ed.num_qp = (g_order == 1) ? 1 : reference_cells[0].num_integration_points();
     ed.num_face_qp = (g_order == 1) ? 1 : reference_faces[0].num_integration_points();
-    ed.dof_base = m_dof_base;
-    ed.flux_base = m_flux_base;
+    ed.dof_base = m_dof_base_model;
+    ed.flux_base = m_flux_base_model;
     ed.num_faces_per_elem = num_faces_per_elem();
 
     populate_differentiation_matrices(ed);
@@ -885,7 +921,7 @@ entity::populate_entity_data(entity_data_cpu& ed, const model& mod) const
             for (size_t i = 0; i < face_keys.size(); i++)
                 for (size_t j = 0; j < cell_keys.size(); j++)
                     if ( face_keys.at(i) == cell_keys.at(j) )
-                        ed.fluxdofs_mine.at(face_base + i) = cell_global_dof_offset(iT) + j;
+                        ed.fluxdofs_mine.at(face_base + i) = cell_model_dof_offset(iT) + j;
 
             element_key fk(pf);
             
@@ -897,7 +933,7 @@ entity::populate_entity_data(entity_data_cpu& ed, const model& mod) const
                 for (size_t i = 0; i < face_keys.size(); i++)
                     for (size_t j = 0; j < ne_cell_keys.size(); j++)
                         if ( face_keys.at(i) == ne_cell_keys.at(j) )
-                            ed.fluxdofs_neigh.at(face_base + i) = ne_e.cell_global_dof_offset(neigh.iT) + j;
+                            ed.fluxdofs_neigh.at(face_base + i) = ne_e.cell_model_dof_offset(neigh.iT) + j;
             }
         }
     }
diff --git a/src/gmsh_io.cpp b/src/gmsh_io.cpp
index 358ea89d1770589e8dbd8e602551ba9f412f9125..e4038779f0b5896a0741c7724e884d921053145b 100644
--- a/src/gmsh_io.cpp
+++ b/src/gmsh_io.cpp
@@ -195,13 +195,14 @@ model::import_gmsh_entities_rank0(void)
     gvp_t gmsh_entities;
     gm::getEntities(gmsh_entities, DIMENSION(3));
 
-    size_t entity_index = 0;
-    size_t dof_base = 0;
-    size_t flux_base = 0;
-    size_t index_base = 0;
-#ifdef USE_MPI
-    size_t global_dof_base = 0;
-#endif /* USE_MPI */
+    size_t entity_index_model = 0;
+    size_t dof_base_model = 0;
+    size_t flux_base_model = 0;
+    size_t index_base_model = 0;
+    size_t entity_index_world = 0;
+    size_t dof_base_world = 0;
+    size_t flux_base_world = 0;
+    size_t index_base_world = 0;
 
     for (auto [dim, tag] : gmsh_entities)
     {
@@ -245,16 +246,22 @@ model::import_gmsh_entities_rank0(void)
             
                 if (entity_is_remote)
                 {
-                    dof_range dr(e, global_dof_base);
-                    global_dof_ranges.emplace(tag, std::move(dr));
-                    global_dof_base += e.num_dofs();
-
+                    /* Update world-offsets */
+                    e.base_world(dof_base_world, flux_base_world, index_base_world);
+                    //e.number(entity_index_world);
+                    dof_base_world += e.num_dofs();
+                    flux_base_world += e.num_fluxes();
+                    index_base_world += e.num_cells();
+
+                    /* Send entity to the appropriate process */
                     assert(tag_partition > 0);
                     int rank = tag_partition - 1;
                     e.mpi_send(rank, MPI_COMM_WORLD);
                     std::cout << "Sending entity " << tag << " to " << rank << " " << e.num_dofs() << std::endl;
                     /* Save locally the entity */
                     remote_entities[tag_partition].push_back( std::move(e) );
+                    all_entities_tags.push_back(tag);
+                    entity_index_world++;
                     continue;
                 }
             }
@@ -267,24 +274,30 @@ model::import_gmsh_entities_rank0(void)
                 if ( eitor != etag_to_entity_offset.end() )
                     throw std::logic_error("Duplicate tag");
 
-                etag_to_entity_offset[ pe.element_tag() ] = std::make_pair(entity_index, i);
+                etag_to_entity_offset[ pe.element_tag() ] = std::make_pair(entity_index_model, i);
             }
 
-            e.base(dof_base, flux_base, index_base);
-            e.number(entity_index);
-            dof_base += e.num_dofs();
-            flux_base += e.num_fluxes();
-            index_base += e.num_cells();
-
+            /* Update model-offsets */
+            e.base_model(dof_base_model, flux_base_model, index_base_model);
+            e.number(entity_index_model);
+            dof_base_model += e.num_dofs();
+            flux_base_model += e.num_fluxes();
+            index_base_model += e.num_cells();
+
+            /* Update world-offsets */
+            e.base_world(dof_base_world, flux_base_world, index_base_world);
+            //e.number(entity_index_world);
+            dof_base_world += e.num_dofs();
+            flux_base_world += e.num_fluxes();
+            index_base_world += e.num_cells();
+
+            update_connectivity(e, entity_index_model);
+            entities.push_back( std::move(e) );
 #ifdef USE_MPI
-            dof_range dr(e, global_dof_base);
-            global_dof_ranges.emplace(tag, std::move(dr));
-            global_dof_base += e.num_dofs();
+            all_entities_tags.push_back(tag);
 #endif /* USE_MPI */
-
-            update_connectivity(e, entity_index);
-            entities.push_back( std::move(e) );
-            entity_index++;
+            entity_index_model++;
+            entity_index_world++;
         }
     }
 
@@ -293,20 +306,21 @@ model::import_gmsh_entities_rank0(void)
     assert( IMPLIES(num_partitions == 1, remote_entities.size() == 0) );
     for (auto& re : remote_entities)
     {
-        entity_index = 0;
-        dof_base = 0;
-        flux_base = 0;
-        index_base = 0;
+        entity_index_model = 0;
+        dof_base_model = 0;
+        flux_base_model = 0;
+        index_base_model = 0;
         for (auto& e : re.second)
         {
-            e.base(dof_base, flux_base, index_base);
-            e.number(entity_index);
-            dof_base += e.num_dofs();
-            flux_base += e.num_fluxes();
-            index_base += e.num_cells();
-            entity_index++;
+            e.base_model(dof_base_model, flux_base_model, index_base_model);
+            e.number(entity_index_model);
+            dof_base_model += e.num_dofs();
+            flux_base_model += e.num_fluxes();
+            index_base_model += e.num_cells();
+            entity_index_model++;
         }
     }
+    std::sort(all_entities_tags.begin(), all_entities_tags.end());
 #endif
 }
 
@@ -337,7 +351,7 @@ model::import_gmsh_entities_rankN(void)
             etag_to_entity_offset[ pe.element_tag() ] = std::make_pair(entity_index, i);
         }
 
-        e.base(dof_base, flux_base, index_base);
+        e.base_model(dof_base, flux_base, index_base);
         e.number(entity_index);
         dof_base += e.num_dofs();
         flux_base += e.num_fluxes();
@@ -585,7 +599,7 @@ model::partition_mesh(int parts)
 }
 
 entity&
-model::entity_global_lookup(int tag)
+model::entity_world_lookup(int tag)
 {
     /* Look for tag between local entities */
     for (auto& e : entities)
@@ -602,7 +616,7 @@ model::entity_global_lookup(int tag)
 }
 
 const entity&
-model::entity_global_lookup(int tag) const
+model::entity_world_lookup(int tag) const
 {
     /* Look for tag between local entities */
     for (const auto& e : entities)
diff --git a/src/maxwell_common.cpp b/src/maxwell_common.cpp
index e731922adb158bb0a7e9cd44ad25469575a22de6..3a8bf17b3de6fe37ec3ecf8342768d69b9ac040c 100644
--- a/src/maxwell_common.cpp
+++ b/src/maxwell_common.cpp
@@ -156,7 +156,7 @@ void export_E_field_nodal(const model& mod, silo& sdb, const vecxd& Ex,
         for (size_t iT = 0; iT < e.num_cells(); iT++)
         {
             auto& pe = e.cell(iT);
-            auto cgofs = e.cell_global_dof_offset(iT);
+            auto cgofs = e.cell_model_dof_offset(iT);
             auto nodetags = pe.node_tags();
             for (size_t iD = 0; iD < 4; iD++)
             {
@@ -196,7 +196,7 @@ void export_H_field_nodal(const model& mod, silo& sdb, const vecxd& Hx,
         for (size_t iT = 0; iT < e.num_cells(); iT++)
         {
             auto& pe = e.cell(iT);
-            auto cgofs = e.cell_global_dof_offset(iT);
+            auto cgofs = e.cell_model_dof_offset(iT);
             auto nodetags = pe.node_tags();
             for (size_t iD = 0; iD < 4; iD++)
             {
@@ -239,7 +239,7 @@ void export_J_field_nodal(const model& mod, silo& sdb, const vecxd& Ex,
             auto& pe = e.cell(iT);
             auto bar = pe.barycenter();
             auto sigma = mpl.sigma(etag, bar);
-            auto cgofs = e.cell_global_dof_offset(iT);
+            auto cgofs = e.cell_model_dof_offset(iT);
             auto nodetags = pe.node_tags();
             for (size_t iD = 0; iD < 4; iD++)
             {
@@ -284,8 +284,8 @@ void export_E_field_zonal(const model& mod, silo& sdb, const vecxd& Ex,
             auto& re = e.cell_refelem(pe);
             vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
             auto num_bf = re.num_basis_functions();
-            auto ofs = e.cell_global_dof_offset(iT);
-            auto gi = e.cell_global_index_by_gmsh(iT);
+            auto ofs = e.cell_model_dof_offset(iT);
+            auto gi = e.cell_model_index_by_gmsh(iT);
             out_Ex[gi] = Ex.segment(ofs, num_bf).dot(phi_bar);
             out_Ey[gi] = Ey.segment(ofs, num_bf).dot(phi_bar);
             out_Ez[gi] = Ez.segment(ofs, num_bf).dot(phi_bar);
@@ -314,8 +314,8 @@ void export_H_field_zonal(const model& mod, silo& sdb, const vecxd& Hx,
             auto& re = e.cell_refelem(pe);
             vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
             auto num_bf = re.num_basis_functions();
-            auto ofs = e.cell_global_dof_offset(iT);
-            auto gi = e.cell_global_index_by_gmsh(iT);
+            auto ofs = e.cell_model_dof_offset(iT);
+            auto gi = e.cell_model_index_by_gmsh(iT);
             out_Hx[gi] = Hx.segment(ofs, num_bf).dot(phi_bar);
             out_Hy[gi] = Hy.segment(ofs, num_bf).dot(phi_bar);
             out_Hz[gi] = Hz.segment(ofs, num_bf).dot(phi_bar);
@@ -347,8 +347,8 @@ void export_J_field_zonal(const model& mod, silo& sdb, const vecxd& Ex,
             auto sigma = mpl.sigma(etag, bar);
             vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
             auto num_bf = re.num_basis_functions();
-            auto ofs = e.cell_global_dof_offset(iT);
-            auto gi = e.cell_global_index_by_gmsh(iT);
+            auto ofs = e.cell_model_dof_offset(iT);
+            auto gi = e.cell_model_index_by_gmsh(iT);
             out_Jx[gi] = sigma*Ex.segment(ofs, num_bf).dot(phi_bar);
             out_Jy[gi] = sigma*Ey.segment(ofs, num_bf).dot(phi_bar);
             out_Jz[gi] = sigma*Ez.segment(ofs, num_bf).dot(phi_bar);
diff --git a/src/maxwell_cpu.cpp b/src/maxwell_cpu.cpp
index d0b49e4db42205ef33cdf8896a22ba21bda2fff4..942b9cd8d0b3e597e11a5ef2769b3ac587f74cb8 100644
--- a/src/maxwell_cpu.cpp
+++ b/src/maxwell_cpu.cpp
@@ -31,7 +31,7 @@ init_matparams(const model& mod, solver_state& state,
             auto epsilon = mpl.epsilon(etag, bar);
             auto mu = mpl.mu(etag, bar);
             auto sigma = mpl.sigma(etag, bar);
-            auto ofs = e.cell_global_dof_offset(iT);
+            auto ofs = e.cell_model_dof_offset(iT);
 
             for (size_t iD = 0; iD < re.num_basis_functions(); iD++)
             {
@@ -600,39 +600,37 @@ gather_field(const model& mod, maxwell::field& in, maxwell::field& out,
 
     if (rank == root)
     {
-        auto& gdr = mod.dof_ranges();
-        for (auto& dr : gdr)
+        auto& all_tags = mod.world_entities_tags();
+        for (auto& tag : all_tags)
         {
-            auto tag = dr.second.tag;
             if (mod.entity_rank(tag) == root)
                 continue;
 
-            auto base = dr.second.base;
-            auto len = dr.second.dof_offsets.size();
-            //std::cout << "GR: " << tag << " " << len << std::endl;
+            auto& e = mod.entity_world_lookup(tag);
+            auto base = e.dof_base_world();
+            auto len = e.num_dofs();
+            std::cout << "GR: " << tag << " " << len << std::endl;
             receive_field(out, mod.entity_rank(tag), tag, base, len, comm);
         }
 
         for (auto& e : mod)
         {
             auto in_base = e.dof_base();
-            auto in_len = e.num_dofs();
-            auto out_base = gdr.at(e.gmsh_tag()).base;
-            auto out_len = gdr.at(e.gmsh_tag()).dof_offsets.size();
-            assert(in_len == out_len);
-            out.Ex.segment(out_base, out_len) = in.Ex.segment(in_base, in_len);
-            out.Ey.segment(out_base, out_len) = in.Ey.segment(in_base, in_len);
-            out.Ez.segment(out_base, out_len) = in.Ez.segment(in_base, in_len);
-            out.Hx.segment(out_base, out_len) = in.Hx.segment(in_base, in_len);
-            out.Hy.segment(out_base, out_len) = in.Hy.segment(in_base, in_len);
-            out.Hz.segment(out_base, out_len) = in.Hz.segment(in_base, in_len);
+            auto len = e.num_dofs();
+            auto out_base = e.dof_base_world();
+            out.Ex.segment(out_base, len) = in.Ex.segment(in_base, len);
+            out.Ey.segment(out_base, len) = in.Ey.segment(in_base, len);
+            out.Ez.segment(out_base, len) = in.Ez.segment(in_base, len);
+            out.Hx.segment(out_base, len) = in.Hx.segment(in_base, len);
+            out.Hy.segment(out_base, len) = in.Hy.segment(in_base, len);
+            out.Hz.segment(out_base, len) = in.Hz.segment(in_base, len);
         }
     }
     else
     {
         for (auto& e : mod)
         {
-            //std::cout << "GS: " << e.gmsh_tag() << " " << e.num_dofs() << std::endl;
+            std::cout << "GS: " << e.gmsh_tag() << " " << e.num_dofs() << std::endl;
             send_field(in, root, e.gmsh_tag(), e.dof_base(), e.num_dofs(), comm);
         }
     }
@@ -648,7 +646,7 @@ export_fields_to_silo_2(const model& mod, maxwell::solver_state& state,
 
     if (rank == 0)
     {
-        size_t gdnum = mod.global_num_dofs();
+        size_t gdnum = mod.num_dofs_world();
         maxwell::field f;
         f.resize(gdnum);
         gather_field(mod, state.emf_next, f, 0, MPI_COMM_WORLD);
@@ -661,29 +659,25 @@ export_fields_to_silo_2(const model& mod, maxwell::solver_state& state,
         sdb.import_mesh_from_gmsh();
         sdb.write_mesh(state.curr_time, state.curr_timestep);
 
-        std::vector<double> out_Ex( mod.global_num_cells() );
-        std::vector<double> out_Ey( mod.global_num_cells() );
-        std::vector<double> out_Ez( mod.global_num_cells() );
+        std::vector<double> out_Ex( mod.num_cells_world() );
+        std::vector<double> out_Ey( mod.num_cells_world() );
+        std::vector<double> out_Ez( mod.num_cells_world() );
 
-        size_t nc = 0;
-        for (auto& dr : mod.dof_ranges())
+        for (auto& tag : mod.world_entities_tags())
         {
-            auto etag = dr.first;
-            std::cout << etag << std::endl;
-            auto& e = mod.entity_global_lookup(etag);
+            auto& e = mod.entity_world_lookup(tag);
             for (size_t iT = 0; iT < e.num_cells(); iT++)
             {
                 auto& pe = e.cell(iT);
                 auto& re = e.cell_refelem(pe);
                 vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
                 auto num_bf = re.num_basis_functions();
-                auto ofs = dr.second.base + num_bf*iT;
-                auto gi = nc+e.cell_local_index_by_gmsh(iT);
+                auto ofs = e.cell_world_dof_offset(iT);
+                auto gi = e.cell_world_index_by_gmsh(iT);
                 out_Ex[gi] = f.Ex.segment(ofs, num_bf).dot(phi_bar);
                 out_Ey[gi] = f.Ey.segment(ofs, num_bf).dot(phi_bar);
                 out_Ez[gi] = f.Ez.segment(ofs, num_bf).dot(phi_bar);
             }
-            nc += e.num_cells();
         }
 
         sdb.write_zonal_variable("Ex", out_Ex);
diff --git a/src/maxwell_solver.cpp b/src/maxwell_solver.cpp
index 7c705643b533b198631008d69e115958b0af8632..9653c9ba196871876dc42f666a703199771b3998 100644
--- a/src/maxwell_solver.cpp
+++ b/src/maxwell_solver.cpp
@@ -139,7 +139,7 @@ integrate_electric_field(const model& mod, const maxwell::solver_state& state,
         auto& e = mod.entity_at(entnum);
         auto& pe = e.cell(ofs);
         auto& re = e.cell_refelem(pe);
-        auto dof_ofs = e.cell_global_dof_offset(ofs);
+        auto dof_ofs = e.cell_model_dof_offset(ofs);
         auto dof_num = re.num_basis_functions();
         vecxd phi = re.basis_functions( refp );
         vecxd locEx = state.emf_curr.Ex.segment(dof_ofs, dof_num);
@@ -175,7 +175,7 @@ reinit_materials(const model& mod, maxwell::solver_state& state, maxwell::parame
             auto epsilon = mpl.epsilon(tag, bar);
             //auto mu = mpl.mu(tag, bar);
             auto sigma = mpl.sigma(tag, bar);
-            auto ofs = e.cell_global_dof_offset(iT);
+            auto ofs = e.cell_model_dof_offset(iT);
 
             for (size_t iD = 0; iD < re.num_basis_functions(); iD++)
             {
diff --git a/tests/test_differentiation.cpp b/tests/test_differentiation.cpp
index 884d8f56f95c364f34e08b780a0abf3302a85c86..6e792c2d979fa1a21c471451ab47fb001b517c5c 100644
--- a/tests/test_differentiation.cpp
+++ b/tests/test_differentiation.cpp
@@ -147,7 +147,7 @@ int test_differentiation_convergence(int geometric_order, int approximation_orde
                 auto& re = e.cell_refelem(pe);
                 matxd mass = e.mass_matrix(iT);
                 auto num_bf = re.num_basis_functions();
-                auto ofs = e.cell_global_dof_offset(iT);
+                auto ofs = e.cell_model_dof_offset(iT);
                 vecxd diff_x = Pdf_dx.segment(ofs, num_bf) - Cdf_dx.segment(ofs, num_bf); 
                 vecxd diff_y = Pdf_dy.segment(ofs, num_bf) - Cdf_dy.segment(ofs, num_bf); 
                 vecxd diff_z = Pdf_dz.segment(ofs, num_bf) - Cdf_dz.segment(ofs, num_bf); 
diff --git a/tests/test_differentiation_gpu.cpp b/tests/test_differentiation_gpu.cpp
index c9afbd298d2a3cf870712fb8bfb71544efaf74a9..ebf9c4f585db9c4e45ea170cb77bbdd7bbe9a627 100644
--- a/tests/test_differentiation_gpu.cpp
+++ b/tests/test_differentiation_gpu.cpp
@@ -205,7 +205,7 @@ int test_differentiation_convergence(int geometric_order, int approximation_orde
                 auto& re = e.cell_refelem(pe);
                 matxd mass = e.mass_matrix(iT);
                 auto num_bf = re.num_basis_functions();
-                auto ofs = e.cell_global_dof_offset(iT);
+                auto ofs = e.cell_model_dof_offset(iT);
                 vecxd diff_x = Pdf_dx.segment(ofs, num_bf) - Cdf_dx.segment(ofs, num_bf); 
                 vecxd diff_y = Pdf_dy.segment(ofs, num_bf) - Cdf_dy.segment(ofs, num_bf); 
                 vecxd diff_z = Pdf_dz.segment(ofs, num_bf) - Cdf_dz.segment(ofs, num_bf); 
diff --git a/tests/test_lifting.cpp b/tests/test_lifting.cpp
index a543f9080da4d427ec9c61a3f47678f526bb281a..022b562d2aa4e071352f3623fd49d56914db39b1 100644
--- a/tests/test_lifting.cpp
+++ b/tests/test_lifting.cpp
@@ -198,7 +198,7 @@ int test_lifting(int geometric_order, int approximation_order)
                     vol += pqp.weight() * (dphi*J.inverse()) * F(pqp.point());
                 }
 
-                auto ofs = e.cell_global_dof_offset(iT);
+                auto ofs = e.cell_model_dof_offset(iT);
                 vecxd excpected_div_F = Pdiv_F.segment(ofs, num_bf);
                 vecxd Plf = LiftF.segment(ofs, num_bf);
                 vecxd Pvol = mass.ldlt().solve(vol);
diff --git a/tests/test_lifting_gpu.cpp b/tests/test_lifting_gpu.cpp
index cd289c8390a77cf7cc933f8672bde94d7861e48e..82bce0d704d591b9aa7f4df317b0dba35278fbd5 100644
--- a/tests/test_lifting_gpu.cpp
+++ b/tests/test_lifting_gpu.cpp
@@ -207,7 +207,7 @@ int test_lifting(int geometric_order, int approximation_order)
                     vol += pqp.weight() * (dphi*J.inverse()) * F(pqp.point());
                 }
 
-                auto ofs = e.cell_global_dof_offset(iT);
+                auto ofs = e.cell_model_dof_offset(iT);
                 vecxd excpected_div_F = Pdiv_F.segment(ofs, num_bf);
                 vecxd Plf = LiftF.segment(ofs, num_bf);
                 vecxd Pvol = mass.ldlt().solve(vol);