Merge pull request #3324 from roystgnr/distmesh_fixes

Poly2TriTriangulator fixes for DistributedMesh

Author: roystgnr

examples/subdomains/subdomains_ex3/subdomains_ex3.C

@@ -90,7 +90,7 @@ int main (int argc, char ** argv)
 
     if (argc == 3 && std::atoi(argv[2]) == 3)
       {
-        libmesh_here();
+        std::cout << "Running in 3D" << std::endl;
         dim=3;
       }
 

include/mesh/mesh_serializer.h

@@ -30,13 +30,26 @@ namespace libMesh
 class MeshBase;
 
 /**
- * Temporarily serialize a DistributedMesh for output; a distributed
- * mesh is allgathered by the MeshSerializer constructor if
- * need_serial is true, then remote elements are deleted again by the
- * destructor.
+ * Temporarily serialize a DistributedMesh for non-distributed-mesh
+ * capable code paths.  A distributed mesh is allgathered by the
+ * MeshSerializer constructor if need_serial is true, in which case
+ * remote elements are deleted again by the destructor.
+ *
+ * Serialization to processor 0 alone (e.g. for serial output from
+ * that processor) can be selected in the constructor.
+ *
+ * If allow_remote_element_removal() is set to true, that will also be
+ * temporarily disabled by the serializer, to be reenabled after
+ * serializer distruction if so; this allows prepare_for_use() to be
+ * called safely from within serialized code.
+ *
+ * If a mesh is explicitly distributed by a `delete_remote_elements()`
+ * call within serialized code, or if allow_remote_element_removal()
+ * is explicitly set to true within serialized code, the behavior is
+ * undefined.
  *
  * \author Roy Stogner
- * \date 2011
+ * \date 2011-2022
  * \brief Temporarily serializes a DistributedMesh for output.
  */
 class MeshSerializer
@@ -47,8 +60,24 @@ class MeshSerializer
   ~MeshSerializer();
 
 private:
+  /*
+   * The mesh which should remain serialized while this serializer
+   * exists
+   */
   MeshBase & _mesh;
+
+  /*
+   * Whether to delete remote elements on the mesh (returning it to a
+   * distributed state) when the serializer is destroyed
+   */
   bool reparallelize;
+
+  /*
+   * Whether to again allow `prepare_for_use` to delete remote
+   * elements on the mesh (returning it to a distributed state) after
+   * the serializer is destroyed
+   */
+  bool resume_allow_remote_element_removal;
 };
 
 } // namespace libMesh

include/mesh/poly2tri_triangulator.h

@@ -27,7 +27,6 @@
 
 // Local Includes
 #include "libmesh/dof_object.h"
-#include "libmesh/mesh_serializer.h"
 #include "libmesh/triangulator_interface.h"
 
 namespace libMesh
@@ -139,11 +138,6 @@ class Poly2TriTriangulator : public TriangulatorInterface
    */
   std::map<const Hole *, std::unique_ptr<ArbitraryHole>> replaced_holes;
 
-  /**
-   * We only operate on serialized meshes.
-   */
-  MeshSerializer _serializer;
-
   /**
    * Keep track of how many mesh nodes are boundary nodes.
    */

include/mesh/triangulator_interface.h

@@ -234,11 +234,17 @@ class TriangulatorInterface
 protected:
   /**
    * Helper function to create PSLG segments from our other
-   * boundary-defining options (nodal ordering, 1D mesh edges, 2D mesh
+   * boundary-defining options (1D mesh edges, 2D mesh
    * boundary sides), if no segments already exist.
    */
   void elems_to_segments();
 
