Rank 0 solver looks OK.

include/maxwell_common.h

@@ -22,7 +22,11 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
         {
             auto gfnum = face_num_base + iF;
             auto bd = bds[gfnum];
-            if (bd.type == face_type::NONE or bd.type == face_type::INTERFACE)
+            if (bd.type == face_type::NONE or bd.type == face_type::INTERFACE
+#ifdef USE_MPI
+                    or bd.type == face_type::INTERPROCESS_BOUNDARY
+#endif /* USE_MPI */
+               )
                 continue;
 
             auto& pf = e.face(iF);
@@ -55,10 +59,10 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
 
                 case boundary_condition::PLANE_WAVE_E: {
                     auto fE = [&](const point_3d& pt) -> vec3d {
-                        return mpl.eval_boundary_source(bd.gmsh_entity, pt, state.curr_time);
+                        return mpl.eval_boundary_source(bd.material_tag(), pt, state.curr_time);
                     };
                     auto fH = [&](const point_3d& pt) -> vec3d {
-                        vec3d H = mpl.eval_boundary_source(bd.gmsh_entity, pt, state.curr_time);
+                        vec3d H = mpl.eval_boundary_source(bd.material_tag(), pt, state.curr_time);
                         return n.cross(H)/Z;
                     };
                     matxd prjE = e.project_on_face(iF, fE);
@@ -106,12 +110,16 @@ eval_interface_sources(const model& mod, const parameter_loader& mpl,
         {
             auto gfnum = face_num_base + iF;
             auto bd = bds[gfnum];
-            if (bd.type == face_type::NONE or bd.type == face_type::BOUNDARY)
+            if (bd.type == face_type::NONE or bd.type == face_type::BOUNDARY
+#ifdef USE_MPI
+                    or bd.type == face_type::INTERPROCESS_BOUNDARY
+#endif /* USE_MPI */
+               )
                 continue;
 
             vec3d n = state.eds[e.number()].normals.row(iF);
 
-            switch ( mpl.interface_type(bd.gmsh_entity) )
+            switch ( mpl.interface_type(bd.material_tag()) )
             {
                 case interface_condition::UNSPECIFIED:
                     break;
@@ -119,7 +127,7 @@ eval_interface_sources(const model& mod, const parameter_loader& mpl,
                 case interface_condition::E_FIELD: {
                     /* This is BS actually... */
                     auto fE = [&](const point_3d& pt) -> vec3d {
-                        return mpl.eval_interface_source(bd.gmsh_entity, pt, state.curr_time);
+                        return mpl.eval_interface_source(bd.material_tag(), pt, state.curr_time);
                     };
 
                     matxd prjE = e.project_on_face(iF, fE);
@@ -134,7 +142,7 @@ eval_interface_sources(const model& mod, const parameter_loader& mpl,
 
                 case interface_condition::SURFACE_CURRENT: {
                     auto fJ = [&](const point_3d& pt) -> vec3d {
-                        vec3d J = mpl.eval_interface_source(bd.gmsh_entity, pt, state.curr_time);
+                        vec3d J = mpl.eval_interface_source(bd.material_tag(), pt, state.curr_time);
                         return n.cross(J);
                     };
 

src/gmsh_io.cpp

@@ -87,6 +87,7 @@ model::make_boundary_to_faces_map(void)
             std::vector<int> parts;
             gm::getPartitions(dim, tag, parts);
             comm_map[tag] = parts;
+            if (pdim == dim)
                 surface_to_parent_map[tag] = ptag;
         }
 
@@ -102,6 +103,10 @@ model::make_boundary_to_faces_map(void)
         }
     }
 
+
+    for (auto& sp : surface_to_parent_map)
+        std::cout << sp.first << " -> " << sp.second << std::endl;
+
     assert( IFF(comm_map.size() == 0, num_partitions == 1) );
 }
 
@@ -129,8 +134,8 @@ model::make_partition_to_entities_map()
                 continue;
 
             std::vector<int> parts;
-            assert(parts.size() == 1);
             gm::getPartitions(dim, tag, parts);
+            assert(parts.size() == 1);
             partition_map[ parts[0] ].push_back(tag);
             partition_inverse_map[tag] = parts[0];
         }
@@ -144,7 +149,7 @@ void
 model::update_connectivity(const entity& e, size_t entity_index)
 {
 #ifdef USE_MPI
-    ASSERT_MPI_RANK_0;
+    //ASSERT_MPI_RANK_0;
 #endif /* USE_MPI */
 
     for (size_t iT = 0; iT < e.num_cells(); iT++)
@@ -313,7 +318,7 @@ void
 model::map_boundaries(void)
 {
 #ifdef USE_MPI
-    ASSERT_MPI_RANK_0;
+    //ASSERT_MPI_RANK_0;
 #endif /* USE_MPI */
 
     /* Make a vector mapping element_key to entity tag */
@@ -366,8 +371,11 @@ model::map_boundaries(void)
             /* and store the gmsh tag in the descriptor. */
             bd.gmsh_entity = (*itor).second;
             if (num_partitions > 1)
+            {
+                auto sitor = surface_to_parent_map.find( (*itor).second );
+                if ( sitor != surface_to_parent_map.end() )
                     bd.parent_entity = surface_to_parent_map.at( (*itor).second );
-
+            }
             /* Finally, put the descriptor in the global array of faces. */
             bnd_descriptors.at(fbase + iF) = bd;
         }

src/maxwell_common.cpp

@@ -29,7 +29,7 @@ init_boundary_conditions(const model& mod, const parameter_loader& mpl,
                 continue;
 
             double bc_coeff = DIRICHLET;
-            switch ( mpl.boundary_type(bd.gmsh_entity) )
+            switch ( mpl.boundary_type(bd.material_tag()) )
             {
                 case boundary_condition::UNSPECIFIED:
                 case boundary_condition::PEC:

src/maxwell_interface.cpp

@@ -395,12 +395,12 @@ parameter_loader::validate_conditions(const model& mod) const
                 break;
 
             case face_type::BOUNDARY:
-                if (not validate_boundary_condition(mod, bd.gmsh_entity) )
+                if (not validate_boundary_condition(mod, bd.material_tag()) )
                     success = false;
                 break;
 
             case face_type::INTERFACE:
-                if (not validate_interface_condition(mod, bd.gmsh_entity) )
+                if (not validate_interface_condition(mod, bd.material_tag()) )
                     success = false;
                 break;
 

src/maxwell_solver.cpp

@@ -221,6 +221,9 @@ void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
         do_sources(mod, state, mpl);
         
 
+        int rank;
+        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+        if (rank == 0)
             if ( (silo_output_rate != 0) and (i%silo_output_rate == 0) )
                 export_fields_to_silo(mod, state, mpl);
 
@@ -306,12 +309,26 @@ int main(int argc, const char *argv[])
     gmsh::open( mpl.sim_gmshmodel() );
 
     model mod( mpl.sim_geomorder(), mpl.sim_approxorder() );
-    
+    /*
     auto [scaled, sf] = mpl.mesh_scalefactor();
     if (scaled)
         mod.build(sf);
     else
         mod.build();
+    */
+
+    if (rank == 0)
+    {
+        mod.generate_mesh();
+        if (size > 1)
+            mod.partition_mesh(size);
+
+        mod.populate_from_gmsh();
+    }
+    else
+    {
+        mod.populate_from_gmsh();
+    }
 
     if ( not mpl.validate_model_params(mod) )
     {

src/silo_io.cpp

@@ -54,6 +54,9 @@ silo::gmsh_get_elements()
     gmsh::model::getEntities(entities, 3);
     for (auto [dim, tag] : entities)
     {
+        if (tag != 4 and tag != 5)
+            continue;
+
         std::vector<int> elemTypes;
         gmsh::model::mesh::getElementTypes(elemTypes, dim, tag);
         for (auto& elemType : elemTypes)