Simplified boundary handling.

include/libgmshdg/gmsh_io.h

@@ -132,25 +132,6 @@ operator<<(std::ostream& os, const comm_descriptor& cd)
 
 #define IMPLIES(a, b) ((not(a)) or (b))
 #define IFF(a,b) (IMPLIES(a,b) and IMPLIES(b,a))
-struct surface_identifier
-{
-    int     surface_entity;
-    int     adjacent_volume_entity;
-
-    bool operator<(const surface_identifier& other) const
-    {
-        return (surface_entity < other.surface_entity) or
-                    ((surface_entity == other.surface_entity) and
-                    (adjacent_volume_entity < other.adjacent_volume_entity));
-    }
-};
-
-inline std::ostream&
-operator<<(std::ostream& os, const surface_identifier& si)
-{
-    os << "(" << si.surface_entity << ", " << si.adjacent_volume_entity << ")";
-    return os;
-}
 
 struct surface_descriptor
 {
@@ -244,8 +225,8 @@ class model
     std::vector<boundary_descriptor>            bnd_descriptors;
     element_connectivity<element_key>           conn;               
 
-    std::map<surface_identifier, surface_descriptor>     boundaries;
-    std::map<surface_identifier, surface_descriptor>     interfaces;
+    std::map<int, surface_descriptor>           boundaries;
+    std::map<int, surface_descriptor>           interfaces;
 
 #ifdef USE_MPI
     interprocess_connectivity<element_key>      ipconn;
@@ -284,7 +265,6 @@ class model
     void    populate_from_gmsh_rank0(void);
 
 #ifdef USE_MPI
-    bool    is_interprocess_boundary(int) const;
     void    populate_from_gmsh_rankN(void);
     void    import_gmsh_entities_rankN(void);
     int     my_partition(void) const;
@@ -315,13 +295,15 @@ public:
 
     void map_boundaries_new(void);
 
-    const std::map<surface_identifier, surface_descriptor>&
+    const std::map<int, surface_descriptor>&
     get_boundaries() const;
 
-    const std::map<surface_identifier, surface_descriptor>&
+    const std::map<int, surface_descriptor>&
     get_interfaces() const;
 
 #ifdef USE_MPI
+    bool    is_interprocess_boundary(int) const;
+
     const std::map<int, comm_descriptor>& comm_descriptors(void) const {
         return ipc_boundary_comm_table;
     }

share/validation/params_twomat.lua

@@ -13,6 +13,9 @@ postpro.cycle_print_rate = 100                  -- console print rate
 
 postpro["E"].mode = "nodal"
 
+debug = {};
+debug.dump_cell_ranks = true;
+
 function setup_materials()
     local bulk = volume_physical_group_entities(1);
     for i,v in ipairs(bulk) do

share/validation/twomat.geo

@@ -9,4 +9,4 @@ Physical Surface("pmc", 11) = {3,4,8,9};
 Physical Surface("pec", 12) = {5,6,10,11};
 Physical Surface("abc", 13) = {7};
 
-MeshSize{:} = 0.075;
+MeshSize{:} = 0.05;

src/libgmshdg/gmsh_io.cpp

@@ -528,37 +528,27 @@ model::map_boundaries_new(void)
             if ( (itor == bfk.end()) or not ((*itor).first == fk) )
                 continue;
             
-            surface_identifier si;
-            si.surface_entity = (*itor).second;
-            si.adjacent_volume_entity = iE;
+            int si = (*itor).second;
 
-            using sm_t = std::map<surface_identifier, surface_descriptor>;
+            using sm_t = std::map<int, surface_descriptor>;
 
             auto add_face = [&](sm_t& m, size_t se, size_t ave) -> void {
                 if ( m.find(si) == m.end() )
                 {
                     m[si].surface_entity = se;
                     m[si].adjacent_volume_entity = ave;
-#ifdef USE_MPI
-                    if (num_partitions > 1)
-                    {
-                        auto sitor = surface_to_parent_map.find( se );
-                        if ( sitor != surface_to_parent_map.end() )
-                            m[si].parent_entity = surface_to_parent_map.at( se );
-                    }
-#endif /* USE_MPI */
                 }
                 m[si].add_face(iF, fbase+iF);
             };
 
-            switch ( face_type(si.surface_entity, fk) )
+            switch ( face_type(si, fk) )
             {
                 case face_type::BOUNDARY:
-                    add_face(boundaries, si.surface_entity, si.adjacent_volume_entity);
+                    add_face(boundaries, si, iE);
                     break;
 
                 case face_type::INTERFACE:
-                    add_face(interfaces, si.surface_entity, si.adjacent_volume_entity);
+                    add_face(interfaces, si, iE);
                     break;
 
                 default:
@@ -584,13 +574,13 @@ model::map_boundaries_new(void)
     }
 }
 
