Interprocess connectivity.

include/connectivity.h

@@ -6,8 +6,13 @@
 #include <map>
 #include <set>
 
+#include "point.hpp"
 #include "sgr.hpp"
 
+#ifdef USE_MPI
+#include "gmsh_mpi.h"
+#endif /* USE_MPI */
+
 /***************************************************************************/
 struct neighbour_descriptor
 {
@@ -22,6 +27,25 @@ struct neighbour_descriptor
 
 std::ostream& operator<<(std::ostream&, const neighbour_descriptor&);
 
+struct neighbour_descriptor_interprocess
+{
+    size_t      iE;         /* Entity offset in data structures */
+    size_t      iT;         /* Element offset in entity in data structures */
+    size_t      via_iF;     /* Element face index (i.e. 0,1,2 for triangles) */
+
+    int         parent_tag; /* Tag of the parent entity */
+    int         tag;        /* Tag of the actual entity */
+    int         partition;  /* Partition on which the element resides */
+
+    point_3d    bar;        /* Barycenter, for materials. */
+
+    neighbour_descriptor_interprocess();
+    bool operator<(const neighbour_descriptor_interprocess&) const;
+
+};
+
+std::ostream& operator<<(std::ostream&, const neighbour_descriptor_interprocess&);
+
 /***************************************************************************/
 struct cell_descriptor
 {
@@ -235,3 +259,75 @@ public:
     }
 };
 
+
+
+
+/***************************************************************************/
+#ifdef USE_MPI
+template<typename FaceKey>
+class interprocess_connectivity
+{
+public:
+    using ConnData  = std::pair<FaceKey, neighbour_descriptor_interprocess>;
+    using NeighPair = neighbour_pair<neighbour_descriptor_interprocess>;
+
+private:
+    std::map<FaceKey, NeighPair>                    face_cell_adj;
+
+public:
+    interprocess_connectivity()
+    {}
+
+    void
+    connect(const FaceKey& fk, const neighbour_descriptor_interprocess& nd)
+    {
+        /* Insert the information that one the two neighbours around
+         * face `e' is `nd'. */
+        auto& nr = face_cell_adj[fk];
+        nr.add_neighbour(nd);
+    }
+
+    /* Find the neighbour of the element (iE,iT) via face fk */
+    std::pair<neighbour_descriptor_interprocess, bool>
+    neighbour_via(size_t iE, size_t iT, const FaceKey& fk) const
+    {
+        auto itor = face_cell_adj.find(fk);
+        if (itor == face_cell_adj.end())
+            return std::make_pair(neighbour_descriptor_interprocess(), false);
+
+        auto nr = (*itor).second;
+        assert( nr.num_neighbours() == 2 );
+
+        auto fst = nr.first();
+        auto snd = nr.second();
+
+        assert( (fst.iE == iE and fst.iT == iT) or
+                (snd.iE == iE and snd.iT == iT)
+              );
+
+        if ( fst.iE == iE and fst.iT == iT)
+            return std::make_pair(snd, true);
+
+        return std::make_pair(fst, true);
+    }
+
+    void clear(void)
+    {
+        face_cell_adj.clear();
+    }
+
+    void dump(void)
+    {
+        for (auto& [key, pair] : face_cell_adj)
+        {
+            std::cout << key << ": " << pair.first() << " - " << pair.second() << std::endl;
+        }
+    }
+
+    void mpi_bcast(int root, MPI_Comm comm)
+    {
+        priv_MPI_Bcast(face_cell_adj, root, comm);
+    }
+};
+
+#endif /* USE_MPI */
\ No newline at end of file

include/gmsh_io.h

@@ -111,6 +111,8 @@ class model
     element_connectivity<element_key>           conn;
 
 #ifdef USE_MPI
+    interprocess_connectivity<element_key>      ipconn;
+
     /* Map from boundary tag to all the partitions it belongs to */
     std::map<int, std::vector<int>>             comm_map;       
 
@@ -145,6 +147,7 @@ class model
     void    import_gmsh_entities_rankN(void);
     int     my_partition(void) const;
     void    make_partition_to_entities_map(void);
+    void    map_interprocess_boundaries();
 #endif /* USE_MPI */
 
     void    make_boundary_to_faces_map(void);
@@ -206,6 +209,12 @@ public:
 
     entity&         entity_world_lookup(int);
     const entity&   entity_world_lookup(int) const;