+  /**
+   * Helper function to create PSLG segments from our node ordering,
+   * up to the maximum node id, if no segments already exist.
+   */
+  void nodes_to_segments(dof_id_type max_node_id);
+
   /**
    * Helper function to add extra points (midpoints of initial
    * segments) to a PSLG triangulation

src/error_estimation/exact_solution.C

@@ -24,6 +24,7 @@
 #include "libmesh/function_base.h"
 #include "libmesh/mesh_base.h"
 #include "libmesh/mesh_function.h"
+#include "libmesh/mesh_serializer.h"
 #include "libmesh/numeric_vector.h"
 #include "libmesh/parallel.h"
 #include "libmesh/quadrature.h"
@@ -475,9 +476,14 @@ void ExactSolution::_compute_error(std::string_view sys_name,
   const unsigned int var_component =
     computed_system.variable_scalar_number(var, 0);
 
-  // Prepare a global solution and a MeshFunction of the coarse system if we need one
+  // Prepare a global solution, a serialized mesh, and a MeshFunction
+  // of the coarse system, if we need them for fine-system integration
   std::unique_ptr<MeshFunction> coarse_values;
-  std::unique_ptr<NumericVector<Number>> comparison_soln = NumericVector<Number>::build(_equation_systems.comm());
+  std::unique_ptr<NumericVector<Number>> comparison_soln =
+    NumericVector<Number>::build(_equation_systems.comm());
+  MeshSerializer
+    serial(const_cast<MeshBase&>(_equation_systems.get_mesh()),
+           _equation_systems_fine);
   if (_equation_systems_fine)
     {
       const System & comparison_system

src/mesh/mesh_serializer.C

@@ -26,7 +26,8 @@ namespace libMesh
 
 MeshSerializer::MeshSerializer(MeshBase & mesh, bool need_serial, bool serial_only_needed_on_proc_0) :
   _mesh(mesh),
-  reparallelize(false)
+  reparallelize(false),
+  resume_allow_remote_element_removal(false)
 {
   libmesh_parallel_only(mesh.comm());
   if (need_serial && !_mesh.is_serial() && !serial_only_needed_on_proc_0) {
@@ -38,6 +39,12 @@ MeshSerializer::MeshSerializer(MeshBase & mesh, bool need_serial, bool serial_on
     // Just waste a bit of space on processor 0 to speed things up
     _mesh.gather_to_zero();
   }
+
+  if (_mesh.allow_remote_element_removal())
+    {
+      resume_allow_remote_element_removal = true;
+      _mesh.allow_remote_element_removal(false);
+    }
 }
 
 
@@ -46,6 +53,9 @@ MeshSerializer::~MeshSerializer()
 {
   if (reparallelize)
     _mesh.delete_remote_elements();
+
+  if (resume_allow_remote_element_removal)
+    _mesh.allow_remote_element_removal(true);
 }
 
 } // namespace libMesh

src/mesh/mesh_triangle_holes.C

@@ -25,6 +25,7 @@
 #include "libmesh/elem.h"
 #include "libmesh/int_range.h"
 #include "libmesh/mesh_base.h"
+#include "libmesh/mesh_serializer.h"
 #include "libmesh/node.h"
 #include "libmesh/parallel_algebra.h" // Packing<Point>
 #include "libmesh/simple_range.h"
@@ -262,6 +263,9 @@ TriangulatorInterface::MeshedHole::MeshedHole(const MeshBase & mesh,
   // pointers as keys.
   libmesh_parallel_only(mesh.comm());
 
+  MeshSerializer serial(const_cast<MeshBase &>(mesh),
+                        /* serial */ true, /* only proc 0 */ true);
+
   if (mesh.processor_id() != 0)
     {
       // Receive what proc 0 will send later

src/mesh/mesh_triangle_interface.C

@@ -60,12 +60,16 @@ void TriangleInterface::triangulate()
   // Will the triangulation have holes?
   const bool have_holes = ((_holes != nullptr) && (!_holes->empty()));
 
+  unsigned int n_hole_points = this->total_hole_points();
+
+  // If we're doing PSLG without segments, construct them from all our
+  // mesh nodes
+  this->nodes_to_segments(_mesh.max_node_id());
+
   // Insert additional new points in between existing boundary points,
   // if that is requested and reasonable
   this->insert_any_extra_boundary_points();
 
-  unsigned int n_hole_points = this->total_hole_points();
-
   // Regardless of whether we added additional points, the set of points to
   // triangulate is now sitting in the mesh.
 
@@ -339,6 +343,8 @@ void TriangleInterface::triangulate()
   }
 
 
+  _mesh.set_mesh_dimension(2);
+
 
   // To the naked eye, a few smoothing iterations usually looks better,
   // so we do this by default unless the user says not to.

src/mesh/poly2tri_triangulator.C

@@ -28,6 +28,7 @@
 #include "libmesh/enum_elem_type.h"
 #include "libmesh/function_base.h"
 #include "libmesh/hashing.h"
+#include "libmesh/mesh_serializer.h"
 #include "libmesh/mesh_smoother_laplace.h"
 #include "libmesh/mesh_triangle_holes.h"
 #include "libmesh/unstructured_mesh.h"
@@ -151,7 +152,6 @@ namespace libMesh
 Poly2TriTriangulator::Poly2TriTriangulator(UnstructuredMesh & mesh,
                                            dof_id_type n_boundary_nodes)
   : TriangulatorInterface(mesh),