-const std::map<surface_identifier, surface_descriptor>&
+const std::map<int, surface_descriptor>&
 model::get_boundaries() const
 {
     return boundaries;
 }
 
-const std::map<surface_identifier, surface_descriptor>&
+const std::map<int, surface_descriptor>&
 model::get_interfaces() const
 {
     return interfaces;
@@ -661,14 +651,7 @@ model::map_boundaries(void)
             }
             /* and store the gmsh tag in the descriptor. */
             bd.gmsh_entity = (*itor).second;
-#ifdef USE_MPI
-            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 );
-            }
-#endif /* USE_MPI */
+
             /* Finally, put the descriptor in the global array of faces. */
             bnd_descriptors.at(fbase + iF) = bd;
         }

src/maxwell/maxwell_common.cpp

@@ -23,7 +23,7 @@ namespace maxwell {
  * coefficients 0, 1 and 2 to select the appropriate boundary condition.
  * See the jump kernels for details about the 0, 1, and 2. */
 void
-init_boundary_conditions(const model& mod, const parameter_loader& mpl,
+init_boundary_conditions_old(const model& mod, const parameter_loader& mpl,
     vecxd& bcc)
 {
     auto bds = mod.boundary_descriptors();
@@ -82,6 +82,56 @@ init_boundary_conditions(const model& mod, const parameter_loader& mpl,
     }
 }
 
+/* Take a model and a parameter loader and return a vector with the
+ * coefficients 0, 1 and 2 to select the appropriate boundary condition.
+ * See the jump kernels for details about the 0, 1, and 2. */
+void
+init_boundary_conditions(const model& mod, const parameter_loader& mpl,
+    vecxd& bcc)
+{
+    bcc = vecxd::Ones( mod.num_fluxes() );
+
+    auto& bds = mod.get_boundaries();
+
+    for (auto& [si, sd] : bds)
+    {
+        const auto& e = mod.entity_at(sd.adjacent_volume_entity);
+        auto mtag = sd.gmsh_tag();
+
+        for (auto& [local, global] : sd.face_numbers)
+        {
+            double bc_coeff = DIRICHLET;
+
+            switch ( mpl.boundary_type(mtag) )
+            {
+                case boundary_condition::UNSPECIFIED:
+                case boundary_condition::PEC:
+                case boundary_condition::E_FIELD:
+                    bc_coeff = DIRICHLET;
+                    break;
+
+                case boundary_condition::IMPEDANCE:
+                case boundary_condition::PLANE_WAVE_E:
+                case boundary_condition::PLANE_WAVE_H:
+                    bc_coeff = ROBIN;
+                    break;
+
+                case boundary_condition::PMC:
+                case boundary_condition::SURFACE_CURRENT:
+                case boundary_condition::H_FIELD:
+                    bc_coeff = NEUMANN;
+                    break;
+            }
+
+            auto& pf = e.face(local);
+            auto& rf = e.face_refelem(pf);
+            auto num_bf = rf.num_basis_functions();
+            for (size_t i = 0; i < num_bf; i++)
+                bcc(global*num_bf + i) = bc_coeff;
+        }
+    }
+}
+
 /* Take a model and a parameter loader and return the four vectors with
  * the Lax-Milgram flux coefficients */
 void