+
+    const interprocess_connectivity<element_key>& ip_connectivity() const
+    {
+        return ipconn;
+    }
+
 #endif /* USE_MPI */
 
     const element_connectivity<element_key>& connectivity() const

src/connectivity.cpp

@@ -24,6 +24,27 @@ std::ostream& operator<<(std::ostream& os, const neighbour_descriptor& nd)
     return os;
 }
 
+
+
+
+neighbour_descriptor_interprocess::neighbour_descriptor_interprocess()
+{}
+
+bool
+neighbour_descriptor_interprocess::operator<(const neighbour_descriptor_interprocess& other) const
+{
+    std::array<size_t, 3>   arr_me{iE, iT, via_iF};
+    std::array<size_t, 3>   arr_other{other.iE, other.iT, other.via_iF};
+    return arr_me < arr_other;
+}
+
+std::ostream& operator<<(std::ostream& os, const neighbour_descriptor_interprocess& nd)
+{
+    os << "(" << nd.iE << ", " << nd.iT << ", " << nd.via_iF << ", PTAG: " << nd.parent_tag;
+    os << ", TAG: " << nd.tag << ", PART: " << nd.partition << ")";
+    return os;
+}
+
 /***************************************************************************/
 cell_descriptor::cell_descriptor()
 {}

src/gmsh_io.cpp

@@ -93,7 +93,9 @@ model::make_boundary_to_faces_map(void)
              * used to determine who communicates with who. */
             std::vector<int> parts;
             gm::getPartitions(dim, tag, parts);
-            comm_map[tag] = parts;
+            if (parts.size() == 2)
+                comm_map[tag] = parts;
+
             if (pdim == dim)
                 surface_to_parent_map[tag] = ptag;
         }
@@ -157,10 +159,6 @@ 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++)
     {
         for (size_t iF = 0; iF < e.num_faces_per_elem(); iF++)
@@ -368,7 +366,7 @@ bool
 model::is_interprocess_boundary(int tag)
 {
     if (num_partitions > 1)
-        return comm_map.at(tag).size() == 2;
+        return comm_map.find(tag) != comm_map.end();
     
     return false;
 }
@@ -468,6 +466,73 @@ model::map_boundaries(void)
 #endif
 }
 
+#ifdef USE_MPI
+void
+model::map_interprocess_boundaries(void)
+{
+    /* Make a vector mapping element_key to entity tag */
+    using bfk_t = std::pair<element_key, int>;
+    std::vector<bfk_t> bfk;
+    for (auto& [tag, keys] : boundary_map)
+    {
+        if ( not is_interprocess_boundary(tag) )
+            continue;
+        
+        for (auto& k : keys)
+            bfk.push_back( std::make_pair(k, tag) );
+    }
+    /* Sort it for lookups */
+    auto comp = [](const bfk_t& a, const bfk_t& b) -> bool {
+        return a.first < b.first;
+    };
+    std::sort(bfk.begin(), bfk.end(), comp);
+
+    auto lbcomp = [](const bfk_t& a, const element_key& b) -> bool {
+        return a.first < b;
+    };
+
+    auto connect = [&](const entity& e, int part) {
+        for (size_t iT = 0; iT < e.num_cells(); iT++)
+        {
+            auto nfpe = e.num_faces_per_elem();
+            for (size_t iF = 0; iF < nfpe; iF++)
+            {
+                auto fnum = iT*nfpe + iF;
+                element_key fk(e.face(fnum));
+
+                auto itor = std::lower_bound(bfk.begin(), bfk.end(), fk, lbcomp);
+                if ( (itor == bfk.end()) or not ((*itor).first == fk) )
+                    continue;
+            
+                neighbour_descriptor_interprocess ndi;
+                ndi.iE = e.number();
+                ndi.iT = iT;
+                ndi.via_iF = iF;
+                ndi.parent_tag = e.material_tag();
+                ndi.tag = e.gmsh_tag();
+                ndi.partition = part;
+                ndi.bar = e.cell(iT).barycenter();
+                ipconn.connect(fk, ndi);
+            }
+        }
+    };
+
+    for (auto& e : entities)
+    {
+        auto part = my_partition();
+        assert(part == 1);
+        connect(e, part);
+    }
+
+    for (auto& [part, res] : remote_entities)
+    {
+        assert(part != 1);
+        for (auto& e : res)
+            connect(e, part);
+    }
+}
+#endif
+
 void
 model::populate_from_gmsh_rank0()
 {
@@ -512,6 +577,11 @@ model::populate_from_gmsh_rank0()
     entities.clear();
     import_gmsh_entities_rank0();
     map_boundaries(); /* This must happen after import_gmsh_entities(). */
+    
+#ifdef USE_MPI
+    map_interprocess_boundaries();
+    ipconn.mpi_bcast(0, MPI_COMM_WORLD);
+#endif /* USE_MPI */
 }
 
 #ifdef USE_MPI
