Merge pull request #4394 from roystgnr/distributed_inffe_fixes

Distributed infinite element mesh building fixes

Author: roystgnr

examples/miscellaneous/miscellaneous_ex15/miscellaneous_ex15.C

@@ -50,6 +50,8 @@
 #include "libmesh/error_vector.h"
 #include "libmesh/mesh_refinement.h"
 #include "libmesh/kelly_error_estimator.h"
+// for distributed mesh editing
+#include "libmesh/parallel_ghost_sync.h"
 
 // Bring in everything from the libMesh namespace
 using namespace libMesh;
@@ -65,7 +67,7 @@ void assemble_func(EquationSystems & es, const std::string & system_name)
   std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
   std::unique_ptr<FEBase> inf_fe (FEBase::build_InfFE(dim, fe_type));
 
-  const std::vector<dof_id_type > charged_objects = es.parameters.get<std::vector<dof_id_type> >("charged_elem_id");
+  const std::set<dof_id_type > charged_objects = es.parameters.get<std::set<dof_id_type> >("charged_elem_id");
   const auto qrule=fe_type.default_quadrature_rule(/*dim = */ 3, /*extra order = */ 3);
 
   fe->attach_quadrature_rule (&*qrule);
@@ -115,7 +117,7 @@ void assemble_func(EquationSystems & es, const std::string & system_name)
 
       Real rho;
       // check if charged_objects contains elem->id(), we add the corresponding charge to the rhs vector
-      if (std::find(charged_objects.begin(), charged_objects.end(), elem->id()) != charged_objects.end())
+      if (charged_objects.count(elem->id()))
         rho=1./elem->volume();
       else
         rho = 0;
@@ -178,6 +180,10 @@ int main (int argc, char** argv)
   // creation of an empty mesh-object
   Mesh mesh(init.comm(), dim);
 
+  // We're going to keep track of elements by id, so we don't want
+  // even a DistributedMesh to renumber them.
+  mesh.allow_renumbering(false);
+
   // fill the meshes with a spherical grid of type HEX27 with radius r
   MeshTools::Generation::build_cube (mesh, /*nx= */5, /*ny=*/5, /*nz=*/3,
                                      /*xmin=*/ -5.2, /*xmax=*/5.2,
@@ -205,31 +211,44 @@ int main (int argc, char** argv)
         elem->neighbor_ptr(0)->subdomain_id()=2;
       }
 
+  // If we're on a distributed mesh then we might have missed some
+  // ghosted base elements with remote infinite elements.
+  if (!mesh.is_serial())
+    {
+      SyncSubdomainIds sync_obj(mesh);
+      Parallel::sync_dofobject_data_by_id
+        (mesh.comm(), mesh.elements_begin(), mesh.elements_end(),
+         sync_obj);
+    }
 
   // Now we set the sources of the field: prism-shaped objects that are
   // determined here by containing certain points:
-  PointLocatorTree pt_lctr(mesh);
-  std::vector<dof_id_type> charged_elem_ids(3);
+  auto p_pt_lctr = mesh.sub_point_locator();
+  auto & pt_lctr = *p_pt_lctr;
+  pt_lctr.enable_out_of_mesh_mode();
+
+  std::set<dof_id_type> charged_elem_ids;
   {
     Point pt_0(-3.,-3.0,-1.5);
     Point pt_1(2.,-2.6,-1.5);
     Point pt_2(2., 3.1, 1.7);
     const Elem * elem_0=pt_lctr(pt_0);
     if (elem_0)
-      charged_elem_ids[0]=elem_0->id();
-    else
-      // this indicates some error which I don't know how to handle.
-      libmesh_not_implemented();
+      charged_elem_ids.insert(elem_0->id());
     const Elem * elem_1=pt_lctr(pt_1);
     if (elem_1)
-      charged_elem_ids[1]=elem_1->id();
-    else
-      libmesh_not_implemented();
+      charged_elem_ids.insert(elem_1->id());
     const Elem * elem_2=pt_lctr(pt_2);
     if (elem_2)
-      charged_elem_ids[2]=elem_2->id();
-    else
-      libmesh_not_implemented();
+      charged_elem_ids.insert(elem_2->id());
+
+    // On a distributed mesh we might not have every point in a
+    // semilocal element on every processor
+    if (!mesh.is_serial())
+      mesh.comm().set_union(charged_elem_ids);
+
+    // But we should have every point on *some* processor
+    libmesh_assert_equal_to(charged_elem_ids.size(), 3);
   }
 
   // Create an equation systems object
@@ -238,7 +257,7 @@ int main (int argc, char** argv)
   // This is the only system added here.
   LinearImplicitSystem & eig_sys = eq_sys.add_system<LinearImplicitSystem> ("Poisson");
 
