Pass only retained nodes to MueLu

Author: amgebauer

src/core/fem/src/discretization/4C_fem_discretization.cpp

@@ -563,12 +563,13 @@ void Core::FE::Discretization::replace_dof_set(
 
 /*----------------------------------------------------------------------*
  *----------------------------------------------------------------------*/
-std::map<int, std::vector<int>>* Core::FE::Discretization::get_all_pbc_coupled_col_nodes()
+const std::map<int, std::vector<int>>* Core::FE::Discretization::get_all_pbc_coupled_col_nodes()
+    const
 {
   // check for pbcs
   for (int nds = 0; nds < num_dof_sets(); nds++)
   {
-    std::shared_ptr<Core::DOFSets::PBCDofSet> pbcdofset =
+    std::shared_ptr<const Core::DOFSets::PBCDofSet> pbcdofset =
         std::dynamic_pointer_cast<Core::DOFSets::PBCDofSet>(dofsets_[nds]);
 
     if (pbcdofset != nullptr)
@@ -584,8 +585,8 @@ std::map<int, std::vector<int>>* Core::FE::Discretization::get_all_pbc_coupled_c
 
 /*----------------------------------------------------------------------*
  *----------------------------------------------------------------------*/
-std::shared_ptr<std::map<int, int>>
-Core::FE::Discretization::get_pbc_slave_to_master_node_connectivity()
+std::shared_ptr<const std::map<int, int>>
+Core::FE::Discretization::get_pbc_slave_to_master_node_connectivity() const
 {
   // check for pbcs
   for (int nds = 0; nds < num_dof_sets(); nds++)

src/core/fem/src/discretization/4C_fem_discretization.hpp

@@ -622,12 +622,12 @@ namespace Core::FE
     /*!
     \brief Get master to slave coupling in case of periodic boundary conditions
     */
-    std::map<int, std::vector<int>>* get_all_pbc_coupled_col_nodes();
+    const std::map<int, std::vector<int>>* get_all_pbc_coupled_col_nodes() const;
 
     /*!
     \brief Get slave to master connectivity in case of periodic boundary conditions
     */
-    std::shared_ptr<std::map<int, int>> get_pbc_slave_to_master_node_connectivity();
+    std::shared_ptr<const std::map<int, int>> get_pbc_slave_to_master_node_connectivity() const;
 
     //@}
 

src/core/fem/src/discretization/4C_fem_discretization_faces.cpp

@@ -239,7 +239,7 @@ void Core::FE::DiscretizationFaces::build_faces(const bool verbose)
   std::map<std::vector<int>, std::shared_ptr<Core::Elements::Element>> faces;
 
   // get pbcs
-  std::map<int, std::vector<int>>* col_pbcmapmastertoslave = get_all_pbc_coupled_col_nodes();
+  const std::map<int, std::vector<int>>* col_pbcmapmastertoslave = get_all_pbc_coupled_col_nodes();
 
   std::map<std::vector<int>, InternalFacesData>::iterator face_it;
   for (face_it = surfmapdata.begin(); face_it != surfmapdata.end(); ++face_it)
@@ -383,11 +383,11 @@ void Core::FE::DiscretizationFaces::build_faces(const bool verbose)
             for (std::size_t rr = 0; rr < mymasternodeids.size(); rr++)
             {
               // this master node has one slave node
-              if (((*col_pbcmapmastertoslave)[mymasternodeids[rr]]).size() == 1)
+              if (((*col_pbcmapmastertoslave).at(mymasternodeids[rr])).size() == 1)
               {
-                myslavenodeids.push_back(((*col_pbcmapmastertoslave)[mymasternodeids[rr]])[0]);
+                myslavenodeids.push_back(((*col_pbcmapmastertoslave).at(mymasternodeids[rr]))[0]);
                 local_pbcmapmastertoslave[mymasternodeids[rr]] =
-                    ((*col_pbcmapmastertoslave)[mymasternodeids[rr]])[0];
+                    ((*col_pbcmapmastertoslave).at(mymasternodeids[rr]))[0];
               }
               // this master node has several slave nodes
               // it is a corner or edge node of two or three pbc sets
@@ -443,7 +443,7 @@ void Core::FE::DiscretizationFaces::build_faces(const bool verbose)
                 //  this special case is marked by flag
                 bool three_sets_edge_node = false;
                 if (numpbcpairs == 3 and
-                    ((*col_pbcmapmastertoslave)[mymasternodeids[rr]]).size() == 3)
+                    ((*col_pbcmapmastertoslave).at(mymasternodeids[rr])).size() == 3)
                   three_sets_edge_node = true;
 
                 // pbc id of master face also for the slave
@@ -459,10 +459,10 @@ void Core::FE::DiscretizationFaces::build_faces(const bool verbose)
                 // check which node fulfills above conditions
                 int actslaveid = -999;
                 for (std::size_t islave = 0;
-                    islave < ((*col_pbcmapmastertoslave)[mymasternodeids[rr]]).size(); islave++)
+                    islave < ((*col_pbcmapmastertoslave).at(mymasternodeids[rr])).size(); islave++)
                 {
                   // get id
-                  actslaveid = ((*col_pbcmapmastertoslave)[mymasternodeids[rr]])[islave];
+                  actslaveid = ((*col_pbcmapmastertoslave).at(mymasternodeids[rr]))[islave];
 
                   // check first criterion -> (i)
                   if ((slavetopbcset[slavepbcid]).find(actslaveid) !=
@@ -547,7 +547,7 @@ void Core::FE::DiscretizationFaces::build_faces(const bool verbose)
 
                 // first, look for masters with more than one slave
                 // then, check whether node id (i.e. actnodeid) is contained in slave list
-                std::map<int, std::vector<int>>::iterator master_it;
+                std::map<int, std::vector<int>>::const_iterator master_it;
                 for (master_it = col_pbcmapmastertoslave->begin();
                     master_it != col_pbcmapmastertoslave->end(); master_it++)
                 {

src/core/fem/src/discretization/4C_fem_discretization_partition.cpp

@@ -569,7 +569,7 @@ void Core::FE::Discretization::extended_ghosting(const Core::LinAlg::Map& elecol
     {
       have_pbc = true;
       // fill content of pbcmap int std::map<int, std::set<int> > in preparation for gather_all
-      std::map<int, std::vector<int>>* tmp = pbcdofset->get_coupled_nodes();
+      const std::map<int, std::vector<int>>* tmp = pbcdofset->get_coupled_nodes();
       for (auto& [gid, nodes] : *tmp) pbcmap[gid].insert(nodes.begin(), nodes.end());
 
       // it is assumed that, if one pbc set is available, all other potential dofsets hold the same

src/core/fem/src/discretization/4C_fem_discretization_utils.cpp

@@ -7,6 +7,7 @@
 
 #include "4C_fem_discretization_utils.hpp"
 
+#include "4C_comm_mpi_utils.hpp"
 #include "4C_fem_discretization.hpp"
 #include "4C_fem_general_element.hpp"
 #include "4C_fem_general_node.hpp"
@@ -48,6 +49,57 @@ std::shared_ptr<Core::LinAlg::MultiVector<double>> Core::FE::extract_node_coordi
   return coordinates;
 }
 
+std::unique_ptr<Core::LinAlg::MultiVector<double>> Core::FE::extract_retained_node_coordinates(
+    const Discretization& discretization, const Core::LinAlg::Map& node_row_map)
+{
+  std::set<int> skipped_nodes = {};
+  if (discretization.get_pbc_slave_to_master_node_connectivity())
+  {
+    const auto& pbcconnectivity = *(discretization.get_pbc_slave_to_master_node_connectivity());
+    for (const auto& [slave, _] : pbcconnectivity)
+    {
+      skipped_nodes.insert(slave);
+    }
+  }
+  std::set<int> fully_redundant_skipped_nodes =
+      Core::Communication::all_reduce(skipped_nodes, node_row_map.get_comm());
+
+  std::vector<int> my_filtered_global_node_ids{};
+  my_filtered_global_node_ids.reserve(node_row_map.num_my_elements());
+  for (int lid = 0; lid < node_row_map.num_my_elements(); ++lid)
+  {
+    const int gid = node_row_map.gid(lid);
+
+    // filter pbc slave nodes
+    if (fully_redundant_skipped_nodes.contains(gid)) continue;
+    my_filtered_global_node_ids.emplace_back(gid);
+  }
+  Core::LinAlg::Map filtered_row_map(
+      node_row_map.num_global_elements() - static_cast<int>(fully_redundant_skipped_nodes.size()),
+      static_cast<int>(my_filtered_global_node_ids.size()), my_filtered_global_node_ids.data(), 0,
+      node_row_map.get_comm());
+
+  std::unique_ptr<Core::LinAlg::MultiVector<double>> coordinates =
+      std::make_unique<Core::LinAlg::MultiVector<double>>(filtered_row_map, 3, true);
+
+  for (int lid = 0; lid < filtered_row_map.num_my_elements(); ++lid)
+  {
+    const int gid = filtered_row_map.gid(lid);
+    if (!discretization.have_global_node(gid)) continue;
+
+    auto x = discretization.g_node(gid)->x();
+    for (size_t dim = 0; dim < 3; ++dim)
+    {
+      if (dim >= discretization.n_dim())
+        coordinates->replace_local_value(lid, dim, 0.0);
+      else
+        coordinates->replace_local_value(lid, dim, x[dim]);
+    }
+  }
+
+  return coordinates;
+}
+
 
 /*----------------------------------------------------------------------------*
  *----------------------------------------------------------------------------*/

src/core/fem/src/discretization/4C_fem_discretization_utils.hpp

@@ -60,6 +60,17 @@ namespace Core::FE
   std::shared_ptr<Core::LinAlg::MultiVector<double>> extract_node_coordinates(
       const Core::FE::Discretization& discretization, const Core::LinAlg::Map& node_row_map);
 
+  /**
+   * Similar to the other extract_node_coordinates function, but allows to specify a node row map
+   * and filters all eliminated nodes (e.g., in periodic boundary conditions).
+   *
+   * @note Eliminated nodes (e.g., when using periodic boundary conditions) are skipped.
+   *
+   * The map has to be a submap of the overall full node rowmap of this discretization.
+   */
+  std::unique_ptr<Core::LinAlg::MultiVector<double>> extract_retained_node_coordinates(
+      const Core::FE::Discretization& discretization, const Core::LinAlg::Map& node_row_map);
+
 
   /** \brief Evaluate the elements of the given discretization and fill the
    *         system matrix and vector

src/core/fem/src/dofset/4C_fem_dofset_pbc.hpp

@@ -70,10 +70,14 @@ namespace Core::DOFSets
     virtual void set_coupled_nodes(std::shared_ptr<std::map<int, std::vector<int>>> couplednodes);
 
     /// Get coupled nodes map (corresponding col format)
-    std::map<int, std::vector<int>>* get_coupled_nodes() { return perbndcouples_.get(); }
+    [[nodiscard]] const std::map<int, std::vector<int>>* get_coupled_nodes() const
+    {
+      return perbndcouples_.get();
+    }
 
     /// Get connectivity map between slave node and its master node
-    virtual std::shared_ptr<std::map<int, int>> get_slave_to_master_node_connectivity()
+    [[nodiscard]] virtual std::shared_ptr<const std::map<int, int>>
+    get_slave_to_master_node_connectivity() const
     {
       return perbnd_slavetomaster_;
     };

src/core/fem/src/general/utils/4C_fem_general_l2_projection.cpp

@@ -135,7 +135,7 @@ std::shared_ptr<Core::LinAlg::MultiVector<double>> Core::FE::compute_nodal_l2_pr
   // build inverse map from slave to master nodes
   std::map<int, int> slavetomastercolnodesmap;
 
-  std::map<int, std::vector<int>>* allcoupledcolnodes = dis.get_all_pbc_coupled_col_nodes();
+  const std::map<int, std::vector<int>>* allcoupledcolnodes = dis.get_all_pbc_coupled_col_nodes();
   if (allcoupledcolnodes)
   {
     for (auto [master_gid, slave_gids] : *allcoupledcolnodes)

src/core/fem/src/general/utils/4C_fem_general_utils_superconvergent_patch_recovery.cpp

@@ -41,7 +41,7 @@ std::shared_ptr<Core::LinAlg::MultiVector<double>> Core::FE::compute_superconver
   // handle pbcs if existing
   // build inverse map from slave to master nodes
   std::map<int, int> slavetomastercolnodesmap;
-  std::map<int, std::vector<int>>* allcoupledcolnodes = dis.get_all_pbc_coupled_col_nodes();
+  const std::map<int, std::vector<int>>* allcoupledcolnodes = dis.get_all_pbc_coupled_col_nodes();
 
   if (allcoupledcolnodes)
   {
@@ -432,8 +432,7 @@ std::shared_ptr<Core::LinAlg::MultiVector<double>> Core::FE::compute_superconver
         // this will result in off processor assembling (boundary node as ghost node)
         if (closestnode->owner() != myrank) continue;
 
-        std::map<int, std::vector<int>>::iterator masternode =
-            allcoupledcolnodes->find(closestnodeid);
+        auto masternode = allcoupledcolnodes->find(closestnodeid);
 
         int numslavenodes = -1;
         if (masternode != allcoupledcolnodes->end())

src/core/linear_solver/src/method/4C_linear_solver_method_parameters.cpp

@@ -86,9 +86,9 @@ void Core::LinearSolver::Parameters::compute_solver_parameters(
   {
     std::shared_ptr<Core::LinAlg::MultiVector<double>> coordinates;
     if (nullspace_node_map == nullptr)
-      coordinates = extract_node_coordinates(dis);
+      coordinates = extract_retained_node_coordinates(dis, *dis.node_row_map());
     else
-      coordinates = extract_node_coordinates(dis, *nullspace_node_map);
+      coordinates = extract_retained_node_coordinates(dis, *nullspace_node_map);
 
     solverlist.set<std::shared_ptr<Core::LinAlg::MultiVector<double>>>("Coordinates", coordinates);
   }