Overset intergrid boundary and hole search

src/Control/Inciter/InputDeck/Grammar.hpp

@@ -338,14 +338,16 @@ namespace grm {
         for (std::size_t i=0; i<n; ++i) bc[i] = b.count( i + 1 );
       };
 
-      // Generate BC mesh masks
+      // Generate mesh masks
       auto& conf = stack.template get< tag::param, eq, tag::mesh >();
       mask( conf.template get< tag::bcdir >() );
       mask( conf.template get< tag::bcsym >() );
       mask( conf.template get< tag::bcinlet >() );
       mask( conf.template get< tag::bcoutlet >() );
       mask( conf.template get< tag::bcfarfield >() );
       mask( conf.template get< tag::bcextrapolate >() );
+      auto& ig = stack.template get<tag::param, eq, tag::intergrid_boundary>();
+      mask( ig.template get< tag::mesh >() );
     }
   };
 
@@ -1454,6 +1456,32 @@ namespace deck {
                                         tk::grm::check_vector,
                                         eq, tag::bc, param > > > {};
 
+  //! Integrid boundary block
+  template< class eq >
+  struct intergrid_boundary :
+         pegtl::if_must<
+           tk::grm::readkw<typename use<kw::intergrid_boundary>::pegtl_string>,
+           tk::grm::start_vector< tag::param, eq, tag::intergrid_boundary,
+                                  tag::sideset >,
+           tk::grm::block<
+             use< kw::end >,
+             tk::grm::parameter_vector< use,
+                                        use< kw::mesh >,
+                                        tk::grm::Store_back,
+                                        tk::grm::noop,
+                                        tk::grm::check_vector,
+                                        eq,
+                                        tag::intergrid_boundary,
+                                        tag::mesh >,
+             tk::grm::parameter_vector< use,
+                                        use< kw::sideset >,
+                                        tk::grm::Store_back_back,
+                                        tk::grm::noop,
+                                        tk::grm::check_vector,
+                                        eq,
+                                        tag::intergrid_boundary,
+                                        tag::sideset > > > {};
+
   //! edgelist ... end block
   struct edgelist :
          tk::grm::vector< use< kw::amr_edgelist >,
@@ -1758,6 +1786,7 @@ namespace deck {
                            farfield_bc< kw::bc_farfield,
                                         tag::compflow,
                                         tag::bcfarfield >,
+                           intergrid_boundary< tag::compflow >,
                            bc< kw::bc_extrapolate, tag::compflow,
                                tag::bcextrapolate >,
                            timedep_bc< tag::compflow >

src/Control/Inciter/InputDeck/InputDeck.hpp

@@ -128,6 +128,7 @@ class InputDeck : public tk::TaggedTuple< InputDeckMembers > {
                                  , kw::diagnostics
                                  , kw::history_output
                                  , kw::mesh
+                                 , kw::intergrid_boundary
                                  , kw::filename
                                  , kw::location
                                  , kw::orientation

src/Control/Inciter/Types.hpp

@@ -315,14 +315,20 @@ using bc = tk::TaggedTuple< brigand::list<
 
 //! Mesh assignment and configuration
 using mesh = tk::TaggedTuple< brigand::list<
-    tag::filename,      std::vector< std::string >
-  , tag::velocity,      std::vector< kw::velocity::info::expect::type >
-  , tag::bcdir,         std::vector< std::size_t >
-  , tag::bcsym,         std::vector< std::size_t >
-  , tag::bcinlet,       std::vector< std::size_t >
-  , tag::bcoutlet,      std::vector< std::size_t >
-  , tag::bcfarfield,    std::vector< std::size_t >
-  , tag::bcextrapolate, std::vector< std::size_t >
+    tag::filename,           std::vector< std::string >
+  , tag::velocity,           std::vector< kw::velocity::info::expect::type >
+  , tag::bcdir,              std::vector< std::size_t >
+  , tag::bcsym,              std::vector< std::size_t >
+  , tag::bcinlet,            std::vector< std::size_t >
+  , tag::bcoutlet,           std::vector< std::size_t >
+  , tag::bcfarfield,         std::vector< std::size_t >
+  , tag::bcextrapolate,      std::vector< std::size_t >
+> >;
+
+//! Integrid boundary configuration
+using intergrid_boundary = tk::TaggedTuple< brigand::list<
+    tag::mesh,    std::vector< std::size_t >
+  , tag::sideset, std::vector< std::vector< std::size_t > >
 > >;
 
 //! Stagnation points parameters storage
@@ -417,6 +423,7 @@ using CompFlowPDEParameters = tk::TaggedTuple< brigand::list<
   , tag::farfield_pressure, kw::pressure::info::expect::type
   , tag::farfield_density,  kw::density::info::expect::type
   , tag::farfield_velocity, std::vector< kw::velocity::info::expect::type >
+  , tag::intergrid_boundary, intergrid_boundary
   , tag::bc,            bc
   , tag::bctimedep,     std::vector< time_dependent_bc >
   , tag::sponge,        SpongeParameters

src/Control/Keywords.hpp

@@ -2848,6 +2848,24 @@ struct cweight_info {
 };
 using cweight = keyword< cweight_info, TAOCPP_PEGTL_STRING("cweight") >;
 
+struct intergrid_boundary_info {
+  static std::string name() { return "intergrid_boundary"; }
+  static std::string shortDescription() { return
+    "Designate an side set an intergrid boundary";
+  }
+  static std::string longDescription() { return
+    R"(This keyword is used to select an side set of a mesh to be used as an
+       intergrid boundary through which solutions on multiple meshes
+       interact.)";
+  }
+  struct expect {
+    using type = std::string;
+    static std::string description() { return "strings"; }
+  };
+};
+using intergrid_boundary =
+  keyword< intergrid_boundary_info, TAOCPP_PEGTL_STRING("intergrid_boundary") >;
+
 struct sideset_info {
   static std::string name() { return "sideset"; }
   static std::string shortDescription() { return

src/Control/Tags.hpp

@@ -253,6 +253,8 @@ struct init_time { static std::string name() { return "init_time"; } };
 struct front_width { static std::string name() { return "front_width"; } };
 struct radius { static std::string name() { return "radius"; } };
 struct sideset { static std::string name() { return "sideset"; } };
+struct intergrid_boundary { static std::string name()
+ { return "intergrid_boundary"; } };
 struct sponge { static std::string name() { return "sponge"; } };
 struct bcinlet { static std::string name() { return "bcinlet"; } };
 struct bcoutlet { static std::string name() { return "bcoutlet"; } };

src/Inciter/CMakeLists.txt

@@ -20,6 +20,7 @@ add_library(Inciter
             FluxCorrector.cpp
             DistFCT.cpp
             DiagReducer.cpp
+            HoleReducer.cpp
             NodeDiagnostics.cpp
             ElemDiagnostics.cpp
             FieldOutput.cpp

src/Inciter/Discretization.cpp

@@ -24,11 +24,13 @@
 #include "QuinoaBuildConfig.hpp"
 #include "ConjugateGradients.hpp"
 #include "ALE.hpp"
+#include "HoleReducer.hpp"
 
 #include "M2MTransfer.hpp"
 
 namespace inciter {
 
+static CkReduction::reducerType HoleMerger;
 static CkReduction::reducerType PDFMerger;
 extern ctr::InputDeck g_inputdeck;
 extern ctr::InputDeck g_inputdeck_defaults;
@@ -98,7 +100,7 @@ Discretization::Discretization(
   m_meshvel( 0, 3 ),
   m_meshvel_converged( true ),
   m_bface( bface ),
-  m_triinpoel( triinpoel ),
+  m_triinpoel( tk::remap( triinpoel, m_lid ) ),
   m_elemblockid( elemblockid )
 // *****************************************************************************
 //  Constructor
@@ -456,6 +458,98 @@ Discretization::bndel() const
   return e;
 }
 
+void
+Discretization::setupHoles( CkCallback cont )
+// *****************************************************************************
+// Communicate holes to background mesh
+//! \param[in] cont Callback to continue with after holes have been aggregated
+// *****************************************************************************
+{
+  m_holcont = cont;
+
+  std::unordered_map< std::size_t, std::vector< tk::real > > hol;
+
+  if (m_meshid == 0) {
+    tk::destroy( m_hol );
+  }
+  else {
+    const auto& ib = g_inputdeck.get< tag::param, tag::compflow,
+                                      tag::intergrid_boundary >();
+    if (ib.get< tag::mesh >()[ m_meshid ]) {    // own meshid only
+      std::size_t nhol = 0;
+      for (const auto& ibs : ib.get< tag::sideset >()) {  // each hole
+        std::unordered_set< std::size_t > is;
+        for (auto s : ibs) is.insert( s );  // collect sidesets of hole
+        if (!is.empty()) {
+          const auto& x = m_coord[0];
+          const auto& y = m_coord[1];
+          const auto& z = m_coord[2];
+          auto& h = hol[ nhol++ ];
+          for (const auto& [ setid, faceids ] : m_bface) {
+            if ( is.count(setid) ) {
+              for (auto f : faceids) {
+                const auto t = m_triinpoel.data() + f*3;
+                h.push_back( x[t[0]] );
+                h.push_back( y[t[0]] );
+                h.push_back( z[t[0]] );
+                h.push_back( x[t[1]] );
+                h.push_back( y[t[1]] );
+                h.push_back( z[t[1]] );
+                h.push_back( x[t[2]] );
+                h.push_back( y[t[2]] );
+                h.push_back( z[t[2]] );
+              }
+            }
+          }
+        }
+      }
+    }
+
+    // all overset meshes send their hole-parts to the bg mesh for assembly
+    // (allreduce to all partitions of mesh 0)
+    for (std::size_t m=1; m<m_disc.size(); ++m) {
+      auto stream = serialize( m_meshid, hol );
+      contribute( stream.first, stream.second.get(), HoleMerger,
+       CkCallback(CkIndex_Discretization::aggregateHoles(nullptr), m_disc[0]) );
+    }
+  }
+}
+
+void
+Discretization::aggregateHoles( CkReductionMsg* msg )
+// *****************************************************************************
+//  Receive hole data from other meshes
+//! \param[in] msg Aggregated hole data
+//! \note Only bg (mesh 0) is supposed to call this.
+// *****************************************************************************
+{
+  std::size_t mid;
+  std::unordered_map< std::size_t, std::vector< tk::real > > inhol;
+
+  PUP::fromMem creator( msg->getData() );
+  creator | mid;
+  creator | inhol;
+  delete msg;
+
+  for (auto&& [hid,data] : inhol) {
+    tk::concat( std::move(data), m_hol[mid][hid] );
+  }
+
+  // back to sender overset mesh so it can continue (enough to this partition)
+  m_disc[mid][thisIndex].holComplete();
+  // bg mesh also complete
+  m_holcont.send();
+}
+
+void
+Discretization::holComplete()
+// *****************************************************************************
+// Hole communication complete
+// *****************************************************************************
+{
+  m_holcont.send();
+}
+
 void
 Discretization::resizePostAMR(
   const tk::UnsMesh::Chunk& chunk,
@@ -524,6 +618,7 @@ Discretization::registerReducers()
 // *****************************************************************************
 {
   PDFMerger = CkReduction::addReducer( tk::mergeUniPDFs );
+  HoleMerger = CkReduction::addReducer( mergeHole );
 }
 
 tk::UnsMesh::Coords

src/Inciter/Discretization.hpp

@@ -176,6 +176,9 @@ class Discretization : public CBase_Discretization {
       const std::unordered_map< std::size_t, std::set< std::size_t > >& nodeblk,
       std::size_t nuserblk );
 
+    //! Communicate holes to background mesh
+    void setupHoles( CkCallback cont );
+
     /** @name Accessors */
     ///@{
     //! Coordinates accessor as const-ref
@@ -307,6 +310,11 @@ class Discretization : public CBase_Discretization {
 
     //! Edge communication map accessor as const-ref
     const tk::EdgeCommMap& EdgeCommMap() const { return m_edgeCommMap; }
+
+    //! Holes data accessor
+    const std::unordered_map< std::size_t,
+      std::unordered_map< std::size_t, std::vector< tk::real > > >& hol() const
+    { return m_hol; }
     //@}
 
     //! Set time step size
@@ -409,6 +417,12 @@ class Discretization : public CBase_Discretization {
     //! Find elements along our mesh chunk boundary
     std::vector< std::size_t > bndel() const;
 
+    //! Receive hole data
+    void aggregateHoles( CkReductionMsg* msg );
+
+    //! Hole communication complete
+    void holComplete();
+
     //! Decide if field output iteration count interval is hit
     bool fielditer() const;
 
@@ -504,6 +518,8 @@ class Discretization : public CBase_Discretization {
 
     //! Mesh ID
     std::size_t m_meshid;
+    //! Callback to continue with after hole communication
+    CkCallback m_holcont;
     //! \brief Charm++ callback of the function to call after a mesh-to-mesh
     //!   solution transfer (to-and-fro) is complete
     CkCallback m_transfer_complete;
@@ -645,6 +661,12 @@ class Discretization : public CBase_Discretization {
     std::vector< std::size_t > m_triinpoel;
     //! Local tet ids associated with mesh block ids
     std::unordered_map< std::size_t, std::set< std::size_t > > m_elemblockid;
+    //! Holes tessellation
+    //! \details: Value: list of triangle node coordinates
+    //!    (x1,y1,z1,x2,y2,z2,x3,y3,z3), ..., inner key, hole id, outer key
+    //!    meshid
+    std::unordered_map< std::size_t,
+      std::unordered_map< std::size_t, std::vector< tk::real > > > m_hol;
 
     //! Generate the Bid data-structure based on the node communication-map
     std::unordered_map< std::size_t, std::size_t > genBid();