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; } 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/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/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/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 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 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 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 d5ecd79d0d3..8c8406f0b43 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 @@ -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,13 +741,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 @@ -738,7 +755,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); @@ -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; diff --git a/src/mesh/triangulator_interface.C b/src/mesh/triangulator_interface.C index 7ae5a288270..af7398469cd 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,39 +170,42 @@ 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); + // 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); - // 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 // 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); + ((n_interpolated-i) * *(Point *)(begin_node) + + (i+1) * *(Point *)(end_node)) / + (n_interpolated + 1); + Node * next_node = _mesh.add_point(new_point, nn++); + 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/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 diff --git a/tests/mesh/mesh_stitch.C b/tests/mesh/mesh_stitch.C index 5d341b62111..a745f31fa00 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) = @@ -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", 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); } 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