-  eq_sys.parameters.set<std::vector<dof_id_type> >("charged_elem_id")=charged_elem_ids;
+  eq_sys.parameters.set<std::set<dof_id_type> >("charged_elem_id")=charged_elem_ids;
 
   //set the complete type of the variable
   FEType fe_type(SECOND, LAGRANGE, FOURTH, JACOBI_20_00, CARTESIAN);
@@ -272,21 +291,31 @@ int main (int argc, char** argv)
       // in the refined mesh, find the elements that describe the
       // sources of gravitational field: Due to refinement, there are
       // successively more elements to account for the same object.
-      std::vector<dof_id_type> charged_elem_child(0);
-      for(unsigned int j=0; j< charged_elem_ids.size(); ++j)
+      std::set<dof_id_type> charged_elem_children;
+      for(auto id : charged_elem_ids)
         {
-          Elem * charged_elem = mesh.elem_ptr(charged_elem_ids[j]);
+          Elem * charged_elem = mesh.query_elem_ptr(id);
+          if (!charged_elem)
+            {
+              libmesh_assert(!mesh.is_serial());
+              continue;
+            }
+
           if (charged_elem->has_children())
             {
               for(auto & child : charged_elem->child_ref_range())
-                charged_elem_child.push_back(child.id());
+                if (!child.is_remote())
+                  charged_elem_children.insert(child.id());
             }
           else
-            charged_elem_child.push_back(charged_elem->id());
+            charged_elem_children.insert(charged_elem->id());
         }
 
-      charged_elem_ids=charged_elem_child;
-      eq_sys.parameters.set<std::vector<dof_id_type> >("charged_elem_id")=charged_elem_ids;
+      charged_elem_ids=charged_elem_children;
+      if (!mesh.is_serial())
+        mesh.comm().set_union(charged_elem_ids);
+
+      eq_sys.parameters.set<std::set<dof_id_type> >("charged_elem_id")=charged_elem_ids;
 
       // re-assemble and than solve again.
       eig_sys.solve();

include/mesh/mesh_base.h

@@ -973,14 +973,20 @@ class MeshBase : public ParallelObject
    * After this routine is called all the elements with a \p nullptr neighbor
    * pointer are guaranteed to be on the boundary.  Thus this routine is
    * useful for automatically determining the boundaries of the domain.
