From 1c40faf25d0ef3b6ea9a7c901c6b657366777e5b Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 24 Jun 2022 16:55:07 -0500 Subject: [PATCH 01/13] Explicit node ids for triangulator interpolation This fixes most of the bugs I see when running with DistributedMesh --- src/mesh/triangulator_interface.C | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/mesh/triangulator_interface.C b/src/mesh/triangulator_interface.C index 7ae5a288270..e58078be04d 100644 --- a/src/mesh/triangulator_interface.C +++ b/src/mesh/triangulator_interface.C @@ -155,7 +155,8 @@ void TriangulatorInterface::insert_any_extra_boundary_points() // Insert a new point on each PSLG at evenly spaced locations // between existing boundary points. - // np=index into new points vector + // np=index into new points vector, also an explicit node index + // (to avoid the stride in automatic DistributedMesh indexing) // n =index into original points vector for (std::size_t np=0, n=0, tops=n_interpolated*original_points.size(); np Date: Fri, 24 Jun 2022 17:31:56 -0500 Subject: [PATCH 02/13] Serialize for poly2tri right when we need it This fixes another few --enable-distmesh unit test bugs for me. --- include/mesh/poly2tri_triangulator.h | 6 ------ src/mesh/poly2tri_triangulator.C | 7 ++++++- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/include/mesh/poly2tri_triangulator.h b/include/mesh/poly2tri_triangulator.h index 6250974202e..be55467d882 100644 --- a/include/mesh/poly2tri_triangulator.h +++ b/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> replaced_holes; - /** - * We only operate on serialized meshes. - */ - MeshSerializer _serializer; - /** * Keep track of how many mesh nodes are boundary nodes. */ diff --git a/src/mesh/poly2tri_triangulator.C b/src/mesh/poly2tri_triangulator.C index d5ecd79d0d3..6439a56d687 100644 --- a/src/mesh/poly2tri_triangulator.C +++ b/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 From e809d62049529816d46aec51e68e0aee37267de7 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Sat, 25 Jun 2022 07:16:29 -0500 Subject: [PATCH 03/13] Don't let DistributedMesh number poly2tri nodes I don't this actually fixed any bugs, in hindsight, but it doesn't hurt to use an a priori numbering. --- src/mesh/poly2tri_triangulator.C | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/src/mesh/poly2tri_triangulator.C b/src/mesh/poly2tri_triangulator.C index 6439a56d687..76b95dda18a 100644 --- a/src/mesh/poly2tri_triangulator.C +++ b/src/mesh/poly2tri_triangulator.C @@ -278,9 +278,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 @@ -460,7 +461,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; } } @@ -582,6 +583,10 @@ 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(); + for (auto & elem : mesh.element_ptr_range()) { // element_ptr_range skips deleted elements ... right? @@ -668,13 +673,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; @@ -729,13 +734,13 @@ 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 @@ -743,7 +748,7 @@ bool Poly2TriTriangulator::insert_refinement_points() } } else - new_node = mesh.add_point(new_pt); + new_node = mesh.add_point(new_pt, nn++); } else libmesh_assert(new_node); From 0d639e84a959c628a4ff728573964af48608efe2 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Sat, 25 Jun 2022 07:17:45 -0500 Subject: [PATCH 04/13] Stop automatic reparallelization inside serializer This lets us call `prepare_for_use()` within serialized sections, as I'd been trying to do in our poly2tri code. This fixes all our unit tests with DistributedMesh in parallel for me. --- include/mesh/mesh_serializer.h | 39 +++++++++++++++++++++++++++++----- src/mesh/mesh_serializer.C | 12 ++++++++++- 2 files changed, 45 insertions(+), 6 deletions(-) diff --git a/include/mesh/mesh_serializer.h b/include/mesh/mesh_serializer.h index a5860380a1f..438ada7ba65 100644 --- a/include/mesh/mesh_serializer.h +++ b/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 diff --git a/src/mesh/mesh_serializer.C b/src/mesh/mesh_serializer.C index 3b81ef5addb..9cb1f8a6cab 100644 --- a/src/mesh/mesh_serializer.C +++ b/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 From bfab3dc93a9f669f790c2ad56d4fbf2007e3fc76 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 28 Jun 2022 16:51:46 -0500 Subject: [PATCH 05/13] Add missing header --- tests/mesh/write_elemset_data.C | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/mesh/write_elemset_data.C b/tests/mesh/write_elemset_data.C index 97c3dfa08ca..b79320f8483 100644 --- a/tests/mesh/write_elemset_data.C +++ b/tests/mesh/write_elemset_data.C @@ -5,6 +5,7 @@ #include "libmesh/mesh.h" #include "libmesh/mesh_generation.h" #include "libmesh/parallel.h" // set_union +#include "libmesh/replicated_mesh.h" #include "libmesh/string_to_enum.h" #include "libmesh/boundary_info.h" #include "libmesh/utility.h" // libmesh_map_find From afc3769bebf6b7247d8ef52e47631fc1b7f25430 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 29 Jun 2022 15:46:42 -0500 Subject: [PATCH 06/13] Use MeshSerializer in MeshedHole --- src/mesh/mesh_triangle_holes.C | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/mesh/mesh_triangle_holes.C b/src/mesh/mesh_triangle_holes.C index 03f5ff2cd8d..4818acc70a3 100644 --- a/src/mesh/mesh_triangle_holes.C +++ b/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 #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(mesh), + /* serial */ true, /* only proc 0 */ true); + if (mesh.processor_id() != 0) { // Receive what proc 0 will send later From 0a2a37b506afc3445e48450974dfa47d99aa1d20 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 29 Jun 2022 15:47:28 -0500 Subject: [PATCH 07/13] Try testBadElemFECombo only if min(n_local_elem) In a real application we're currently happy just to get a sane error message from one rank before we die, but here we really need to be able to recover afterward, not lose sync. --- tests/base/dof_map_test.C | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/tests/base/dof_map_test.C b/tests/base/dof_map_test.C index 8f096b36ece..328d88b8649 100644 --- a/tests/base/dof_map_test.C +++ b/tests/base/dof_map_test.C @@ -4,6 +4,8 @@ #include #include +#include + #include "test_comm.h" #include "libmesh_cppunit.h" @@ -123,6 +125,7 @@ public: void testBadElemFECombo() { LOG_UNIT_TEST; + Mesh mesh(*TestCommWorld); EquationSystems es(mesh); @@ -131,6 +134,15 @@ public: MeshTools::Generation::build_square (mesh,4,4,-1., 1.,-1., 1., QUAD4); + // We need at least one element per processor to make sure + // everyone throws and we don't get out of sync before going on to + // future tests. + dof_id_type min_local_elem = mesh.n_local_elem(); + mesh.comm().min(min_local_elem); + + if (!min_local_elem) + return; + // We can't just CPPUNIT_ASSERT_THROW, because we want to make // sure we were thrown from the right place with the right error // message! @@ -144,12 +156,14 @@ public: threw_desired_exception = true; } catch (...) { - CPPUNIT_ASSERT(false); + CPPUNIT_ASSERT_MESSAGE("Unexpected exception type thrown", false); } - CPPUNIT_ASSERT(threw_desired_exception); + // If we have more than 4*4 processors, or a poor partitioner, we + // might not get an exception on every processor + mesh.comm().max(threw_desired_exception); - CPPUNIT_ASSERT_THROW_MESSAGE("Incompatible Elem/FE combo not detected", es.init(), libMesh::LogicError); + CPPUNIT_ASSERT(threw_desired_exception); } #endif From 73a2f32f1babd80ae8773619cef8104245235a8e Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 30 Jun 2022 14:52:35 -0500 Subject: [PATCH 08/13] Fix mesh stitching unit test setup The DistributedMesh invocation no longer fails on 9+ ranks. --- tests/mesh/mesh_stitch.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/mesh/mesh_stitch.C b/tests/mesh/mesh_stitch.C index 5d341b62111..d07e1d5d785 100644 --- a/tests/mesh/mesh_stitch.C +++ b/tests/mesh/mesh_stitch.C @@ -45,7 +45,7 @@ public: const std::string & boundary_name_prefix) { BoundaryInfo & boundary_info = mesh.get_boundary_info(); - const auto & mesh_boundary_ids = boundary_info.get_boundary_ids(); + const auto & mesh_boundary_ids = boundary_info.get_global_boundary_ids(); for (auto rit = mesh_boundary_ids.rbegin(); rit != mesh_boundary_ids.rend(); ++rit) { boundary_info.sideset_name(*rit + boundary_id_offset) = From 937dfb04bd83bfe8d0d6219af40f2b8e3cf55f19 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 30 Jun 2022 14:53:53 -0500 Subject: [PATCH 09/13] Clarifying comment --- tests/mesh/mesh_stitch.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/mesh/mesh_stitch.C b/tests/mesh/mesh_stitch.C index d07e1d5d785..a745f31fa00 100644 --- a/tests/mesh/mesh_stitch.C +++ b/tests/mesh/mesh_stitch.C @@ -86,6 +86,8 @@ public: const auto & nbi = bi.get_node_boundary_ids(); CPPUNIT_ASSERT_EQUAL(expected_size, nbi.size()); + // We expect that the "zero_right" and "one_left" boundaries have + // disappeared after being stitched together. std::set expected_names = {{"zero_left", "zero_top", "zero_front", From 613426e97737cf4e3515ddc5481995e4aa4ab402 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 30 Jun 2022 16:04:53 -0500 Subject: [PATCH 10/13] Fix ExactSolution w/ coarse DistMesh + fine mesh This isn't a good fix, because it precludes users carefully setting up ghosting manually between meshes... but since we're already serializing the solution here, we might as well serialize the mesh too. This fixes --enable-distmesh runs of miscellaneous_ex10 on 3+ processors. --- src/error_estimation/exact_solution.C | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/error_estimation/exact_solution.C b/src/error_estimation/exact_solution.C index 69771802b5a..b3780f4116d 100644 --- a/src/error_estimation/exact_solution.C +++ b/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 coarse_values; - std::unique_ptr> comparison_soln = NumericVector::build(_equation_systems.comm()); + std::unique_ptr> comparison_soln = + NumericVector::build(_equation_systems.comm()); + MeshSerializer + serial(const_cast(_equation_systems.get_mesh()), + _equation_systems_fine); if (_equation_systems_fine) { const System & comparison_system From e98bc0cc0687c4d0c4c2cce943debd362164b953 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 6 Jul 2022 09:19:07 -0500 Subject: [PATCH 11/13] More Poly2TriTriangulator robustness fixes Our interpretation of number of interpolation points didn't match our documentation or intentions, so some unit tests here had to be changed along with the code. Other than that, this just refactors us to be more robust to use cases like incoming meshes with weird node numberings. Poly2Tri is now passing all the tests I can throw at this in both MOOSE and libMesh. --- include/mesh/triangulator_interface.h | 8 ++- src/mesh/mesh_triangle_interface.C | 10 +++- src/mesh/poly2tri_triangulator.C | 17 ++++-- src/mesh/triangulator_interface.C | 77 ++++++++++++++++++--------- tests/mesh/mesh_triangulation.C | 16 +++--- 5 files changed, 87 insertions(+), 41 deletions(-) diff --git a/include/mesh/triangulator_interface.h b/include/mesh/triangulator_interface.h index 0f73e86b01c..b3b25ba3a3f 100644 --- a/include/mesh/triangulator_interface.h +++ b/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 diff --git a/src/mesh/mesh_triangle_interface.C b/src/mesh/mesh_triangle_interface.C index bf97de9bbc0..ba47ccc293a 100644 --- a/src/mesh/mesh_triangle_interface.C +++ b/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. diff --git a/src/mesh/poly2tri_triangulator.C b/src/mesh/poly2tri_triangulator.C index 76b95dda18a..8c8406f0b43 100644 --- a/src/mesh/poly2tri_triangulator.C +++ b/src/mesh/poly2tri_triangulator.C @@ -189,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(); @@ -204,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) @@ -307,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, @@ -586,6 +592,7 @@ bool Poly2TriTriangulator::insert_refinement_points() // 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()) { @@ -855,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; diff --git a/src/mesh/triangulator_interface.C b/src/mesh/triangulator_interface.C index e58078be04d..643d4d7162f 100644 --- a/src/mesh/triangulator_interface.C +++ b/src/mesh/triangulator_interface.C @@ -128,6 +128,35 @@ void TriangulatorInterface::elems_to_segments() +void TriangulatorInterface::nodes_to_segments(dof_id_type max_node_id) +{ + // Don't try to override manually specified segments, or try to add + // segments if we're doing a convex hull + if (!this->segments.empty() || _triangulation_type != PSLG) + return; + + for (auto node_it = _mesh.nodes_begin(), + node_end = _mesh.nodes_end(); + node_it != node_end;) + { + Node * node = *node_it; + + // If we're out of boundary nodes, the rest are going to be + // Steiner points or hole points + if (node->id() >= max_node_id) + break; + + ++node_it; + + Node * next_node = (node_it == node_end) ? + *_mesh.nodes_begin() : *node_it; + + this->segments.emplace_back(node->id(), next_node->id()); + } +} + + + void TriangulatorInterface::insert_any_extra_boundary_points() { // If the initial PSLG is really simple, e.g. an L-shaped domain or @@ -141,40 +170,38 @@ void TriangulatorInterface::insert_any_extra_boundary_points() const int n_interpolated = this->get_interpolate_boundary_points(); if ((_triangulation_type==PSLG) && n_interpolated) { - // Make a copy of the original points from the Mesh - std::vector original_points; - original_points.reserve (_mesh.n_nodes()); - for (auto & node : _mesh.node_ptr_range()) - original_points.push_back(*node); + std::vector> old_segments = + std::move(this->segments); - // Clear out the mesh - _mesh.clear(); + // We expect to have converted any elems and/or nodes into + // segments by now. + libmesh_assert(!old_segments.empty()); - // Make sure the new Mesh will be 2D - _mesh.set_mesh_dimension(2); + this->segments.clear(); - // Insert a new point on each PSLG at evenly spaced locations + // Insert a new point on each segment at evenly spaced locations // between existing boundary points. - // np=index into new points vector, also an explicit node index - // (to avoid the stride in automatic DistributedMesh indexing) + // np=index into new points vector // n =index into original points vector - for (std::size_t np=0, n=0, tops=n_interpolated*original_points.size(); npid(); + for (auto i : make_range(n_interpolated)) { - libmesh_assert_less(n-1, original_points.size()); - const std::size_t next_n = n % original_points.size(); + // new points are equispaced along the original segments const Point new_point = - (offset*original_points[next_n] + - (n_interpolated-offset)*original_points[n-1])/n_interpolated; - _mesh.add_point(new_point, np); + ((n_interpolated-i) * *(Point *)(begin_node) + + (i+1) * *(Point *)(end_node)) / + (n_interpolated + 1); + Node * next_node = _mesh.add_point(new_point); + this->segments.emplace_back(current_id, + next_node->id()); + current_id = next_node->id(); } + this->segments.emplace_back(current_id, + end_node->id()); } } } diff --git a/tests/mesh/mesh_triangulation.C b/tests/mesh/mesh_triangulation.C index d1785c95535..0c725f70109 100644 --- a/tests/mesh/mesh_triangulation.C +++ b/tests/mesh/mesh_triangulation.C @@ -32,7 +32,7 @@ public: #ifdef LIBMESH_HAVE_POLY2TRI CPPUNIT_TEST( testPoly2Tri ); CPPUNIT_TEST( testPoly2TriInterp ); - CPPUNIT_TEST( testPoly2TriInterp3 ); + CPPUNIT_TEST( testPoly2TriInterp2 ); CPPUNIT_TEST( testPoly2TriHoles ); CPPUNIT_TEST( testPoly2TriMeshedHoles ); CPPUNIT_TEST( testPoly2TriEdges ); @@ -54,7 +54,7 @@ public: #ifdef LIBMESH_HAVE_TRIANGLE CPPUNIT_TEST( testTriangle ); CPPUNIT_TEST( testTriangleInterp ); - CPPUNIT_TEST( testTriangleInterp3 ); + CPPUNIT_TEST( testTriangleInterp2 ); CPPUNIT_TEST( testTriangleHoles ); CPPUNIT_TEST( testTriangleMeshedHoles ); CPPUNIT_TEST( testTriangleSegments ); @@ -397,17 +397,17 @@ public: Mesh mesh(*TestCommWorld); TriangleInterface triangle(mesh); - testTriangulatorInterp(mesh, triangle, 2, 6); + testTriangulatorInterp(mesh, triangle, 1, 6); } - void testTriangleInterp3() + void testTriangleInterp2() { LOG_UNIT_TEST; Mesh mesh(*TestCommWorld); TriangleInterface triangle(mesh); - testTriangulatorInterp(mesh, triangle, 3, 10); + testTriangulatorInterp(mesh, triangle, 2, 10); } @@ -459,17 +459,17 @@ public: Mesh mesh(*TestCommWorld); Poly2TriTriangulator p2t_tri(mesh); - testTriangulatorInterp(mesh, p2t_tri, 2, 6); + testTriangulatorInterp(mesh, p2t_tri, 1, 6); } - void testPoly2TriInterp3() + void testPoly2TriInterp2() { LOG_UNIT_TEST; Mesh mesh(*TestCommWorld); Poly2TriTriangulator p2t_tri(mesh); - testTriangulatorInterp(mesh, p2t_tri, 3, 10); + testTriangulatorInterp(mesh, p2t_tri, 2, 10); } From 896f7d2cdc91aec55a48c09299c28826bb97ca37 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 5 Jul 2022 16:09:31 -0500 Subject: [PATCH 12/13] insert_any_extra_boundary_points contiguously This fixes Triangle unit tests on distributed meshes for me. ... which means we really ought to be finding whatever it is in our Triangle interface that doesn't support non-contiguous node numbering, but I couldn't find it on a quick glance, so "don't trigger the bug in my other code" is the best I can do right now. --- src/mesh/triangulator_interface.C | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/mesh/triangulator_interface.C b/src/mesh/triangulator_interface.C index 643d4d7162f..af7398469cd 100644 --- a/src/mesh/triangulator_interface.C +++ b/src/mesh/triangulator_interface.C @@ -170,6 +170,10 @@ void TriangulatorInterface::insert_any_extra_boundary_points() const int n_interpolated = this->get_interpolate_boundary_points(); if ((_triangulation_type==PSLG) && n_interpolated) { + // If we were lucky enough to start with contiguous node ids, + // let's keep them that way. + dof_id_type nn = _mesh.max_node_id(); + std::vector> old_segments = std::move(this->segments); @@ -195,7 +199,7 @@ void TriangulatorInterface::insert_any_extra_boundary_points() ((n_interpolated-i) * *(Point *)(begin_node) + (i+1) * *(Point *)(end_node)) / (n_interpolated + 1); - Node * next_node = _mesh.add_point(new_point); + Node * next_node = _mesh.add_point(new_point, nn++); this->segments.emplace_back(current_id, next_node->id()); current_id = next_node->id(); From 1e7ffca9c145b034e981259f6ad4b55b9ed36f07 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 5 Jul 2022 16:43:29 -0500 Subject: [PATCH 13/13] Don't spam cerr just because we're in 3D --- examples/subdomains/subdomains_ex3/subdomains_ex3.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/subdomains/subdomains_ex3/subdomains_ex3.C b/examples/subdomains/subdomains_ex3/subdomains_ex3.C index 019b932bb40..d75fc0dd85c 100644 --- a/examples/subdomains/subdomains_ex3/subdomains_ex3.C +++ b/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; }