@@ -539,6 +609,7 @@ model::populate_from_gmsh_rankN()
     entities.clear();
     import_gmsh_entities_rankN();
     map_boundaries();
+    ipconn.mpi_bcast(0, MPI_COMM_WORLD);
 }
 #endif /* USE_MPI */
 

src/maxwell_common.cpp

@@ -2,9 +2,10 @@
 #include "maxwell_interface.h"
 #include "maxwell_common.h"
 
-#define DIRICHLET   2.0
-#define ROBIN       1.0
-#define NEUMANN     0.0
+#define DIRICHLET       2.0
+#define ROBIN           1.0
+#define INTERPROCESS    1.0
+#define NEUMANN         0.0
 
 namespace maxwell {
 
@@ -29,25 +30,32 @@ init_boundary_conditions(const model& mod, const parameter_loader& mpl,
                 continue;
 
             double bc_coeff = DIRICHLET;
-            switch ( mpl.boundary_type(bd.material_tag()) )
+            if (bd.type == face_type::INTERPROCESS_BOUNDARY)
             {
-                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;
+                bc_coeff = INTERPROCESS;
+            }
+            else
+            {
+                switch ( mpl.boundary_type(bd.material_tag()) )
+                {
+                    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(iF);
@@ -73,6 +81,9 @@ init_lax_milgram(const model& mod, const parameter_loader& mpl,
     pbH = vecxd::Zero( mod.num_fluxes() );
 
     auto& conn = mod.connectivity();
+#ifdef USE_MPI
+    auto& ipconn = mod.ip_connectivity();
+#endif /* USE_MPI */
     for (auto& e : mod)
     {
         for (size_t iT = 0; iT < e.num_cells(); iT++)
@@ -120,16 +131,39 @@ init_lax_milgram(const model& mod, const parameter_loader& mpl,
                         paH(glob_ofs+iD) = aH;
                         pbH(glob_ofs+iD) = bH;
                     }
+                    continue;
                 }
-                else
+#ifdef USE_MPI
+                auto [ip_neigh, has_ip_neigh] = ipconn.neighbour_via(e.number(), iT, fk);
+                if (has_ip_neigh)
                 {
+                    auto ne_eps = mpl.epsilon(ip_neigh.parent_tag, ip_neigh.bar);
+                    auto ne_mu = mpl.mu(ip_neigh.parent_tag, ip_neigh.bar);
+
+                    auto ne_Z = std::sqrt(ne_mu/ne_eps);
+                    auto ne_Y = 1./ne_Z;
+
+                    auto aE = ne_Z/(my_Z + ne_Z);
+                    auto bE = 1./(my_Z + ne_Z);
+                    auto aH = ne_Y/(my_Y + ne_Y);
+                    auto bH = 1./(my_Y + ne_Y);
+
                     for (size_t iD = 0; iD < rf.num_basis_functions(); iD++)
                     {
-                        paE(glob_ofs+iD) = 0.5;
-                        pbE(glob_ofs+iD) = 0.5/my_Z;
-                        paH(glob_ofs+iD) = 0.5;
-                        pbH(glob_ofs+iD) = 0.5/my_Y;
+                        paE(glob_ofs+iD) = aE;
+                        pbE(glob_ofs+iD) = bE;
+                        paH(glob_ofs+iD) = aH;
+                        pbH(glob_ofs+iD) = bH;
                     }
+                    continue;
+                }
+#endif /* USE_MPI */
+                for (size_t iD = 0; iD < rf.num_basis_functions(); iD++)
+                {
+                    paE(glob_ofs+iD) = 0.5;
+                    pbE(glob_ofs+iD) = 0.5/my_Z;
+                    paH(glob_ofs+iD) = 0.5;
+                    pbH(glob_ofs+iD) = 0.5/my_Y;
                 }
             }
         }