-    _serializer(_mesh),
     _n_boundary_nodes(n_boundary_nodes),
     _refine_bdy_allowed(true)
 {
@@ -164,6 +164,11 @@ Poly2TriTriangulator::~Poly2TriTriangulator() = default;
 // Primary function responsible for performing the triangulation
 void Poly2TriTriangulator::triangulate()
 {
+  // We only operate on serialized meshes.  And it's not safe to
+  // serialize earlier, because it would then be possible for the user
+  // to re-parallelize the mesh in between there and here.
+  MeshSerializer serializer(_mesh);
+
   // We don't yet support every set of Triangulator options in the
   // poly2tri implementation
 
@@ -184,6 +189,14 @@ void Poly2TriTriangulator::triangulate()
       _elem_type != TRI7)
     libmesh_not_implemented();
 
+  // If we have no explicit segments defined, we may get them from
+  // mesh elements
+  this->elems_to_segments();
+
+  // If we *still* have no explicit segments defined, we get them from
+  // the order of nodes.
+  this->nodes_to_segments(_n_boundary_nodes);
+
   // Insert additional new points in between existing boundary points,
   // if that is requested and reasonable
   this->insert_any_extra_boundary_points();
@@ -199,6 +212,8 @@ void Poly2TriTriangulator::triangulate()
   if (_markers)
     libmesh_not_implemented();
 
+  _mesh.set_mesh_dimension(2);
+
   // To the naked eye, a few smoothing iterations usually looks better,
   // so we do this by default unless the user says not to.
   if (this->_smooth_after_generating)
@@ -273,9 +288,10 @@ void Poly2TriTriangulator::triangulate_current_points()
 
   // Unless we're using an explicit segments list, we assume node ids
   // are contiguous here.
+  dof_id_type nn = _mesh.max_node_id();
   if (this->segments.empty())
     libmesh_error_msg_if
-      (_mesh.n_nodes() != _mesh.max_node_id(),
+      (_mesh.n_nodes() != nn,
        "Poly2TriTriangulator needs contiguous node ids or explicit segments!");
 
   // And if we have more nodes than outer boundary points, the rest
@@ -301,10 +317,6 @@ void Poly2TriTriangulator::triangulate_current_points()
   // Prepare poly2tri points for our nodes, sorted into outer boundary
   // points and interior Steiner points.
 
-  // If we have no explicit segments defined, we may get them from
-  // mesh elements
-  this->elems_to_segments();
-
   if (this->segments.empty())
     {
       // If we have no segments even after taking elems into account,
@@ -455,7 +467,7 @@ void Poly2TriTriangulator::triangulate_current_points()
             }
           else
             {
-              Node * node = _mesh.add_point(p);
+              Node * node = _mesh.add_point(p, nn++);
               point_node_map[pt] = node;
             }
         }
@@ -577,6 +589,11 @@ bool Poly2TriTriangulator::insert_refinement_points()
 
   BoundaryInfo & boundary_info = _mesh.get_boundary_info();
 
+  // In cases where we've been working with contiguous node id ranges;
+  // let's keep it that way.
+  dof_id_type nn = _mesh.max_node_id();
+  dof_id_type ne = _mesh.max_elem_id();
+
   for (auto & elem : mesh.element_ptr_range())
     {
       // element_ptr_range skips deleted elements ... right?
@@ -663,13 +680,13 @@ bool Poly2TriTriangulator::insert_refinement_points()
                                                    side))
                 {
                   new_pt = ray_start;
-                  new_node = mesh.add_point(new_pt);
+                  new_node = mesh.add_point(new_pt, nn++);
                   boundary_refine(side);
                 }
               else
                 {
                   new_pt = cavity_elem->vertex_average();
-                  new_node = mesh.add_point(new_pt);
+                  new_node = mesh.add_point(new_pt, nn++);
                   // This was going to be a side refinement but it's
                   // now an internal refinement
                   side = invalid_uint;
@@ -724,21 +741,21 @@ bool Poly2TriTriangulator::insert_refinement_points()
                   // Let's just try bisecting for now
                   new_pt = (cavity_elem->point(side) +
                             cavity_elem->point((side+1)%3)) / 2;
-                  new_node = mesh.add_point(new_pt);
+                  new_node = mesh.add_point(new_pt, nn++);
                   boundary_refine(side);
                 }
               else // Do the best we can under these restrictions
                 {
                   new_pt = cavity_elem->vertex_average();
-                  new_node = mesh.add_point(new_pt);
+                  new_node = mesh.add_point(new_pt, nn++);
 
                   // This was going to be a side refinement but it's
                   // now an internal refinement
                   side = invalid_uint;
                 }
             }
           else
-            new_node = mesh.add_point(new_pt);
+            new_node = mesh.add_point(new_pt, nn++);
         }
       else
         libmesh_assert(new_node);
@@ -845,7 +862,7 @@ bool Poly2TriTriangulator::insert_refinement_points()
                       continue;
                     }
 
-                  auto new_elem = Elem::build(TRI3);
+                  auto new_elem = Elem::build_with_id(TRI3, ne++);
                   new_elem->set_node(0) = new_node;
                   new_elem->set_node(1) = node_CW;
                   new_elem->set_node(2) = node_CCW;