diff --git a/examples/miscellaneous/miscellaneous_ex15/miscellaneous_ex15.C b/examples/miscellaneous/miscellaneous_ex15/miscellaneous_ex15.C index ad68b7ba9b..6394f16648 100644 --- a/examples/miscellaneous/miscellaneous_ex15/miscellaneous_ex15.C +++ b/examples/miscellaneous/miscellaneous_ex15/miscellaneous_ex15.C @@ -50,6 +50,8 @@ #include "libmesh/error_vector.h" #include "libmesh/mesh_refinement.h" #include "libmesh/kelly_error_estimator.h" +// for distributed mesh editing +#include "libmesh/parallel_ghost_sync.h" // Bring in everything from the libMesh namespace using namespace libMesh; @@ -65,7 +67,7 @@ void assemble_func(EquationSystems & es, const std::string & system_name) std::unique_ptr fe (FEBase::build(dim, fe_type)); std::unique_ptr inf_fe (FEBase::build_InfFE(dim, fe_type)); - const std::vector charged_objects = es.parameters.get >("charged_elem_id"); + const std::set charged_objects = es.parameters.get >("charged_elem_id"); const auto qrule=fe_type.default_quadrature_rule(/*dim = */ 3, /*extra order = */ 3); fe->attach_quadrature_rule (&*qrule); @@ -115,7 +117,7 @@ void assemble_func(EquationSystems & es, const std::string & system_name) Real rho; // check if charged_objects contains elem->id(), we add the corresponding charge to the rhs vector - if (std::find(charged_objects.begin(), charged_objects.end(), elem->id()) != charged_objects.end()) + if (charged_objects.count(elem->id())) rho=1./elem->volume(); else rho = 0; @@ -178,6 +180,10 @@ int main (int argc, char** argv) // creation of an empty mesh-object Mesh mesh(init.comm(), dim); + // We're going to keep track of elements by id, so we don't want + // even a DistributedMesh to renumber them. + mesh.allow_renumbering(false); + // fill the meshes with a spherical grid of type HEX27 with radius r MeshTools::Generation::build_cube (mesh, /*nx= */5, /*ny=*/5, /*nz=*/3, /*xmin=*/ -5.2, /*xmax=*/5.2, @@ -205,31 +211,44 @@ int main (int argc, char** argv) elem->neighbor_ptr(0)->subdomain_id()=2; } + // If we're on a distributed mesh then we might have missed some + // ghosted base elements with remote infinite elements. + if (!mesh.is_serial()) + { + SyncSubdomainIds sync_obj(mesh); + Parallel::sync_dofobject_data_by_id + (mesh.comm(), mesh.elements_begin(), mesh.elements_end(), + sync_obj); + } // Now we set the sources of the field: prism-shaped objects that are // determined here by containing certain points: - PointLocatorTree pt_lctr(mesh); - std::vector charged_elem_ids(3); + auto p_pt_lctr = mesh.sub_point_locator(); + auto & pt_lctr = *p_pt_lctr; + pt_lctr.enable_out_of_mesh_mode(); + + std::set charged_elem_ids; { Point pt_0(-3.,-3.0,-1.5); Point pt_1(2.,-2.6,-1.5); Point pt_2(2., 3.1, 1.7); const Elem * elem_0=pt_lctr(pt_0); if (elem_0) - charged_elem_ids[0]=elem_0->id(); - else - // this indicates some error which I don't know how to handle. - libmesh_not_implemented(); + charged_elem_ids.insert(elem_0->id()); const Elem * elem_1=pt_lctr(pt_1); if (elem_1) - charged_elem_ids[1]=elem_1->id(); - else - libmesh_not_implemented(); + charged_elem_ids.insert(elem_1->id()); const Elem * elem_2=pt_lctr(pt_2); if (elem_2) - charged_elem_ids[2]=elem_2->id(); - else - libmesh_not_implemented(); + charged_elem_ids.insert(elem_2->id()); + + // On a distributed mesh we might not have every point in a + // semilocal element on every processor + if (!mesh.is_serial()) + mesh.comm().set_union(charged_elem_ids); + + // But we should have every point on *some* processor + libmesh_assert_equal_to(charged_elem_ids.size(), 3); } // Create an equation systems object @@ -238,7 +257,7 @@ int main (int argc, char** argv) // This is the only system added here. LinearImplicitSystem & eig_sys = eq_sys.add_system ("Poisson"); - eq_sys.parameters.set >("charged_elem_id")=charged_elem_ids; + eq_sys.parameters.set >("charged_elem_id")=charged_elem_ids; //set the complete type of the variable FEType fe_type(SECOND, LAGRANGE, FOURTH, JACOBI_20_00, CARTESIAN); @@ -272,21 +291,31 @@ int main (int argc, char** argv) // in the refined mesh, find the elements that describe the // sources of gravitational field: Due to refinement, there are // successively more elements to account for the same object. - std::vector charged_elem_child(0); - for(unsigned int j=0; j< charged_elem_ids.size(); ++j) + std::set charged_elem_children; + for(auto id : charged_elem_ids) { - Elem * charged_elem = mesh.elem_ptr(charged_elem_ids[j]); + Elem * charged_elem = mesh.query_elem_ptr(id); + if (!charged_elem) + { + libmesh_assert(!mesh.is_serial()); + continue; + } + if (charged_elem->has_children()) { for(auto & child : charged_elem->child_ref_range()) - charged_elem_child.push_back(child.id()); + if (!child.is_remote()) + charged_elem_children.insert(child.id()); } else - charged_elem_child.push_back(charged_elem->id()); + charged_elem_children.insert(charged_elem->id()); } - charged_elem_ids=charged_elem_child; - eq_sys.parameters.set >("charged_elem_id")=charged_elem_ids; + charged_elem_ids=charged_elem_children; + if (!mesh.is_serial()) + mesh.comm().set_union(charged_elem_ids); + + eq_sys.parameters.set >("charged_elem_id")=charged_elem_ids; // re-assemble and than solve again. eig_sys.solve(); diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h index b1c40a44ee..818e3c05e5 100644 --- a/include/mesh/mesh_base.h +++ b/include/mesh/mesh_base.h @@ -837,14 +837,20 @@ class MeshBase : public ParallelObject * After this routine is called all the elements with a \p nullptr neighbor * pointer are guaranteed to be on the boundary. Thus this routine is * useful for automatically determining the boundaries of the domain. - * If reset_remote_elements is left to false, remote neighbor links are not - * reset and searched for in the local mesh. If reset_current_list is - * left as true, then any existing links will be reset before initiating - * the algorithm, while honoring the value of the reset_remote_elements - * flag. + * + * If \p reset_remote_elements is left to false, remote neighbor + * links are not reset and searched for in the local mesh. + * + * If \p reset_current_list is left as true, then any existing links + * will be reset before initiating the algorithm, while honoring the + * value of the \p reset_remote_elements flag. + * + * If \p assert_valid is left as true, then in dbg mode extensive + * consistency checking is performed before returning. */ virtual void find_neighbors (const bool reset_remote_elements = false, - const bool reset_current_list = true) = 0; + const bool reset_current_list = true, + const bool assert_valid = true) = 0; /** * Removes any orphaned nodes, nodes not connected to any elements. diff --git a/include/mesh/unstructured_mesh.h b/include/mesh/unstructured_mesh.h index be3c5798f4..53aa758499 100644 --- a/include/mesh/unstructured_mesh.h +++ b/include/mesh/unstructured_mesh.h @@ -282,7 +282,8 @@ class UnstructuredMesh : public MeshBase * Other functions from MeshBase requiring re-definition. */ virtual void find_neighbors (const bool reset_remote_elements = false, - const bool reset_current_list = true) override; + const bool reset_current_list = true, + const bool assert_valid = true) override; #ifdef LIBMESH_ENABLE_AMR /** diff --git a/src/mesh/inf_elem_builder.C b/src/mesh/inf_elem_builder.C index 2c919566d7..226b999161 100644 --- a/src/mesh/inf_elem_builder.C +++ b/src/mesh/inf_elem_builder.C @@ -35,6 +35,8 @@ #include "libmesh/mesh_base.h" #include "libmesh/remote_elem.h" +#include "timpi/parallel_implementation.h" + namespace libMesh { @@ -387,10 +389,14 @@ void InfElemBuilder::build_inf_elem(const Point & origin, } // neighbor(s) == nullptr - - - - + // On a distributed mesh it might be some other processor who sees + // the farthest node. + if (!this->_mesh.is_serial()) + { + unsigned int rank; + this->_mesh.comm().maxloc(max_r, rank); + this->_mesh.comm().broadcast(max_r_node, rank); + } // If a boundary side has one node on the outer boundary, // all points of this side are on the outer boundary. @@ -405,8 +411,14 @@ void InfElemBuilder::build_inf_elem(const Point & origin, // Only the case of no outer boundary is to be excluded. onodes.insert(max_r_node); + // If we're not on a serial mesh, we'll need to synchronize that + // onodes list too. + bool did_parallel_update; + do { + did_parallel_update = false; + auto face_it = faces.begin(); auto face_end = faces.end(); unsigned int facesfound=0; @@ -448,7 +460,16 @@ void InfElemBuilder::build_inf_elem(const Point & origin, face_it = faces.begin(); } } + + if (!this->_mesh.is_serial()) + { + auto my_onodes_size = onodes.size(); + this->_mesh.comm().set_union(onodes); + did_parallel_update = (onodes.size() > my_onodes_size); + this->_mesh.comm().max(did_parallel_update); + } } + while (did_parallel_update); #ifdef DEBUG @@ -486,18 +507,24 @@ void InfElemBuilder::build_inf_elem(const Point & origin, // for each boundary node, add an outer_node with // double distance from origin. - for (const auto & dof : onodes) + for (const auto & n : onodes) { - Point p = (Point(this->_mesh.point(dof)) * 2) - origin; + if (!this->_mesh.query_node_ptr(n)) + { + libmesh_assert(!_mesh.is_serial()); + continue; + } + + Point p = (Point(this->_mesh.point(n)) * 2) - origin; if (_mesh.is_serial()) { // Add with a default id in serial - outer_nodes[dof]=this->_mesh.add_point(p); + outer_nodes[n]=this->_mesh.add_point(p); } else { // Pick a unique id in parallel - Node & bnode = _mesh.node_ref(dof); + Node & bnode = _mesh.node_ref(n); dof_id_type new_id = bnode.id() + old_max_node_id; std::unique_ptr new_node = Node::build(p, new_id); new_node->processor_id() = bnode.processor_id(); @@ -505,7 +532,7 @@ void InfElemBuilder::build_inf_elem(const Point & origin, new_node->set_unique_id(old_max_unique_id + bnode.id()); #endif - outer_nodes[dof] = + outer_nodes[n] = this->_mesh.add_node(std::move(new_node)); } } @@ -665,6 +692,32 @@ void InfElemBuilder::build_inf_elem(const Point & origin, this->_mesh.add_elem(std::move(el)); } // for + // If we're not in serial, we might not have all our remote_elem + // neighbors set up properly on the new elements: two new infinite + // elements which should be neighbors may not be rooted in finite + // elmeents which are neighbors, and we were only looking for + // remote_elem links on each root finite element. We'll fix the + // problem by first finding all the neighbors we can and then + // recognizing that any missing neighbors on infinite elements must + // be remote. + if (!_mesh.is_serial()) + { + _mesh.find_neighbors(/*reset_remote_elements=*/ false, + /*reset_current_list=*/ true, + /*assert_valid=*/ false); + + for (auto & p : ofaces) + { + Elem & belem = this->_mesh.elem_ref(p.first); + Elem * inf_elem = belem.neighbor_ptr(p.second); + libmesh_assert(inf_elem); + + for (auto s : make_range(inf_elem->n_sides())) + if (!inf_elem->neighbor_ptr(s)) + inf_elem->set_neighbor + (s, const_cast(remote_elem)); + } + } #ifdef DEBUG _mesh.libmesh_assert_valid_parallel_ids(); diff --git a/src/mesh/unstructured_mesh.C b/src/mesh/unstructured_mesh.C index 935780b5dd..d7b952b724 100644 --- a/src/mesh/unstructured_mesh.C +++ b/src/mesh/unstructured_mesh.C @@ -916,7 +916,8 @@ UnstructuredMesh::~UnstructuredMesh () void UnstructuredMesh::find_neighbors (const bool reset_remote_elements, - const bool reset_current_list) + const bool reset_current_list, + const bool assert_valid) { // We might actually want to run this on an empty mesh // (e.g. the boundary mesh for a nonexistent bcid!) @@ -1303,9 +1304,14 @@ void UnstructuredMesh::find_neighbors (const bool reset_remote_elements, #ifdef DEBUG - MeshTools::libmesh_assert_valid_neighbors(*this, - !reset_remote_elements); - MeshTools::libmesh_assert_valid_amr_interior_parents(*this); + if (assert_valid) + { + MeshTools::libmesh_assert_valid_neighbors(*this, + !reset_remote_elements); + MeshTools::libmesh_assert_valid_amr_interior_parents(*this); + } +#else + libmesh_ignore(assert_valid); #endif }