Distributed mesh fixes

include/geom/elem.h

@@ -290,13 +290,19 @@ class Elem : public ReferenceCountedObject<Elem>,
   virtual dof_id_type key () const;
 
   /**
-   * \returns \p true if two elements are equivalent, false otherwise.
-   * This is true if the elements are connected to identical global
-   * nodes, regardless of how those nodes might be numbered local
-   * to the elements.
+   * \returns \p true if two elements are equivalent, \p false
+   * otherwise.  This is true if the elements are connected to
+   * identical global nodes, regardless of how those nodes might be
+   * numbered local to the elements.
    */
   bool operator == (const Elem & rhs) const;
 
+  /**
+   * \returns \p false if two elements are equivalent, \p true
+   * otherwise.
+   */
+  bool operator != (const Elem & rhs) const;
+
   /**
    * \returns \p true if two elements have equal topologies, false
    * otherwise.
@@ -1585,6 +1591,21 @@ class Elem : public ReferenceCountedObject<Elem>,
                             unsigned int side,
                             bool reset = true);
 
+  /**
+   * Same as the \p total_family_tree() member, but only adds elements
+   * which are next to \p side.
+   */
+  void total_family_tree_by_side (std::vector<const Elem *> & family,
+                                  unsigned int side,
+                                  bool reset = true) const;
+
+  /**
+   * Non-const version of function above; fills a vector of non-const pointers.
+   */
+  void total_family_tree_by_side (std::vector<Elem *> & family,
+                                  unsigned int side,
+                                  bool reset = true);
+
   /**
    * Same as the \p active_family_tree() member, but only adds elements
    * which are next to \p side.
@@ -2579,6 +2600,14 @@ subdomain_id_type & Elem::subdomain_id ()
 
 
 
+inline
+bool Elem::operator != (const Elem & rhs) const
+{
+  return !(*this == rhs);
+}
+
+
+
 inline
 const Elem * Elem::neighbor_ptr (unsigned int i) const
 {

include/geom/elem_internal.h

@@ -144,6 +144,35 @@ family_tree_by_side (T elem,
 
 
 
+template <class T>
+void
+total_family_tree_by_side (T elem,
+                           std::vector<T> & family,
+                           unsigned int s,
+                           bool reset)
+{
+  // Clear the vector if the flag reset tells us to.
+  if (reset)
+    family.clear();
+
+  libmesh_assert_less (s, elem->n_sides());
+
+  // Add this element to the family tree.
+  family.push_back(elem);
+
+  // Recurse into the elements children, if it has them.
+  // Do not clear the vector any more.
+  if (elem->has_children())
+    {
+      const unsigned int nc = elem->n_children();
+      for (unsigned int c = 0; c != nc; c++)
+        if (!elem->child_ptr(c)->is_remote() && elem->is_child_on_side(c, s))
+          ElemInternal::total_family_tree_by_side (elem->child_ptr(c), family, s, false);
+    }
+}
+
+
+
 template <class T>
 void
 active_family_tree_by_side (T elem,

src/geom/elem.C

@@ -1430,25 +1430,22 @@ void Elem::make_links_to_me_local(unsigned int n, unsigned int nn)
   libmesh_assert(neigh);
   libmesh_assert(!neigh->is_remote());
 
+  const unsigned int this_level = this->level();
+  const unsigned int neigh_level = neigh->level();
+
   // We never have neighbors more refined than us
-  libmesh_assert_less_equal (neigh->level(), this->level());
+  libmesh_assert_less_equal (neigh_level, this_level);
 
   // We never have subactive neighbors of non subactive elements
   libmesh_assert(!neigh->subactive() || this->subactive());
 
   // If we have a neighbor less refined than us then it must not
   // have any more refined descendants we could have pointed to
   // instead.
-  libmesh_assert((neigh->level() == this->level()) ||
+  libmesh_assert((neigh_level == this_level) ||
                  (neigh->active() && !this->subactive()) ||
                  (!neigh->has_children() && this->subactive()));
 
-  // If neigh is at our level, then its family might have
-  // remote_elem neighbor links which need to point to us
-  // instead, but if not, then we're done.
-  if (neigh->level() != this->level())
-    return;
-
   // What side of neigh are we on?  nn.
   //
   // We can't use the usual Elem method because we're in the middle of
@@ -1459,20 +1456,27 @@ void Elem::make_links_to_me_local(unsigned int n, unsigned int nn)
   // not-technically-neighbors until they're
   // not-even-geometrically-neighbors.
 
-  // Find any elements that ought to point to elem
+  // Find any elements that might need to point to elem
   std::vector<Elem *> neigh_family;
-#ifdef LIBMESH_ENABLE_AMR
-  if (this->active())
-    neigh->family_tree_by_side(neigh_family, nn);
-  else
-#endif
-    neigh_family.push_back(neigh);
+  neigh->total_family_tree_by_side(neigh_family, nn);
+
+  // Pull objects out of the loop to reduce heap operations
+  std::unique_ptr<Elem> my_side, neigh_side;
 
   // And point them to elem
   for (auto & neigh_family_member : neigh_family)
     {
       // Only subactive elements point to other subactive elements
-      if (this->subactive() && !neigh_family_member->subactive())
+      const bool member_subactive = neigh_family_member->subactive();
+      if (this->subactive() && !member_subactive)
+        continue;
+
+      // neigh (and possibly some of its family) might be at a lower
+      // level than us, with a neighbor at that lower level.  We
+      // would exit early if neigh was at a lower level, except we do
+      // want to handle neighbor links from subactive descendents.
+      const unsigned int member_level = neigh_family_member->level();
+      if (member_level < this_level)
         continue;
 
       // Ideally, the neighbor link ought to either be correct
@@ -1485,6 +1489,7 @@ void Elem::make_links_to_me_local(unsigned int n, unsigned int nn)
 #ifdef LIBMESH_ENABLE_AMR
       libmesh_assert((neigh_family_member->neighbor_ptr(nn) &&
                       (neigh_family_member->neighbor_ptr(nn)->active() ||
+                       this->is_ancestor_of(neigh_family_member->neighbor_ptr(nn)) ||
                        neigh_family_member->neighbor_ptr(nn)->is_ancestor_of(this))) ||
                      (neigh_family_member->neighbor_ptr(nn) == remote_elem) ||
                      ((this->refinement_flag() == JUST_REFINED) &&
@@ -1495,6 +1500,32 @@ void Elem::make_links_to_me_local(unsigned int n, unsigned int nn)
                      (neigh_family_member->neighbor_ptr(nn) == remote_elem));
 #endif
 
+      // Some of neigh's family might be at a higher (finer) level
+      // than us, and might need to point back to one of our children
+      // rather than to us.
+      if (member_level > this_level)
+        {
+          if (this->ancestor())
+            continue;
+
+          if (member_subactive && this->has_children())
+            continue;
+        }
+
+      // If neigh is at a coarser level than us, some of neigh's
+      // subactive family just won't line up with us!
+      if (neigh_level < this_level &&
+          member_level > neigh_level)
+        {
+          libmesh_assert(member_subactive);
+
+          this->side_ptr(my_side, n);
+          neigh_family_member->side_ptr(neigh_side, nn);
+
+          if (*my_side != *neigh_side)
+            continue;
+        }
+
       neigh_family_member->set_neighbor(nn, this);
     }
 }
@@ -2134,6 +2165,24 @@ void Elem:: family_tree_by_side (std::vector<Elem *> & family,
 
 
 
+void Elem::total_family_tree_by_side (std::vector<const Elem *> & family,
+                                      unsigned int side,
+                                      bool reset) const
+{
+  ElemInternal::total_family_tree_by_side(this, family, side, reset);
+}
+
+
+
+void Elem::total_family_tree_by_side (std::vector<Elem *> & family,
+                                      unsigned int side,
+                                      bool reset)
+{
+  ElemInternal::total_family_tree_by_side(this, family, side, reset);
+}
+
+
+
 void Elem::active_family_tree_by_side (std::vector<const Elem *> & family,
                                        unsigned int side,
                                        bool reset) const

src/geom/elem_refinement.C

@@ -140,6 +140,23 @@ void Elem::refine (MeshRefinement & mesh_refinement)
         }
     }
 
+  // Get any new remote neighbor links right; even find_neighbors
+  // relies on those
+  for (unsigned int s : this->side_index_range())
+    {
+      if (this->neighbor_ptr(s) != remote_elem)
+        continue;
+
+      for (unsigned int c = 0; c != nc; c++)
+        {
+          Elem * current_child = this->child_ptr(c);
+          if (current_child != remote_elem &&
+              this->is_child_on_side(c, s))
+            current_child->set_neighbor
+              (s, const_cast<RemoteElem *>(remote_elem));
+        }
+    }
+
   // Un-set my refinement flag now
   this->set_refinement_flag(Elem::INACTIVE);