MeshTools::n_connected_components calculations

include/mesh/mesh_tools.h

@@ -280,6 +280,19 @@ unsigned int n_p_levels (const MeshBase & mesh);
  */
 unsigned int paranoid_n_levels(const MeshBase & mesh);
 
+/**
+ * \returns The number of topologically connected components in the
+ * mesh, where a local connection is defined as two nodes connected to
+ * the same element, two elements connected to the same node, and/or
+ * two nodes connected by a constraint row; thus a connection can be
+ * made at a single vertex, not just at an element side.
+ *
+ * To count sufficiently weak constraint row coefficients as a
+ * lack-of-connection, set a constraint_tol greater than 0.
+ */
+dof_id_type n_connected_components(const MeshBase & mesh,
+                                   Real constraint_tol = 0);
+
 /**
  * Builds a set of node IDs for nodes which belong to non-subactive
  * elements.  Non-subactive elements are those which are either active

src/apps/calculator.C

@@ -30,6 +30,7 @@
 #include "libmesh/exact_solution.h"
 #include "libmesh/getpot.h"
 #include "libmesh/mesh.h"
+#include "libmesh/mesh_tools.h"
 #include "libmesh/newton_solver.h"
 #include "libmesh/numeric_vector.h"
 #include "libmesh/parsed_fem_function.h"
@@ -58,7 +59,9 @@ void usage_error(const char * progname)
                << " --dim d               mesh dimension            [default: autodetect]\n"
                << " --inmesh    filename  input mesh file\n"
                << " --inmat     filename  input constraint matrix   [default: none]\n"
-               << " --insoln    filename  input solution file\n"
+               << " --mattol    filename  constraint tolerance when testing mesh connectivity\n"
+               << "                                                 [default: 0]\n"
+               << " --insoln    filename  input solution file\n     [default: none]\n"
                << " --calc      func      function to calculate\n"
                << " --insys     sysnum    input system number       [default: 0]\n"
                << " --outsoln   filename  output solution file      [default: out_<insoln>]\n"
@@ -192,6 +195,16 @@ int main(int argc, char ** argv)
   libMesh::out << "Mesh:" << std::endl;
   old_mesh.print_info();
 
+  // If we're not using a distributed mesh, this is cheap info to add
+  if (old_mesh.is_serial_on_zero())
+    {
+      const Real mat_tol = cl.follow(Real(0), "--mattol");
+
+      const dof_id_type n_components =
+        MeshTools::n_connected_components(old_mesh, mat_tol);
+      libMesh::out << "Mesh has " << n_components << " connected components." << std::endl;
+    }
+
   const std::string solnname = cl.follow(std::string(""), "--insoln");
 
   // Load the old solution from --insoln filename, if that's been

src/error_estimation/exact_solution.C

@@ -659,6 +659,11 @@ void ExactSolution::_compute_error(std::string_view sys_name,
       // Build a quadrature rule.
       q_rules[dim] = fe_type.default_quadrature_rule (dim, _extra_order);
 
+      // Disallow rules with negative weights.  That will use more
+      // quadrature points, but we're going to be taking square roots
+      // of element integral results here!
+      q_rules[dim]->allow_rules_with_negative_weights = false;
+
       // Construct finite element object
       fe_ptrs[dim] = FEGenericBase<OutputShape>::build(dim, fe_type);
 

src/mesh/mesh_base.C

@@ -2083,8 +2083,20 @@ MeshBase::copy_constraint_rows(const SparseMatrix<T> & constraint_operator)
       if (indices.size() == 1 &&
           values[0] == T(1))
         {
-          existing_unconstrained_nodes.insert(i);
-          existing_unconstrained_columns.emplace(indices[0],i);
+          // If we have multiple simple Ui=Uj constraints, let the
+          // first one be our "unconstrained" node and let the others
+          // be constrained to it.
+          if (existing_unconstrained_columns.find(indices[0]) !=
+              existing_unconstrained_columns.end())
+            {
+              const auto j = indices[0];
+              columns[j].emplace_back(i, 1);
+            }
+          else
+            {
+              existing_unconstrained_nodes.insert(i);
+              existing_unconstrained_columns.emplace(indices[0],i);
+            }
         }
       else
         for (auto jj : index_range(indices))

src/mesh/mesh_tools.C

@@ -23,6 +23,7 @@
 #include "libmesh/libmesh_logging.h"
 #include "libmesh/mesh_base.h"
 #include "libmesh/mesh_communication.h"
+#include "libmesh/mesh_serializer.h"
 #include "libmesh/mesh_tools.h"
 #include "libmesh/node_range.h"
 #include "libmesh/parallel.h"
@@ -836,6 +837,118 @@ unsigned int paranoid_n_levels(const MeshBase & mesh)
 
 
 
+dof_id_type n_connected_components(const MeshBase & mesh,
+                                   Real constraint_tol)
+{
+  // Yes, I'm being lazy.  This is for mesh analysis before a
+  // simulation, not anything going in any loops.
+  if (!mesh.is_serial_on_zero())
+    libmesh_not_implemented();
+
+  dof_id_type n_components = 0;
+
+  if (mesh.processor_id())
+  {
+    mesh.comm().broadcast(n_components);
+    return n_components;
+  }
+
+  // All nodes in a set here are connected (at least indirectly) to
+  // all other nodes in the same set, but have not yet been discovered
+  // to be connected to nodes in other sets.
+  std::vector<std::unordered_set<const Node *>> components;
+
+  auto find_component = [&components](const Node * n) {
+    std::unordered_set<const Node *> * component = nullptr;
+
+    for (auto & c: components)
+      if (c.find(n) != c.end())
+        {
+          libmesh_assert(component == nullptr);
+          component = &c;
+        }
+
+    return component;
+  };
+
+  auto add_to_component =
+    [&find_component]
+    (std::unordered_set<const Node *> & component, const Node * n) {
+
+    auto current_component = find_component(n);
+    // We may already know we're in the desired component
+    if (&component == current_component)
+      return;
+
+    // If we're unknown, we should be in the desired component
+    else if (!current_component)
+      component.insert(n);
+
+    // If we think we're in another component, it should actually be
+    // part of the desired component
+    else
+      {
+        component.merge(*current_component);
+        libmesh_assert(current_component->empty());
+      }
+  };
+
+  auto & constraint_rows = mesh.get_constraint_rows();
+
+  for (const auto & elem : mesh.element_ptr_range())
+    {
+      const Node * first_node = elem->node_ptr(0);
+
+      auto component = find_component(first_node);
+
+      // If we didn't find one, make a new one, reusing an existing
+      // slot if possible or growing our vector if necessary
+      if (!component)
+        for (auto & c: components)
+          if (c.empty())
+            component = &c;
+
+      if (!component)
+        component = &components.emplace_back();
+
+      for (const Node & n : elem->node_ref_range())
+        {
+          add_to_component(*component, &n);
+
+          auto it = constraint_rows.find(&n);
+          if (it == constraint_rows.end())
+            continue;
+
+          for (const auto & [pr, val] : it->second)
+            {
+              // Ignore too-trivial constraint coefficients if
+              // we get a non-default-0 constraint_tol
+              if (std::abs(val) < constraint_tol)
+                continue;
+
+              const Elem * spline_elem = pr.first;
+              libmesh_assert(spline_elem == mesh.elem_ptr(spline_elem->id()));
+
+              const Node * spline_node =
+                spline_elem->node_ptr(pr.second);
+
+              add_to_component(*component, spline_node);
+            }
+        }
+    }
+
+  for (auto & component : components)
+    if (!component.empty())
+      ++n_components;
+
+  // We calculated this on proc 0; now let everyone else know too
+  mesh.comm().broadcast(n_components);
+
+  return n_components;
+}
+
+
+
 void get_not_subactive_node_ids(const MeshBase & mesh,
                                 std::set<dof_id_type> & not_subactive_node_ids)
 {

tests/Makefile.am

@@ -61,6 +61,7 @@ unit_tests_sources = \
   mesh/boundary_info.C \
   mesh/boundary_points.C \
   mesh/checkpoint.C \
+  mesh/connected_components.C \
   mesh/contains_point.C \
   mesh/distributed_mesh_test.C \
   mesh/extra_integers.C \