Entities are distributed.

include/entity.h

@@ -56,7 +56,7 @@ class entity
     void populate_normals(entity_data_cpu&) const;
 
 public:
-    entity() = delete;
+    entity() {}// = delete;
     entity(const entity&) = delete;
     entity(entity&&) = default;
     entity& operator=(const entity&) = delete;

include/gmsh_io.h

@@ -36,7 +36,9 @@ enum class face_type : int
     NONE = 0,
     INTERFACE,
     BOUNDARY,
+#ifdef USE_MPI
     INTERPROCESS_BOUNDARY
+#endif
 };
 
 struct boundary_descriptor
@@ -84,6 +86,9 @@ class model
     /* Map from the partition number to the entities belonging to that partition */
     std::map<int, std::vector<int>>             partition_map;
 
+    /* Map from 3D entity tag to partition */
+    std::map<int, int>                          partition_inverse_map;
+
     std::vector<entity>                         entities;
     std::vector<boundary_descriptor>            bnd_descriptors;
     element_connectivity<element_key>           conn;
@@ -92,13 +97,16 @@ class model
     std::map<size_t, entofs_pair>               etag_to_entity_offset;
 
     void    map_boundaries(void);
-    void    import_gmsh_entities(void);
+    void    import_gmsh_entities_rank0(void);
     void    update_connectivity(const entity&, size_t);
 
     void    populate_from_gmsh_rank0(void);
 
 #ifdef USE_MPI
+    bool    is_interprocess_boundary(int);
     void    populate_from_gmsh_rankN(void);
+    void    import_gmsh_entities_rankN(void);
+    int     my_partition(void) const;
 #endif
 
     void    make_boundary_to_faces_map(void);

include/gmsh_mpi.h

@@ -14,7 +14,7 @@
         assert(__rank == 0);                    \
     }
 #else
-    #define ASSERT_MPI_RANK_0 (void) 0
+    #define ASSERT_MPI_RANK_0 ((void) 0)
 #endif
 
 template<typename T>

src/gmsh_io.cpp

@@ -57,7 +57,9 @@ model::~model()
 void
 model::make_boundary_to_faces_map(void)
 {
+#ifdef USE_MPI
     ASSERT_MPI_RANK_0;
+#endif /* USE_MPI */
 
     gvp_t gmsh_entities;
 
@@ -105,7 +107,9 @@ model::make_boundary_to_faces_map(void)
 void
 model::make_partition_to_entities_map()
 {
+#ifdef USE_MPI
     ASSERT_MPI_RANK_0;
+#endif /* USE_MPI */
 
     gvp_t gmsh_entities;
 
@@ -127,6 +131,7 @@ model::make_partition_to_entities_map()
             assert(parts.size() == 1);
             gm::getPartitions(dim, tag, parts);
             partition_map[ parts[0] ].push_back(tag);
+            partition_inverse_map[tag] = parts[0];
         }
     }
 
@@ -137,7 +142,9 @@ model::make_partition_to_entities_map()
 void
 model::update_connectivity(const entity& e, size_t entity_index)
 {
+#ifdef USE_MPI
     ASSERT_MPI_RANK_0;
+#endif /* USE_MPI */
 
     for (size_t iT = 0; iT < e.num_cells(); iT++)
     {
@@ -151,9 +158,25 @@ model::update_connectivity(const entity& e, size_t entity_index)
     }
 }
 
+#ifdef USE_MPI
+int
+model::my_partition(void) const
+{
+    /* Here we are assuming that GMSH enumerates mesh partitions
+     * sequentially and starting from 1. */
+    int rank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+    return rank+1;
+}
+#endif /* USE_MPI */
+
 void