-   * If reset_remote_elements is left to false, remote neighbor links are not
-   * reset and searched for in the local mesh.  If reset_current_list is
-   * left as true, then any existing links will be reset before initiating
-   * the algorithm, while honoring the value of the reset_remote_elements
-   * flag.
+   *
+   * If \p reset_remote_elements is left to false, remote neighbor
+   * links are not reset and searched for in the local mesh.
+   *
+   * If \p reset_current_list is left as true, then any existing links
+   * will be reset before initiating the algorithm, while honoring the
+   * value of the \p reset_remote_elements flag.
+   *
+   * If \p assert_valid is left as true, then in dbg mode extensive
+   * consistency checking is performed before returning.
    */
   virtual void find_neighbors (const bool reset_remote_elements = false,
-                               const bool reset_current_list    = true) = 0;
+                               const bool reset_current_list    = true,
+                               const bool assert_valid          = true) = 0;
 
   /**
    * Removes any orphaned nodes, nodes not connected to any elements.

include/mesh/unstructured_mesh.h

@@ -293,7 +293,8 @@ class UnstructuredMesh : public MeshBase
    * Other functions from MeshBase requiring re-definition.
    */
   virtual void find_neighbors (const bool reset_remote_elements = false,
-                               const bool reset_current_list    = true) override;
+                               const bool reset_current_list    = true,
+                               const bool assert_valid          = true) override;
 
 #ifdef LIBMESH_ENABLE_AMR
   /**

src/mesh/inf_elem_builder.C

@@ -35,6 +35,8 @@
 #include "libmesh/mesh_base.h"
 #include "libmesh/remote_elem.h"
 
+#include "timpi/parallel_implementation.h"
+
 namespace libMesh
 {
 
@@ -387,10 +389,14 @@ void InfElemBuilder::build_inf_elem(const Point & origin,
 
         } // neighbor(s) == nullptr
 
-
-
-
-
+  // On a distributed mesh it might be some other processor who sees
+  // the farthest node.
+  if (!this->_mesh.is_serial())
+    {
+      unsigned int rank;
+      this->_mesh.comm().maxloc(max_r, rank);
+      this->_mesh.comm().broadcast(max_r_node, rank);
+    }
 
   //  If a boundary side has one node on the outer boundary,
   //  all points of this side are on the outer boundary.
@@ -405,8 +411,14 @@ void InfElemBuilder::build_inf_elem(const Point & origin,
     // Only the case of no outer boundary is to be excluded.
     onodes.insert(max_r_node);
 
+  // If we're not on a serial mesh, we'll need to synchronize that
+  // onodes list too.
+  bool did_parallel_update;
 
+  do
   {
+    did_parallel_update = false;
+
     auto face_it = faces.begin();
     auto face_end = faces.end();
     unsigned int facesfound=0;
@@ -448,7 +460,16 @@ void InfElemBuilder::build_inf_elem(const Point & origin,
           face_it    = faces.begin();
         }
     }
+
+    if (!this->_mesh.is_serial())
+      {
+        auto my_onodes_size = onodes.size();
+        this->_mesh.comm().set_union(onodes);
+        did_parallel_update = (onodes.size() > my_onodes_size);
+        this->_mesh.comm().max(did_parallel_update);
+      }
   }
+  while (did_parallel_update);
 
 
 #ifdef DEBUG
@@ -486,26 +507,32 @@ void InfElemBuilder::build_inf_elem(const Point & origin,
 
   // for each boundary node, add an outer_node with
   // double distance from origin.
-  for (const auto & dof : onodes)
+  for (const auto & n : onodes)
     {
-      Point p = (Point(this->_mesh.point(dof)) * 2) - origin;
+      if (!this->_mesh.query_node_ptr(n))
+        {
+          libmesh_assert(!_mesh.is_serial());
+          continue;
+        }
+
+      Point p = (Point(this->_mesh.point(n)) * 2) - origin;
       if (_mesh.is_serial())
         {
           // Add with a default id in serial
-          outer_nodes[dof]=this->_mesh.add_point(p);
+          outer_nodes[n]=this->_mesh.add_point(p);
         }
       else
         {
           // Pick a unique id in parallel
-          Node & bnode = _mesh.node_ref(dof);
+          Node & bnode = _mesh.node_ref(n);
           dof_id_type new_id = bnode.id() + old_max_node_id;
           std::unique_ptr<Node> new_node = Node::build(p, new_id);
           new_node->processor_id() = bnode.processor_id();
 #ifdef LIBMESH_ENABLE_UNIQUE_ID
           new_node->set_unique_id(old_max_unique_id + bnode.id());
 #endif
 
-          outer_nodes[dof] =
+          outer_nodes[n] =
             this->_mesh.add_node(std::move(new_node));
         }
     }
@@ -665,6 +692,32 @@ void InfElemBuilder::build_inf_elem(const Point & origin,
       this->_mesh.add_elem(std::move(el));
     } // for
 
+  // If we're not in serial, we might not have all our remote_elem
+  // neighbors set up properly on the new elements: two new infinite
+  // elements which should be neighbors may not be rooted in finite
+  // elmeents which are neighbors, and we were only looking for
+  // remote_elem links on each root finite element.  We'll fix the
+  // problem by first finding all the neighbors we can and then
+  // recognizing that any missing neighbors on infinite elements must
+  // be remote.
+  if (!_mesh.is_serial())
+    {
+      _mesh.find_neighbors(/*reset_remote_elements=*/ false,
+                           /*reset_current_list=*/    true,
+                           /*assert_valid=*/          false);
+
+      for (auto & p : ofaces)
+        {
+          Elem & belem = this->_mesh.elem_ref(p.first);
+          Elem * inf_elem = belem.neighbor_ptr(p.second);
+          libmesh_assert(inf_elem);
+
+          for (auto s : make_range(inf_elem->n_sides()))
+            if (!inf_elem->neighbor_ptr(s))
+              inf_elem->set_neighbor
+                (s, const_cast<RemoteElem *>(remote_elem));
+        }
+    }
 
 #ifdef DEBUG
   _mesh.libmesh_assert_valid_parallel_ids();

src/mesh/unstructured_mesh.C

@@ -941,7 +941,8 @@ UnstructuredMesh::~UnstructuredMesh ()
 
 
 void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
-                                       const bool reset_current_list)
+                                       const bool reset_current_list,
+                                       const bool assert_valid)
 {
   // We might actually want to run this on an empty mesh
   // (e.g. the boundary mesh for a nonexistent bcid!)
@@ -1326,9 +1327,14 @@ void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
 #endif // AMR
 
 #ifdef DEBUG
-  MeshTools::libmesh_assert_valid_neighbors(*this,
-                                            !reset_remote_elements);
-  MeshTools::libmesh_assert_valid_amr_interior_parents(*this);
+  if (assert_valid)
+    {
+      MeshTools::libmesh_assert_valid_neighbors(*this,
+                                                !reset_remote_elements);
+      MeshTools::libmesh_assert_valid_amr_interior_parents(*this);
+    }
+#else
+  libmesh_ignore(assert_valid);
 #endif
 
   this->_preparation.has_neighbor_ptrs = true;