-model::import_gmsh_entities(void)
+model::import_gmsh_entities_rank0(void)
 {
+#ifdef USE_MPI
+    ASSERT_MPI_RANK_0;
+#endif /* USE_MPI */
+
     gvp_t gmsh_entities;
     gm::getEntities(gmsh_entities, DIMENSION(3));
 
@@ -189,6 +212,24 @@ model::import_gmsh_entities(void)
             
             e.sort_by_orientation();
 
+#ifdef USE_MPI
+            if (num_partitions > 1)
+            {
+                auto mp = my_partition();
+                auto tag_partition = partition_inverse_map.at(tag);
+                bool entity_is_remote = ( tag_partition != mp );
+            
+                if (entity_is_remote)
+                {
+                    assert(tag_partition > 0);
+                    int rank = tag_partition - 1;
+                    e.mpi_send(rank, MPI_COMM_WORLD);
+                    std::cout << "Sending entity " << tag << " to " << rank << std::endl;
+                    continue;
+                }
+            }
+            std::cout << "Keeping entity " << tag << std::endl;
+#endif /* USE_MPI */
             for (size_t i = 0; i < e.num_cells(); i++)
             {
                 const auto& pe = e.cell(i);
@@ -211,12 +252,65 @@ model::import_gmsh_entities(void)
             entity_index++;
         }
     }
+    std::cout << dof_base << " " << flux_base << std::endl;
+}
+
+#ifdef USE_MPI
+void
+model::import_gmsh_entities_rankN(void)
+{
+    int mp = my_partition();
+
+    size_t entity_index = 0;
+    size_t dof_base = 0;
+    size_t flux_base = 0;
+    size_t index_base = 0;
+
+    for (auto& rtag : partition_map.at(mp))
+    {
+        std::cout << "Receiving " << rtag << ", rank " << mp-1 << std::endl;
+        entity e;
+        e.mpi_recv(0, MPI_COMM_WORLD);
+
+        for (size_t i = 0; i < e.num_cells(); i++)
+        {
+            const auto& pe = e.cell(i);
+            auto eitor = etag_to_entity_offset.find( pe.element_tag() );
+            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);
+        }
+
+        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_connectivity(e, entity_index);
+
+        entities.push_back( std::move(e) );
+        entity_index++;
+    }
+
+    std::cout << dof_base << " " << flux_base << std::endl;
 }
 
+bool
+model::is_interprocess_boundary(int tag)
+{
+    return comm_map.at(tag).size() == 2;
+}
+
+#endif /* USE_MPI */
+
 void
 model::map_boundaries(void)
 {
+#ifdef USE_MPI
     ASSERT_MPI_RANK_0;
+#endif /* USE_MPI */
 
     /* Make a vector mapping element_key to entity tag */
     using bfk_t = std::pair<element_key, int>;
@@ -253,9 +347,18 @@ model::map_boundaries(void)
             boundary_descriptor bd;
             /* determine if it is an interface or a boundary */
             if (conn.num_neighbours(fk) == 1)
-                bd.type = face_type::BOUNDARY;
+            {
+#ifdef USE_MPI
+                if ( is_interprocess_boundary( (*itor).second ) )
+                    bd.type = face_type::INTERPROCESS_BOUNDARY;
+                else
+#endif /* USE_MPI */
+                    bd.type = face_type::BOUNDARY;
+            }
             else
+            {
                 bd.type = face_type::INTERFACE;
+            }
             /* and store the gmsh tag in the descriptor. */
             bd.gmsh_entity = (*itor).second;
 
@@ -277,7 +380,9 @@ model::map_boundaries(void)
             case face_type::NONE: normal++; break;
             case face_type::INTERFACE: interface++; break;
             case face_type::BOUNDARY: boundary++; break;
+#ifdef USE_MPI
             case face_type::INTERPROCESS_BOUNDARY: ipc_boundary++; break;
+#endif /* USE_MPI */
         }
     }
 
@@ -288,7 +393,9 @@ model::map_boundaries(void)
 void
 model::populate_from_gmsh_rank0()
 {
+#ifdef USE_MPI
     ASSERT_MPI_RANK_0;
+#endif /* USE_MPI */
 
     gmm::setOrder( geometric_order );
     make_boundary_to_faces_map();
@@ -321,7 +428,7 @@ model::populate_from_gmsh_rank0()
 #endif /* USE_MPI */
 
     entities.clear();
-    import_gmsh_entities();
+    import_gmsh_entities_rank0();
     map_boundaries(); /* This must happen after import_gmsh_entities(). */
 }
 
@@ -344,6 +451,10 @@ model::populate_from_gmsh_rankN()
             std::cout << v << " ";
         std::cout << std::endl;
     }
+
+    entities.clear();
+    import_gmsh_entities_rankN();
+    map_boundaries();
 }
 #endif /* USE_MPI */