Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // library configuration
21 : #include "libmesh/libmesh_config.h"
22 :
23 : // Local includes
24 : #include "libmesh/boundary_info.h"
25 : #include "libmesh/libmesh_logging.h"
26 : #include "libmesh/elem.h"
27 : #include "libmesh/ghost_point_neighbors.h"
28 : #include "libmesh/mesh_base.h"
29 : #include "libmesh/mesh_communication.h"
30 : #include "libmesh/mesh_serializer.h"
31 : #include "libmesh/mesh_tools.h"
32 : #include "libmesh/parallel.h"
33 : #include "libmesh/parallel_algebra.h"
34 : #include "libmesh/parallel_fe_type.h"
35 : #include "libmesh/partitioner.h"
36 : #include "libmesh/point_locator_base.h"
37 : #include "libmesh/sparse_matrix.h"
38 : #include "libmesh/threads.h"
39 : #include "libmesh/enum_elem_type.h"
40 : #include "libmesh/enum_point_locator_type.h"
41 : #include "libmesh/enum_to_string.h"
42 : #include "libmesh/point_locator_nanoflann.h"
43 : #include "libmesh/elem_side_builder.h"
44 : #include "libmesh/elem_range.h"
45 : #include "libmesh/node_range.h"
46 :
47 : // C++ includes
48 : #include <algorithm> // for std::min
49 : #include <map> // for std::multimap
50 : #include <memory>
51 : #include <sstream> // for std::ostringstream
52 : #include <unordered_map>
53 :
54 : #include "libmesh/periodic_boundaries.h"
55 : #include "libmesh/periodic_boundary.h"
56 :
57 : namespace libMesh
58 : {
59 :
60 :
61 :
62 : // ------------------------------------------------------------
63 : // MeshBase class member functions
64 33396 : MeshBase::MeshBase (const Parallel::Communicator & comm_in,
65 33396 : unsigned char d) :
66 : ParallelObject (comm_in),
67 33396 : boundary_info (new BoundaryInfo(*this)), // BoundaryInfo has protected ctor, can't use std::make_unique
68 14080 : _n_parts (1),
69 14080 : _default_mapping_type(LAGRANGE_MAP),
70 14080 : _default_mapping_data(0),
71 14080 : _preparation (),
72 14080 : _element_stored_range (),
73 14080 : _const_active_local_element_stored_range (),
74 14080 : _point_locator (),
75 14080 : _count_lower_dim_elems_in_point_locator(true),
76 14080 : _partitioner (),
77 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
78 14080 : _next_unique_id(DofObject::invalid_unique_id),
79 : #endif
80 14080 : _interior_mesh(this),
81 14080 : _skip_noncritical_partitioning(false),
82 43054 : _skip_all_partitioning(libMesh::on_command_line("--skip-partitioning")),
83 14080 : _skip_renumber_nodes_and_elements(false),
84 14080 : _skip_find_neighbors(false),
85 14080 : _skip_detect_interior_parents(false),
86 14080 : _allow_remote_element_removal(true),
87 14080 : _allow_node_and_elem_unique_id_overlap(false),
88 14080 : _spatial_dimension(d),
89 33396 : _default_ghosting(std::make_unique<GhostPointNeighbors>(*this)),
90 148478 : _point_locator_close_to_point_tol(0.)
91 : {
92 23738 : _elem_dims.insert(d);
93 33396 : _ghosting_functors.push_back(_default_ghosting.get());
94 9658 : libmesh_assert_less_equal (LIBMESH_DIM, 3);
95 9658 : libmesh_assert_greater_equal (LIBMESH_DIM, d);
96 9658 : libmesh_assert (libMesh::initialized());
97 33396 : }
98 :
99 :
100 :
101 2632 : MeshBase::MeshBase (const MeshBase & other_mesh) :
102 : ParallelObject (other_mesh),
103 2632 : boundary_info (new BoundaryInfo(*this)), // BoundaryInfo has protected ctor, can't use std::make_unique
104 2632 : _n_parts (other_mesh._n_parts),
105 2632 : _default_mapping_type(other_mesh._default_mapping_type),
106 2632 : _default_mapping_data(other_mesh._default_mapping_data),
107 1040 : _preparation (other_mesh._preparation),
108 1040 : _element_stored_range (),
109 1040 : _const_active_local_element_stored_range (),
110 1040 : _point_locator (),
111 2632 : _count_lower_dim_elems_in_point_locator(other_mesh._count_lower_dim_elems_in_point_locator),
112 1040 : _partitioner (),
113 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
114 2632 : _next_unique_id(other_mesh._next_unique_id),
115 : #endif
116 : // If the other mesh interior_parent pointers just go back to
117 : // itself, so should we
118 2632 : _interior_mesh((other_mesh._interior_mesh == &other_mesh) ?
119 : this : other_mesh._interior_mesh),
120 2632 : _skip_noncritical_partitioning(other_mesh._skip_noncritical_partitioning),
121 2632 : _skip_all_partitioning(other_mesh._skip_all_partitioning),
122 2632 : _skip_renumber_nodes_and_elements(other_mesh._skip_renumber_nodes_and_elements),
123 2632 : _skip_find_neighbors(other_mesh._skip_find_neighbors),
124 2632 : _skip_detect_interior_parents(other_mesh._skip_detect_interior_parents),
125 2632 : _allow_remote_element_removal(other_mesh._allow_remote_element_removal),
126 2632 : _allow_node_and_elem_unique_id_overlap(other_mesh._allow_node_and_elem_unique_id_overlap),
127 804 : _elem_dims(other_mesh._elem_dims),
128 804 : _elem_default_orders(other_mesh._elem_default_orders),
129 2632 : _supported_nodal_order(other_mesh._supported_nodal_order),
130 804 : _mesh_subdomains(other_mesh._mesh_subdomains),
131 804 : _elemset_codes_inverse_map(other_mesh._elemset_codes_inverse_map),
132 804 : _all_elemset_ids(other_mesh._all_elemset_ids),
133 2632 : _spatial_dimension(other_mesh._spatial_dimension),
134 2632 : _default_ghosting(std::make_unique<GhostPointNeighbors>(*this)),
135 12892 : _point_locator_close_to_point_tol(other_mesh._point_locator_close_to_point_tol)
136 : {
137 804 : const GhostingFunctor * const other_default_ghosting = other_mesh._default_ghosting.get();
138 :
139 8556 : for (GhostingFunctor * const gf : other_mesh._ghosting_functors)
140 : {
141 : // If the other mesh is using default ghosting, then we will use our own
142 : // default ghosting
143 5924 : if (gf == other_default_ghosting)
144 : {
145 2632 : _ghosting_functors.push_back(_default_ghosting.get());
146 2632 : continue;
147 : }
148 :
149 5568 : std::shared_ptr<GhostingFunctor> clone_gf = gf->clone();
150 : // Some subclasses of GhostingFunctor might not override the
151 : // clone function yet. If this is the case, GhostingFunctor will
152 : // return nullptr by default. The clone function should be overridden
153 : // in all derived classes. This following code ("else") is written
154 : // for API upgrade. That will allow users gradually to update their code.
155 : // Once the API upgrade is done, we will come back and delete "else."
156 3292 : if (clone_gf)
157 : {
158 3292 : clone_gf->set_mesh(this);
159 5568 : add_ghosting_functor(clone_gf);
160 : }
161 : else
162 : {
163 : libmesh_deprecated();
164 0 : add_ghosting_functor(*gf);
165 : }
166 : }
167 :
168 2632 : if (other_mesh._partitioner.get())
169 4460 : _partitioner = other_mesh._partitioner->clone();
170 :
171 : #ifdef LIBMESH_ENABLE_PERIODIC
172 : // Deep copy of all periodic boundaries
173 2632 : if (other_mesh._disjoint_neighbor_boundary_pairs)
174 : {
175 8 : _disjoint_neighbor_boundary_pairs = std::make_unique<PeriodicBoundaries>();
176 :
177 21 : for (const auto & [id, pb] : *other_mesh._disjoint_neighbor_boundary_pairs)
178 14 : if (pb)
179 24 : (*_disjoint_neighbor_boundary_pairs)[id] = pb->clone();
180 : }
181 : #endif
182 :
183 : // _elemset_codes stores pointers to entries in _elemset_codes_inverse_map,
184 : // so it is not possible to simply copy it directly from other_mesh
185 2632 : for (const auto & [set, code] : _elemset_codes_inverse_map)
186 0 : _elemset_codes.emplace(code, &set);
187 2632 : }
188 :
189 42 : MeshBase& MeshBase::operator= (MeshBase && other_mesh)
190 : {
191 12 : LOG_SCOPE("operator=(&&)", "MeshBase");
192 :
193 : // Move assign as a ParallelObject.
194 12 : this->ParallelObject::operator=(other_mesh);
195 :
196 42 : _n_parts = other_mesh.n_partitions();
197 42 : _default_mapping_type = other_mesh.default_mapping_type();
198 42 : _default_mapping_data = other_mesh.default_mapping_data();
199 42 : _preparation = other_mesh._preparation;
200 12 : _element_stored_range = std::move(other_mesh._element_stored_range);
201 12 : _const_active_local_element_stored_range = std::move(other_mesh._const_active_local_element_stored_range);
202 12 : _point_locator = std::move(other_mesh._point_locator);
203 42 : _count_lower_dim_elems_in_point_locator = other_mesh.get_count_lower_dim_elems_in_point_locator();
204 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
205 42 : _next_unique_id = other_mesh.next_unique_id();
206 : #endif
207 : // If the other mesh interior_parent pointers just go back to
208 : // itself, so should we
209 42 : _interior_mesh = (other_mesh._interior_mesh == &other_mesh) ?
210 : this : other_mesh._interior_mesh;
211 42 : _skip_noncritical_partitioning = other_mesh.skip_noncritical_partitioning();
212 42 : _skip_all_partitioning = other_mesh.skip_partitioning();
213 42 : _skip_renumber_nodes_and_elements = !(other_mesh.allow_renumbering());
214 42 : _skip_find_neighbors = !(other_mesh.allow_find_neighbors());
215 42 : _skip_detect_interior_parents = !(other_mesh.allow_detect_interior_parents());
216 42 : _allow_remote_element_removal = other_mesh.allow_remote_element_removal();
217 42 : _allow_node_and_elem_unique_id_overlap = other_mesh.allow_node_and_elem_unique_id_overlap();
218 12 : _block_id_to_name = std::move(other_mesh._block_id_to_name);
219 12 : _elem_dims = std::move(other_mesh.elem_dimensions());
220 12 : _elem_default_orders = std::move(other_mesh.elem_default_orders());
221 42 : _supported_nodal_order = other_mesh.supported_nodal_order();
222 12 : _mesh_subdomains = other_mesh._mesh_subdomains;
223 12 : _elemset_codes = std::move(other_mesh._elemset_codes);
224 12 : _elemset_codes_inverse_map = std::move(other_mesh._elemset_codes_inverse_map);
225 12 : _all_elemset_ids = std::move(other_mesh._all_elemset_ids);
226 42 : _spatial_dimension = other_mesh.spatial_dimension();
227 42 : _elem_integer_names = std::move(other_mesh._elem_integer_names);
228 42 : _elem_integer_default_values = std::move(other_mesh._elem_integer_default_values);
229 42 : _node_integer_names = std::move(other_mesh._node_integer_names);
230 42 : _node_integer_default_values = std::move(other_mesh._node_integer_default_values);
231 42 : _point_locator_close_to_point_tol = other_mesh.get_point_locator_close_to_point_tol();
232 :
233 : #ifdef LIBMESH_ENABLE_PERIODIC
234 : // Deep copy of all periodic boundaries:
235 : // We must clone each PeriodicBoundaryBase in the source map,
236 : // since unique_ptr cannot be copied and we need independent instances
237 42 : if (other_mesh._disjoint_neighbor_boundary_pairs)
238 : {
239 0 : _disjoint_neighbor_boundary_pairs = std::make_unique<PeriodicBoundaries>();
240 :
241 0 : for (const auto & [id, pb] : *other_mesh._disjoint_neighbor_boundary_pairs)
242 0 : if (pb)
243 0 : (*_disjoint_neighbor_boundary_pairs)[id] = pb->clone();
244 : }
245 : #endif
246 :
247 : // This relies on our subclasses *not* invalidating pointers when we
248 : // do their portion of the move assignment later!
249 12 : boundary_info = std::move(other_mesh.boundary_info);
250 12 : boundary_info->set_mesh(*this);
251 :
252 : #ifdef DEBUG
253 : // Make sure that move assignment worked for pointers
254 12 : for (const auto & [set, code] : _elemset_codes_inverse_map)
255 : {
256 0 : auto it = _elemset_codes.find(code);
257 0 : libmesh_assert_msg(it != _elemset_codes.end(),
258 : "Elemset code " << code << " not found in _elmset_codes container.");
259 0 : libmesh_assert_equal_to(it->second, &set);
260 : }
261 : #endif
262 :
263 : // We're *not* really done at this point, but we have the problem
264 : // that some of our data movement might be expecting subclasses data
265 : // movement to happen first. We'll let subclasses handle that by
266 : // calling our post_dofobject_moves()
267 54 : return *this;
268 : }
269 :
270 :
271 2095 : bool MeshBase::operator== (const MeshBase & other_mesh) const
272 : {
273 1062 : LOG_SCOPE("operator==()", "MeshBase");
274 :
275 2095 : bool is_equal = this->locally_equals(other_mesh);
276 2095 : this->comm().min(is_equal);
277 3157 : return is_equal;
278 : }
279 :
280 :
281 2095 : bool MeshBase::locally_equals (const MeshBase & other_mesh) const
282 : {
283 : // Check whether (almost) everything in the base is equal
284 : //
285 : // We don't check _next_unique_id here, because it's expected to
286 : // change in a DistributedMesh prepare_for_use(); it's conceptually
287 : // "mutable".
288 : //
289 : // We use separate if statements instead of logical operators here,
290 : // to make it easy to see the failing condition when using a
291 : // debugger to figure out why a MeshTools::valid_is_prepared(mesh)
292 : // is failing.
293 2095 : if (_n_parts != other_mesh._n_parts)
294 0 : return false;
295 2095 : if (_default_mapping_type != other_mesh._default_mapping_type)
296 0 : return false;
297 2095 : if (_default_mapping_data != other_mesh._default_mapping_data)
298 0 : return false;
299 2095 : if (_preparation != other_mesh._preparation)
300 0 : return false;
301 3565 : if (_count_lower_dim_elems_in_point_locator !=
302 2095 : other_mesh._count_lower_dim_elems_in_point_locator)
303 0 : return false;
304 :
305 : // We should either both have our own interior parents or both not;
306 : // but if we both don't then we can't really assert anything else
307 : // because pointing at the same interior mesh is fair but so is
308 : // pointing at two different copies of "the same" interior mesh.
309 3565 : if ((_interior_mesh == this) !=
310 2095 : (other_mesh._interior_mesh == &other_mesh))
311 0 : return false;
312 :
313 2095 : if (_skip_noncritical_partitioning != other_mesh._skip_noncritical_partitioning)
314 0 : return false;
315 2095 : if (_skip_all_partitioning != other_mesh._skip_all_partitioning)
316 0 : return false;
317 2095 : if (_skip_renumber_nodes_and_elements != other_mesh._skip_renumber_nodes_and_elements)
318 0 : return false;
319 2095 : if (_skip_find_neighbors != other_mesh._skip_find_neighbors)
320 0 : return false;
321 2095 : if (_skip_detect_interior_parents != other_mesh._skip_detect_interior_parents)
322 0 : return false;
323 2095 : if (_allow_remote_element_removal != other_mesh._allow_remote_element_removal)
324 0 : return false;
325 2095 : if (_allow_node_and_elem_unique_id_overlap != other_mesh._allow_node_and_elem_unique_id_overlap)
326 0 : return false;
327 2095 : if (_spatial_dimension != other_mesh._spatial_dimension)
328 0 : return false;
329 2095 : if (_point_locator_close_to_point_tol != other_mesh._point_locator_close_to_point_tol)
330 0 : return false;
331 2095 : if (_block_id_to_name != other_mesh._block_id_to_name)
332 0 : return false;
333 2095 : if (_elem_dims != other_mesh._elem_dims)
334 0 : return false;
335 2095 : if (_elem_default_orders != other_mesh._elem_default_orders)
336 0 : return false;
337 2095 : if (_supported_nodal_order != other_mesh._supported_nodal_order)
338 0 : return false;
339 2095 : if (_mesh_subdomains != other_mesh._mesh_subdomains)
340 6 : return false;
341 2074 : if (_all_elemset_ids != other_mesh._all_elemset_ids)
342 0 : return false;
343 2074 : if (_elem_integer_names != other_mesh._elem_integer_names)
344 0 : return false;
345 2476 : if (_elem_integer_default_values != other_mesh._elem_integer_default_values)
346 0 : return false;
347 2074 : if (_node_integer_names != other_mesh._node_integer_names)
348 0 : return false;
349 2476 : if (_node_integer_default_values != other_mesh._node_integer_default_values)
350 0 : return false;
351 2074 : if (static_cast<bool>(_default_ghosting) != static_cast<bool>(other_mesh._default_ghosting))
352 0 : return false;
353 2074 : if (static_cast<bool>(_partitioner) != static_cast<bool>(other_mesh._partitioner))
354 0 : return false;
355 2074 : if (*boundary_info != *other_mesh.boundary_info)
356 6 : return false;
357 :
358 : // First check whether the "existence" of the two pointers differs (one present, one absent)
359 3103 : if (static_cast<bool>(_disjoint_neighbor_boundary_pairs) !=
360 1050 : static_cast<bool>(other_mesh._disjoint_neighbor_boundary_pairs))
361 0 : return false;
362 : // If both exist, compare the contents (Weak Test: just compare sizes like `_ghosting_functors`)
363 2057 : if (_disjoint_neighbor_boundary_pairs &&
364 1054 : (_disjoint_neighbor_boundary_pairs->size() != other_mesh._disjoint_neighbor_boundary_pairs->size()))
365 0 : return false;
366 :
367 : const constraint_rows_type & other_rows =
368 1050 : other_mesh.get_constraint_rows();
369 17308 : for (const auto & [node, row] : this->_constraint_rows)
370 : {
371 15255 : const dof_id_type node_id = node->id();
372 15255 : const Node * other_node = other_mesh.query_node_ptr(node_id);
373 15255 : if (!other_node)
374 0 : return false;
375 :
376 6416 : auto it = other_rows.find(other_node);
377 15255 : if (it == other_rows.end())
378 0 : return false;
379 :
380 6416 : const auto & other_row = it->second;
381 21671 : if (row.size() != other_row.size())
382 0 : return false;
383 :
384 56494 : for (auto i : index_range(row))
385 : {
386 18096 : const auto & [elem_pair, coef] = row[i];
387 18096 : const auto & [other_elem_pair, other_coef] = other_row[i];
388 18096 : libmesh_assert(elem_pair.first);
389 18096 : libmesh_assert(other_elem_pair.first);
390 41239 : if (elem_pair.first->id() !=
391 64382 : other_elem_pair.first->id() ||
392 41239 : elem_pair.second !=
393 59335 : other_elem_pair.second ||
394 41239 : coef != other_coef)
395 0 : return false;
396 : }
397 : }
398 :
399 2053 : for (const auto & [elemset_code, elemset_ptr] : this->_elemset_codes)
400 0 : if (const auto it = other_mesh._elemset_codes.find(elemset_code);
401 0 : it == other_mesh._elemset_codes.end() || *elemset_ptr != *it->second)
402 0 : return false;
403 :
404 : // FIXME: we have no good way to compare ghosting functors, since
405 : // they're in a vector of pointers, and we have no way *at all*
406 : // to compare ghosting functors, since they don't have operator==
407 : // defined and we encourage users to subclass them. We can check if
408 : // we have the same number, is all.
409 3499 : if (_ghosting_functors.size() !=
410 1446 : other_mesh._ghosting_functors.size())
411 0 : return false;
412 :
413 : // Same deal for partitioners. We tested that we both have one or
414 : // both don't, but are they equivalent? Let's guess "yes".
415 :
416 : // Now let the subclasses decide whether everything else is equal
417 2053 : return this->subclass_locally_equals(other_mesh);
418 : }
419 :
420 :
421 46474 : MeshBase::~MeshBase()
422 : {
423 36028 : this->MeshBase::clear();
424 :
425 10462 : libmesh_exceptionless_assert (!libMesh::closed());
426 81388 : }
427 :
428 :
429 :
430 2967994 : unsigned int MeshBase::mesh_dimension() const
431 : {
432 2967994 : if (!_elem_dims.empty())
433 2967994 : return cast_int<unsigned int>(*_elem_dims.rbegin());
434 0 : return 0;
435 : }
436 :
437 :
438 :
439 0 : void MeshBase::set_elem_dimensions(std::set<unsigned char> elem_dims)
440 : {
441 : #ifdef DEBUG
442 : // In debug mode, we call cache_elem_data() and then make sure
443 : // the result actually agrees with what the user specified.
444 0 : parallel_object_only();
445 :
446 0 : this->cache_elem_data();
447 0 : libmesh_assert_msg(_elem_dims == elem_dims, \
448 : "Specified element dimensions does not match true element dimensions!");
449 : #endif
450 :
451 0 : _elem_dims = std::move(elem_dims);
452 0 : }
453 :
454 :
455 :
456 175 : void MeshBase::add_elemset_code(dof_id_type code, MeshBase::elemset_type id_set)
457 : {
458 : // Populate inverse map, stealing id_set's resources
459 125 : auto [it1, inserted1] = _elemset_codes_inverse_map.emplace(std::move(id_set), code);
460 :
461 : // Reference to the newly inserted (or previously existing) id_set
462 175 : const auto & inserted_id_set = it1->first;
463 :
464 : // Keep track of all elemset ids ever added for O(1) n_elemsets()
465 : // performance. Only need to do this if we didn't know about this
466 : // id_set before...
467 175 : if (inserted1)
468 175 : _all_elemset_ids.insert(inserted_id_set.begin(), inserted_id_set.end());
469 :
470 : // Take the address of the newly emplaced set to use in
471 : // _elemset_codes, avoid duplicating std::set storage
472 175 : auto [it2, inserted2] = _elemset_codes.emplace(code, &inserted_id_set);
473 :
474 : // Throw an error if this code already exists with a pointer to a
475 : // different set of ids.
476 175 : libmesh_error_msg_if(!inserted2 && it2->second != &inserted_id_set,
477 : "The elemset code " << code << " already exists with a different id_set.");
478 175 : }
479 :
480 :
481 :
482 1431 : unsigned int MeshBase::n_elemsets() const
483 : {
484 1431 : return _all_elemset_ids.size();
485 : }
486 :
487 371 : void MeshBase::get_elemsets(dof_id_type elemset_code, MeshBase::elemset_type & id_set_to_fill) const
488 : {
489 : // If we don't recognize this elemset_code, hand back an empty set
490 102 : id_set_to_fill.clear();
491 :
492 371 : if (const auto it = _elemset_codes.find(elemset_code);
493 102 : it != _elemset_codes.end())
494 162 : id_set_to_fill.insert(it->second->begin(), it->second->end());
495 371 : }
496 :
497 84 : dof_id_type MeshBase::get_elemset_code(const MeshBase::elemset_type & id_set) const
498 : {
499 24 : auto it = _elemset_codes_inverse_map.find(id_set);
500 84 : return (it == _elemset_codes_inverse_map.end()) ? DofObject::invalid_id : it->second;
501 : }
502 :
503 837 : std::vector<dof_id_type> MeshBase::get_elemset_codes() const
504 : {
505 240 : std::vector<dof_id_type> ret;
506 837 : ret.reserve(_elemset_codes.size());
507 914 : for (const auto & pr : _elemset_codes)
508 77 : ret.push_back(pr.first);
509 837 : return ret;
510 : }
511 :
512 28 : void MeshBase::change_elemset_code(dof_id_type old_code, dof_id_type new_code)
513 : {
514 : // Look up elemset ids for old_code
515 8 : auto it = _elemset_codes.find(old_code);
516 :
517 : // If we don't have the old_code, then do nothing. Alternatively, we
518 : // could throw an error since trying to change an elemset code you
519 : // don't have could indicate there's a problem...
520 28 : if (it == _elemset_codes.end())
521 0 : return;
522 :
523 : // Make copy of the set of elemset ids. We are not changing these,
524 : // only updating the elemset code it corresponds to.
525 28 : elemset_type id_set_copy = *(it->second);
526 :
527 : // Look up the corresponding entry in the inverse map. Note: we want
528 : // the iterator because we are going to remove it.
529 8 : auto inverse_it = _elemset_codes_inverse_map.find(id_set_copy);
530 28 : libmesh_error_msg_if(inverse_it == _elemset_codes_inverse_map.end(),
531 : "Expected _elemset_codes_inverse_map entry for elemset code " << old_code);
532 :
533 : // Erase entry from inverse map
534 20 : _elemset_codes_inverse_map.erase(inverse_it);
535 :
536 : // Erase entry from forward map
537 20 : _elemset_codes.erase(it);
538 :
539 : // Add new code with original set of ids.
540 48 : this->add_elemset_code(new_code, id_set_copy);
541 :
542 : // We can't update any actual elemset codes if there is no extra integer defined for it.
543 28 : if (!this->has_elem_integer("elemset_code"))
544 0 : return;
545 :
546 : // Get index of elemset_code extra integer
547 28 : unsigned int elemset_index = this->get_elem_integer_index("elemset_code");
548 :
549 : // Loop over all elems and update code
550 : Threads::parallel_for
551 28 : (this->element_stored_range(),
552 472 : [elemset_index, old_code, new_code](const ElemRange & range)
553 : {
554 518 : for (Elem * elem : range)
555 : {
556 : dof_id_type elemset_code =
557 490 : elem->get_extra_integer(elemset_index);
558 :
559 490 : if (elemset_code == old_code)
560 175 : elem->set_extra_integer(elemset_index, new_code);
561 : }
562 28 : });
563 : }
564 :
565 28 : void MeshBase::change_elemset_id(elemset_id_type old_id, elemset_id_type new_id)
566 : {
567 : // Early return if we don't have old_id
568 8 : if (!_all_elemset_ids.count(old_id))
569 0 : return;
570 :
571 : // Throw an error if the new_id is already used
572 8 : libmesh_error_msg_if(_all_elemset_ids.count(new_id),
573 : "Cannot change elemset id " << old_id <<
574 : " to " << new_id << ", " << new_id << " already exists.");
575 :
576 : // We will build up a new version of the inverse map so we can iterate over
577 : // the current one without invalidating anything.
578 16 : std::map<MeshBase::elemset_type, dof_id_type> new_elemset_codes_inverse_map;
579 84 : for (const auto & [id_set, elemset_code] : _elemset_codes_inverse_map)
580 : {
581 32 : auto id_set_copy = id_set;
582 16 : if (id_set_copy.count(old_id))
583 : {
584 : // Remove old_id, insert new_id
585 8 : id_set_copy.erase(old_id);
586 20 : id_set_copy.insert(new_id);
587 : }
588 :
589 : // Store in new version of map
590 40 : new_elemset_codes_inverse_map.emplace(id_set_copy, elemset_code);
591 : }
592 :
593 : // Swap existing map with newly-built one
594 8 : _elemset_codes_inverse_map.swap(new_elemset_codes_inverse_map);
595 :
596 : // Reconstruct _elemset_codes map
597 8 : _elemset_codes.clear();
598 84 : for (const auto & [id_set, elemset_code] : _elemset_codes_inverse_map)
599 56 : _elemset_codes.emplace(elemset_code, &id_set);
600 :
601 : // Update _all_elemset_ids
602 8 : _all_elemset_ids.erase(old_id);
603 20 : _all_elemset_ids.insert(new_id);
604 : }
605 :
606 27878 : unsigned int MeshBase::spatial_dimension () const
607 : {
608 36582 : return cast_int<unsigned int>(_spatial_dimension);
609 : }
610 :
611 :
612 :
613 36416 : void MeshBase::set_spatial_dimension(unsigned char d)
614 : {
615 : // The user can set the _spatial_dimension however they wish,
616 : // libMesh will only *increase* the spatial dimension, however,
617 : // never decrease it.
618 36416 : _spatial_dimension = d;
619 36416 : }
620 :
621 :
622 :
623 492 : unsigned int MeshBase::add_elem_integer(std::string name,
624 : bool allocate_data,
625 : dof_id_type default_value)
626 : {
627 664 : for (auto i : index_range(_elem_integer_names))
628 136 : if (_elem_integer_names[i] == name)
629 : {
630 8 : libmesh_assert_less(i, _elem_integer_default_values.size());
631 28 : _elem_integer_default_values[i] = default_value;
632 28 : return i;
633 : }
634 :
635 136 : libmesh_assert_equal_to(_elem_integer_names.size(),
636 : _elem_integer_default_values.size());
637 464 : _elem_integer_names.push_back(std::move(name));
638 464 : _elem_integer_default_values.push_back(default_value);
639 464 : if (allocate_data)
640 388 : this->size_elem_extra_integers();
641 600 : return _elem_integer_names.size()-1;
642 : }
643 :
644 :
645 :
646 4021 : std::vector<unsigned int> MeshBase::add_elem_integers(const std::vector<std::string> & names,
647 : bool allocate_data,
648 : const std::vector<dof_id_type> * default_values)
649 : {
650 1208 : libmesh_assert(!default_values || default_values->size() == names.size());
651 1208 : libmesh_assert_equal_to(_elem_integer_names.size(), _elem_integer_default_values.size());
652 :
653 2416 : std::unordered_map<std::string, std::size_t> name_indices;
654 4105 : for (auto i : index_range(_elem_integer_names))
655 84 : name_indices[_elem_integer_names[i]] = i;
656 :
657 5213 : std::vector<unsigned int> returnval(names.size());
658 :
659 1208 : bool added_an_integer = false;
660 4286 : for (auto i : index_range(names))
661 : {
662 156 : const std::string & name = names[i];
663 265 : if (const auto it = name_indices.find(name);
664 78 : it != name_indices.end())
665 : {
666 42 : returnval[i] = it->second;
667 54 : _elem_integer_default_values[it->second] =
668 42 : default_values ? (*default_values)[i] : DofObject::invalid_id;
669 : }
670 : else
671 : {
672 289 : returnval[i] = _elem_integer_names.size();
673 223 : name_indices[name] = returnval[i];
674 223 : _elem_integer_names.push_back(name);
675 : _elem_integer_default_values.push_back
676 259 : (default_values ? (*default_values)[i] : DofObject::invalid_id);
677 66 : added_an_integer = true;
678 : }
679 : }
680 :
681 4021 : if (allocate_data && added_an_integer)
682 105 : this->size_elem_extra_integers();
683 :
684 5229 : return returnval;
685 : }
686 :
687 :
688 :
689 326 : unsigned int MeshBase::get_elem_integer_index(std::string_view name) const
690 : {
691 410 : for (auto i : index_range(_elem_integer_names))
692 317 : if (_elem_integer_names[i] == name)
693 326 : return i;
694 :
695 0 : libmesh_error_msg("Unknown elem integer " << name);
696 : return libMesh::invalid_uint;
697 : }
698 :
699 :
700 :
701 2599 : bool MeshBase::has_elem_integer(std::string_view name) const
702 : {
703 2697 : for (auto & entry : _elem_integer_names)
704 239 : if (entry == name)
705 40 : return true;
706 :
707 632 : return false;
708 : }
709 :
710 :
711 :
712 2284 : unsigned int MeshBase::add_node_integer(std::string name,
713 : bool allocate_data,
714 : dof_id_type default_value)
715 : {
716 4441 : for (auto i : index_range(_node_integer_names))
717 1440 : if (_node_integer_names[i] == name)
718 : {
719 16 : libmesh_assert_less(i, _node_integer_default_values.size());
720 56 : _node_integer_default_values[i] = default_value;
721 56 : return i;
722 : }
723 :
724 702 : libmesh_assert_equal_to(_node_integer_names.size(),
725 : _node_integer_default_values.size());
726 2228 : _node_integer_names.push_back(std::move(name));
727 2228 : _node_integer_default_values.push_back(default_value);
728 2228 : if (allocate_data)
729 720 : this->size_node_extra_integers();
730 2930 : return _node_integer_names.size()-1;
731 : }
732 :
733 :
734 :
735 4007 : std::vector<unsigned int> MeshBase::add_node_integers(const std::vector<std::string> & names,
736 : bool allocate_data,
737 : const std::vector<dof_id_type> * default_values)
738 : {
739 1204 : libmesh_assert(!default_values || default_values->size() == names.size());
740 1204 : libmesh_assert_equal_to(_node_integer_names.size(), _node_integer_default_values.size());
741 :
742 2408 : std::unordered_map<std::string, std::size_t> name_indices;
743 4063 : for (auto i : index_range(_node_integer_names))
744 56 : name_indices[_node_integer_names[i]] = i;
745 :
746 5195 : std::vector<unsigned int> returnval(names.size());
747 :
748 1204 : bool added_an_integer = false;
749 4335 : for (auto i : index_range(names))
750 : {
751 192 : const std::string & name = names[i];
752 328 : if (const auto it = name_indices.find(name);
753 96 : it != name_indices.end())
754 : {
755 0 : returnval[i] = it->second;
756 0 : _node_integer_default_values[it->second] =
757 0 : default_values ? (*default_values)[i] : DofObject::invalid_id;
758 : }
759 : else
760 : {
761 424 : returnval[i] = _node_integer_names.size();
762 328 : name_indices[name] = returnval[i];
763 328 : _node_integer_names.push_back(name);
764 : _node_integer_default_values.push_back
765 396 : (default_values ? (*default_values)[i] : DofObject::invalid_id);
766 96 : added_an_integer = true;
767 : }
768 : }
769 :
770 4007 : if (allocate_data && added_an_integer)
771 114 : this->size_node_extra_integers();
772 :
773 5211 : return returnval;
774 : }
775 :
776 :
777 :
778 70 : unsigned int MeshBase::get_node_integer_index(std::string_view name) const
779 : {
780 196 : for (auto i : index_range(_node_integer_names))
781 176 : if (_node_integer_names[i] == name)
782 70 : return i;
783 :
784 0 : libmesh_error_msg("Unknown node integer " << name);
785 : return libMesh::invalid_uint;
786 : }
787 :
788 :
789 :
790 112 : bool MeshBase::has_node_integer(std::string_view name) const
791 : {
792 318 : for (auto & entry : _node_integer_names)
793 318 : if (entry == name)
794 32 : return true;
795 :
796 0 : return false;
797 : }
798 :
799 :
800 :
801 4345 : void MeshBase::remove_orphaned_nodes ()
802 : {
803 2584 : LOG_SCOPE("remove_orphaned_nodes()", "MeshBase");
804 :
805 : // Will hold the set of nodes that are currently connected to elements
806 1292 : std::unordered_set<Node *> connected_nodes;
807 :
808 : // Loop over the elements. Find which nodes are connected to at
809 : // least one of them.
810 657636 : for (const auto & element : this->element_ptr_range())
811 4359396 : for (auto & n : element->node_ref_range())
812 3762756 : connected_nodes.insert(&n);
813 :
814 1745084 : for (const auto & node : this->node_ptr_range())
815 723942 : if (!connected_nodes.count(node))
816 9733 : this->delete_node(node);
817 :
818 4345 : _preparation.has_removed_orphaned_nodes = true;
819 4345 : }
820 :
821 :
822 :
823 : #ifdef LIBMESH_ENABLE_DEPRECATED
824 0 : void MeshBase::prepare_for_use (const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
825 : {
826 : libmesh_deprecated();
827 :
828 : // We only respect the users wish if they tell us to skip renumbering. If they tell us not to
829 : // skip renumbering but someone previously called allow_renumbering(false), then the latter takes
830 : // precedence
831 0 : if (skip_renumber_nodes_and_elements)
832 0 : this->allow_renumbering(false);
833 :
834 : // We always accept the user's value for skip_find_neighbors, in contrast to skip_renumber
835 0 : const bool old_allow_find_neighbors = this->allow_find_neighbors();
836 0 : this->allow_find_neighbors(!skip_find_neighbors);
837 :
838 0 : this->prepare_for_use();
839 :
840 0 : this->allow_find_neighbors(old_allow_find_neighbors);
841 0 : }
842 :
843 0 : void MeshBase::prepare_for_use (const bool skip_renumber_nodes_and_elements)
844 : {
845 : libmesh_deprecated();
846 :
847 : // We only respect the users wish if they tell us to skip renumbering. If they tell us not to
848 : // skip renumbering but someone previously called allow_renumbering(false), then the latter takes
849 : // precedence
850 0 : if (skip_renumber_nodes_and_elements)
851 0 : this->allow_renumbering(false);
852 :
853 0 : this->prepare_for_use();
854 0 : }
855 : #endif // LIBMESH_ENABLE_DEPRECATED
856 :
857 :
858 :
859 :
860 44391 : void MeshBase::prepare_for_use ()
861 : {
862 : // Mark everything as unprepared, except for those things we've been
863 : // told we don't need to prepare, for backwards compatibility
864 44391 : this->clear_point_locator();
865 44391 : this->clear_stored_ranges();
866 44391 : _preparation = false;
867 44391 : _preparation.has_neighbor_ptrs = _skip_find_neighbors;
868 44391 : _preparation.has_removed_remote_elements = !_allow_remote_element_removal;
869 :
870 44391 : this->complete_preparation();
871 44391 : }
872 :
873 :
874 44692 : void MeshBase::complete_preparation()
875 : {
876 26228 : LOG_SCOPE("complete_preparation()", "MeshBase");
877 :
878 13114 : parallel_object_only();
879 :
880 13114 : libmesh_assert(this->comm().verify(this->is_serial()));
881 :
882 : // If we don't go into this method with valid constraint rows, we're
883 : // only going to be able to make that worse.
884 : #ifdef DEBUG
885 13114 : MeshTools::libmesh_assert_valid_constraint_rows(*this);
886 : #endif
887 :
888 : // A distributed mesh may have processors with no elements (or
889 : // processors with no elements of higher dimension, if we ever
890 : // support mixed-dimension meshes), but we want consistent
891 : // mesh_dimension anyways.
892 : //
893 : // cache_elem_data() should get the elem_dimensions() and
894 : // mesh_dimension() correct later, and we don't need it earlier.
895 :
896 :
897 : // Renumber the nodes and elements so that they in contiguous
898 : // blocks. By default, _skip_renumber_nodes_and_elements is false.
899 : //
900 : // Instances where you if prepare_for_use() should not renumber the nodes
901 : // and elements include reading in e.g. an xda/r or gmv file. In
902 : // this case, the ordering of the nodes may depend on an accompanying
903 : // solution, and the node ordering cannot be changed.
904 :
905 :
906 : // Mesh modification operations might not leave us with consistent
907 : // id counts, or might leave us with orphaned nodes we're no longer
908 : // using, but our partitioner might need that consistency and/or
909 : // might be confused by orphaned nodes.
910 44692 : if (!_skip_renumber_nodes_and_elements)
911 : {
912 40109 : if (!_preparation.has_removed_orphaned_nodes ||
913 42 : !_preparation.has_synched_id_counts)
914 40067 : this->renumber_nodes_and_elements();
915 : }
916 : else
917 : {
918 4583 : if (!_preparation.has_removed_orphaned_nodes)
919 4345 : this->remove_orphaned_nodes();
920 4583 : if (!_preparation.has_synched_id_counts)
921 4345 : this->update_parallel_id_counts();
922 : }
923 :
924 : // Let all the elements find their neighbors
925 44692 : if (!_skip_find_neighbors && !_preparation.has_neighbor_ptrs)
926 43209 : this->find_neighbors();
927 :
928 : // The user may have set boundary conditions. We require that the
929 : // boundary conditions were set consistently. Because we examine
930 : // neighbors when evaluating non-raw boundary condition IDs, this
931 : // assert is only valid when our neighbor links are in place.
932 : #ifdef DEBUG
933 13114 : MeshTools::libmesh_assert_valid_boundary_ids(*this);
934 : #endif
935 :
936 : // Search the mesh for all the dimensions of the elements
937 : // and cache them.
938 44692 : if (!_preparation.has_cached_elem_data)
939 44412 : this->cache_elem_data();
940 :
941 : // Search the mesh for elements that have a neighboring element
942 : // of dim+1 and set that element as the interior parent
943 44692 : if (!_preparation.has_interior_parent_ptrs)
944 : {
945 44391 : if (!_skip_detect_interior_parents)
946 44391 : this->detect_interior_parents();
947 : else
948 : {
949 : // We must set the flag that says "interior parent pointers have been set up"
950 : // even though we skip detect_interior_parents().
951 0 : _preparation.has_interior_parent_ptrs = true;
952 : }
953 : }
954 :
955 : // Fix up node unique ids in case mesh generation code didn't take
956 : // exceptional care to do so.
957 : // MeshCommunication().make_node_unique_ids_parallel_consistent(*this);
958 :
959 : // We're going to still require that mesh generation code gets
960 : // element unique ids consistent.
961 : #if defined(DEBUG) && defined(LIBMESH_ENABLE_UNIQUE_ID)
962 13114 : MeshTools::libmesh_assert_valid_unique_ids(*this);
963 : #endif
964 :
965 : // Allow our GhostingFunctor objects to reinit if necessary.
966 : // Do this before partitioning and redistributing, and before
967 : // deleting remote elements.
968 44692 : if (!_preparation.has_reinit_ghosting_functors)
969 44391 : this->reinit_ghosting_functors();
970 :
971 : // Partition the mesh unless *all* partitioning is to be skipped.
972 : // If only noncritical partitioning is to be skipped, the
973 : // partition() call will still check for orphaned nodes.
974 44692 : if (!skip_partitioning() && !_preparation.is_partitioned)
975 12918 : this->partition();
976 1286 : else if (!this->n_unpartitioned_elem() &&
977 196 : !this->n_unpartitioned_nodes())
978 643 : _preparation.is_partitioned = true;
979 :
980 : // If we're using DistributedMesh, we'll probably want it
981 : // parallelized.
982 44692 : if (this->_allow_remote_element_removal &&
983 43019 : !_preparation.has_removed_remote_elements)
984 42676 : this->delete_remote_elements();
985 : else
986 2016 : _preparation.has_removed_remote_elements = true;
987 :
988 : // Much of our boundary info may have been for now-remote parts of the mesh,
989 : // in which case we don't want to keep local copies of data meant to be
990 : // local. On the other hand we may have deleted, or the user may have added in
991 : // a distributed fashion, boundary data that is meant to be global. So we
992 : // handle both of those scenarios here
993 44692 : if (!_preparation.has_boundary_id_sets)
994 43215 : this->get_boundary_info().regenerate_id_sets();
995 :
996 44692 : if (!_skip_renumber_nodes_and_elements)
997 40109 : this->renumber_nodes_and_elements();
998 :
999 : // The mesh is now prepared for use, with the possible exception of
1000 : // partitioning that was supposed to be skipped, and it should know
1001 : // it.
1002 : #ifndef NDEBUG
1003 13114 : Preparation completed_preparation = _preparation;
1004 13114 : if (skip_partitioning())
1005 184 : completed_preparation.is_partitioned = true;
1006 13114 : libmesh_assert(completed_preparation);
1007 : #endif
1008 :
1009 : #ifdef DEBUG
1010 13114 : MeshTools::libmesh_assert_valid_boundary_ids(*this);
1011 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1012 13114 : MeshTools::libmesh_assert_valid_unique_ids(*this);
1013 : #endif
1014 : #endif
1015 44692 : }
1016 :
1017 : void
1018 44391 : MeshBase::reinit_ghosting_functors()
1019 : {
1020 100376 : for (auto & gf : _ghosting_functors)
1021 : {
1022 16572 : libmesh_assert(gf);
1023 55985 : gf->mesh_reinit();
1024 : }
1025 :
1026 44391 : _preparation.has_reinit_ghosting_functors = true;
1027 44391 : }
1028 :
1029 103549 : void MeshBase::clear ()
1030 : {
1031 : // Reset the number of partitions
1032 103549 : _n_parts = 1;
1033 :
1034 : // Reset the preparation flags
1035 103549 : _preparation = false;
1036 :
1037 : // Clear boundary information
1038 103549 : if (boundary_info)
1039 103423 : boundary_info->clear();
1040 :
1041 : // Clear cached element data
1042 30389 : _elem_dims.clear();
1043 30389 : _elem_default_orders.clear();
1044 103549 : _supported_nodal_order = MAXIMUM;
1045 :
1046 30389 : _elemset_codes.clear();
1047 30389 : _elemset_codes_inverse_map.clear();
1048 :
1049 30389 : _constraint_rows.clear();
1050 :
1051 : // Clear our point locator.
1052 103549 : this->clear_point_locator();
1053 103549 : this->clear_stored_ranges();
1054 103549 : }
1055 :
1056 :
1057 33241 : bool MeshBase::is_prepared() const
1058 : {
1059 33241 : return static_cast<bool>(_preparation);
1060 : }
1061 :
1062 :
1063 203 : void MeshBase::unset_is_prepared()
1064 : {
1065 203 : _preparation = false;
1066 203 : this->clear_point_locator();
1067 203 : this->clear_stored_ranges();
1068 203 : }
1069 :
1070 :
1071 99549 : void MeshBase::add_ghosting_functor(GhostingFunctor & ghosting_functor)
1072 : {
1073 : // We used to implicitly support duplicate inserts to std::set
1074 : #ifdef LIBMESH_ENABLE_DEPRECATED
1075 : _ghosting_functors.erase
1076 42193 : (std::remove(_ghosting_functors.begin(),
1077 : _ghosting_functors.end(),
1078 99549 : &ghosting_functor),
1079 86034 : _ghosting_functors.end());
1080 : #endif
1081 :
1082 : // We shouldn't have two copies of the same functor
1083 28678 : libmesh_assert(std::find(_ghosting_functors.begin(),
1084 : _ghosting_functors.end(),
1085 : &ghosting_functor) ==
1086 : _ghosting_functors.end());
1087 :
1088 99549 : _ghosting_functors.push_back(&ghosting_functor);
1089 99549 : }
1090 :
1091 :
1092 :
1093 96271 : void MeshBase::remove_ghosting_functor(GhostingFunctor & ghosting_functor)
1094 : {
1095 40939 : auto raw_it = std::find(_ghosting_functors.begin(),
1096 123937 : _ghosting_functors.end(), &ghosting_functor);
1097 :
1098 : // The DofMap has a "to_mesh" parameter that tells it to avoid
1099 : // registering a new functor with the mesh, but it doesn't keep
1100 : // track of which functors weren't added, so we'll support "remove a
1101 : // functor that isn't there" just like we did with set::erase
1102 : // before.
1103 96271 : if (raw_it != _ghosting_functors.end())
1104 96180 : _ghosting_functors.erase(raw_it);
1105 :
1106 : // We shouldn't have had two copies of the same functor
1107 27666 : libmesh_assert(std::find(_ghosting_functors.begin(),
1108 : _ghosting_functors.end(),
1109 : &ghosting_functor) ==
1110 : _ghosting_functors.end());
1111 :
1112 96271 : if (const auto it = _shared_functors.find(&ghosting_functor);
1113 27666 : it != _shared_functors.end())
1114 5 : _shared_functors.erase(it);
1115 96271 : }
1116 :
1117 :
1118 :
1119 3672 : void MeshBase::subdomain_ids (std::set<subdomain_id_type> & ids, const bool global /* = true */) const
1120 : {
1121 : // This requires an inspection on every processor
1122 1298 : if (global)
1123 1298 : parallel_object_only();
1124 :
1125 34 : struct SBDInserter {
1126 : std::set<subdomain_id_type> my_ids;
1127 :
1128 1298 : SBDInserter () {}
1129 36 : SBDInserter (SBDInserter &, Threads::split) {}
1130 :
1131 3778 : void operator()(const ConstElemRange & range) {
1132 696331 : for (const Elem * elem : range)
1133 692553 : my_ids.insert(elem->subdomain_id());
1134 3778 : }
1135 :
1136 36 : void join(SBDInserter & other) {
1137 36 : my_ids.merge(other.my_ids);
1138 36 : }
1139 : };
1140 :
1141 2596 : SBDInserter inserter;
1142 3672 : Threads::parallel_reduce(this->active_local_element_stored_range(), inserter);
1143 :
1144 1298 : ids.swap(inserter.my_ids);
1145 :
1146 3672 : if (global)
1147 : {
1148 : // Only include the unpartitioned elements if the user requests the global IDs.
1149 : // In the case of the local subdomain IDs, it doesn't make sense to include the
1150 : // unpartitioned elements because said elements do not have a sense of locality.
1151 9520 : for (const auto & elem : this->active_unpartitioned_element_ptr_range())
1152 4550 : ids.insert(elem->subdomain_id());
1153 :
1154 : // Some subdomains may only live on other processors
1155 3672 : this->comm().set_union(ids);
1156 : }
1157 3672 : }
1158 :
1159 :
1160 :
1161 42268 : void MeshBase::redistribute()
1162 : {
1163 : // We now have all elements and nodes redistributed; our ghosting
1164 : // functors should be ready to redistribute and/or recompute any
1165 : // cached data they use too.
1166 45452 : for (auto & gf : as_range(this->ghosting_functors_begin(),
1167 122811 : this->ghosting_functors_end()))
1168 51675 : gf->redistribute();
1169 42268 : }
1170 :
1171 :
1172 :
1173 29125 : void MeshBase::update_post_partitioning()
1174 : {
1175 : // A range over all elements might still be fine here, but any range
1176 : // over local elements is obsolete if our partitioner changed the
1177 : // definition of "local".
1178 9768 : _const_active_local_element_stored_range.reset(nullptr);
1179 29125 : }
1180 :
1181 :
1182 :
1183 2031 : subdomain_id_type MeshBase::n_subdomains() const
1184 : {
1185 : // This requires an inspection on every processor
1186 596 : parallel_object_only();
1187 :
1188 1192 : std::set<subdomain_id_type> ids;
1189 :
1190 2031 : this->subdomain_ids (ids);
1191 :
1192 2627 : return cast_int<subdomain_id_type>(ids.size());
1193 : }
1194 :
1195 :
1196 :
1197 0 : subdomain_id_type MeshBase::n_local_subdomains() const
1198 : {
1199 0 : std::set<subdomain_id_type> ids;
1200 :
1201 0 : this->subdomain_ids (ids, /* global = */ false);
1202 :
1203 0 : return cast_int<subdomain_id_type>(ids.size());
1204 : }
1205 :
1206 :
1207 :
1208 :
1209 15285 : dof_id_type MeshBase::n_nodes_on_proc (const processor_id_type proc_id) const
1210 : {
1211 : // We're either counting a processor's nodes or unpartitioned
1212 : // nodes
1213 5752 : libmesh_assert (proc_id < this->n_processors() ||
1214 : proc_id == DofObject::invalid_processor_id);
1215 :
1216 24818 : return static_cast<dof_id_type>(std::distance (this->pid_nodes_begin(proc_id),
1217 40103 : this->pid_nodes_end (proc_id)));
1218 : }
1219 :
1220 :
1221 :
1222 59377 : dof_id_type MeshBase::n_elem_on_proc (const processor_id_type proc_id) const
1223 : {
1224 : // We're either counting a processor's elements or unpartitioned
1225 : // elements
1226 18680 : libmesh_assert (proc_id < this->n_processors() ||
1227 : proc_id == DofObject::invalid_processor_id);
1228 :
1229 100074 : return static_cast<dof_id_type>(std::distance (this->pid_elements_begin(proc_id),
1230 159451 : this->pid_elements_end (proc_id)));
1231 : }
1232 :
1233 :
1234 :
1235 5921 : dof_id_type MeshBase::n_active_elem_on_proc (const processor_id_type proc_id) const
1236 : {
1237 2164 : libmesh_assert_less (proc_id, this->n_processors());
1238 9678 : return static_cast<dof_id_type>(std::distance (this->active_pid_elements_begin(proc_id),
1239 15599 : this->active_pid_elements_end (proc_id)));
1240 : }
1241 :
1242 :
1243 :
1244 0 : dof_id_type MeshBase::n_sub_elem () const
1245 : {
1246 0 : dof_id_type ne=0;
1247 :
1248 0 : for (const auto & elem : this->element_ptr_range())
1249 0 : ne += elem->n_sub_elem();
1250 :
1251 0 : return ne;
1252 : }
1253 :
1254 :
1255 :
1256 4210 : dof_id_type MeshBase::n_active_sub_elem () const
1257 : {
1258 1575 : dof_id_type ne=0;
1259 :
1260 2296668 : for (const auto & elem : this->active_element_ptr_range())
1261 2290883 : ne += elem->n_sub_elem();
1262 :
1263 4210 : return ne;
1264 : }
1265 :
1266 :
1267 :
1268 1424 : std::string MeshBase::get_info(const unsigned int verbosity /* = 0 */, const bool global /* = true */) const
1269 : {
1270 2264 : std::ostringstream oss;
1271 :
1272 1424 : oss << " Mesh Information:" << '\n';
1273 :
1274 1424 : if (!_elem_dims.empty())
1275 : {
1276 1424 : oss << " elem_dimensions()={";
1277 1424 : std::copy(_elem_dims.begin(),
1278 420 : --_elem_dims.end(), // --end() is valid if the set is non-empty
1279 1004 : std::ostream_iterator<unsigned int>(oss, ", "));
1280 1424 : oss << cast_int<unsigned int>(*_elem_dims.rbegin());
1281 1424 : oss << "}\n";
1282 : }
1283 :
1284 1424 : if (!_elem_default_orders.empty())
1285 : {
1286 1424 : oss << " elem_default_orders()={";
1287 1424 : std::transform(_elem_default_orders.begin(),
1288 420 : --_elem_default_orders.end(),
1289 1424 : std::ostream_iterator<std::string>(oss, ", "),
1290 8 : [](Order o)
1291 20 : { return Utility::enum_to_string<Order>(o); });
1292 1424 : oss << Utility::enum_to_string<Order>(*_elem_default_orders.rbegin());
1293 1424 : oss << "}\n";
1294 : }
1295 :
1296 1424 : oss << " supported_nodal_order()=" << this->supported_nodal_order() << '\n'
1297 2264 : << " spatial_dimension()=" << this->spatial_dimension() << '\n'
1298 1844 : << " n_nodes()=" << this->n_nodes() << '\n'
1299 1424 : << " n_local_nodes()=" << this->n_local_nodes() << '\n'
1300 1844 : << " n_elem()=" << this->n_elem() << '\n'
1301 2428 : << " n_local_elem()=" << this->n_local_elem() << '\n';
1302 : #ifdef LIBMESH_ENABLE_AMR
1303 2008 : oss << " n_active_elem()=" << this->n_active_elem() << '\n';
1304 : #endif
1305 1424 : if (global)
1306 1844 : oss << " n_subdomains()=" << static_cast<std::size_t>(this->n_subdomains()) << '\n';
1307 : else
1308 0 : oss << " n_local_subdomains()= " << static_cast<std::size_t>(this->n_local_subdomains()) << '\n';
1309 1844 : oss << " n_elemsets()=" << static_cast<std::size_t>(this->n_elemsets()) << '\n';
1310 1424 : if (!_elemset_codes.empty())
1311 0 : oss << " n_elemset_codes=" << _elemset_codes.size() << '\n';
1312 1424 : oss << " n_partitions()=" << static_cast<std::size_t>(this->n_partitions()) << '\n'
1313 2264 : << " n_processors()=" << static_cast<std::size_t>(this->n_processors()) << '\n'
1314 1844 : << " n_threads()=" << static_cast<std::size_t>(libMesh::n_threads()) << '\n'
1315 2264 : << " processor_id()=" << static_cast<std::size_t>(this->processor_id()) << '\n'
1316 1844 : << " is_prepared()=" << (this->is_prepared() ? "true" : "false") << '\n'
1317 3272 : << " is_replicated()=" << (this->is_replicated() ? "true" : "false") << '\n';
1318 :
1319 1424 : if (verbosity > 0)
1320 : {
1321 28 : if (global)
1322 : {
1323 8 : libmesh_parallel_only(this->comm());
1324 36 : if (this->processor_id() != 0)
1325 12 : oss << "\n Detailed global get_info() (verbosity > 0) is reduced and output to only rank 0.";
1326 : }
1327 :
1328 : // Helper for printing element types
1329 224 : const auto elem_type_helper = [](const std::set<int> & elem_types) {
1330 336 : std::stringstream ss;
1331 448 : for (auto it = elem_types.begin(); it != elem_types.end();)
1332 : {
1333 392 : ss << Utility::enum_to_string((ElemType)*it);
1334 224 : if (++it != elem_types.end())
1335 0 : ss << ", ";
1336 : }
1337 280 : return ss.str();
1338 112 : };
1339 :
1340 : // Helper for whether or not the given DofObject is to be included. If we're doing
1341 : // a global reduction, we also count unpartitioned objects on rank 0.
1342 248680 : const auto include_object = [this, &global](const DofObject & dof_object) {
1343 355712 : return this->processor_id() == dof_object.processor_id() ||
1344 48840 : (global &&
1345 31744 : this->processor_id() == 0 &&
1346 223852 : dof_object.processor_id() == DofObject::invalid_processor_id);
1347 20 : };
1348 :
1349 8 : Real volume = 0;
1350 :
1351 : // Add bounding box information
1352 28 : const auto bbox = global ? MeshTools::create_bounding_box(*this) : MeshTools::create_local_bounding_box(*this);
1353 28 : if (!global || this->processor_id() == 0)
1354 : oss << "\n " << (global ? "" : "Local ") << "Mesh Bounding Box:\n"
1355 16 : << " Minimum: " << bbox.min() << "\n"
1356 12 : << " Maximum: " << bbox.max() << "\n"
1357 24 : << " Delta: " << (bbox.max() - bbox.min()) << "\n";
1358 :
1359 : // Obtain the global or local element types
1360 16 : std::set<int> elem_types;
1361 16432 : for (const Elem * elem : this->active_local_element_ptr_range())
1362 16396 : elem_types.insert(elem->type());
1363 28 : if (global)
1364 : {
1365 : // Pick up unpartitioned elems on rank 0
1366 36 : if (this->processor_id() == 0)
1367 28 : for (const Elem * elem : this->active_unpartitioned_element_ptr_range())
1368 8 : elem_types.insert(elem->type());
1369 :
1370 28 : this->comm().set_union(elem_types);
1371 : }
1372 :
1373 : // Add element types
1374 28 : if (!global || this->processor_id() == 0)
1375 : oss << "\n " << (global ? "" : "Local ") << "Mesh Element Type(s):\n "
1376 44 : << elem_type_helper(elem_types) << "\n";
1377 :
1378 : // Reduce the nodeset ids
1379 16 : auto nodeset_ids = this->get_boundary_info().get_node_boundary_ids();
1380 28 : if (global)
1381 28 : this->comm().set_union(nodeset_ids);
1382 :
1383 : // Accumulate local information for each nodeset
1384 72 : struct NodesetInfo
1385 : {
1386 : std::size_t num_nodes = 0;
1387 : BoundingBox bbox;
1388 : };
1389 16 : std::map<boundary_id_type, NodesetInfo> nodeset_info_map;
1390 28420 : for (const auto & [node, id] : this->get_boundary_info().get_nodeset_map())
1391 : {
1392 28392 : if (!include_object(*node))
1393 8112 : continue;
1394 :
1395 16224 : NodesetInfo & info = nodeset_info_map[id];
1396 :
1397 16224 : ++info.num_nodes;
1398 :
1399 16224 : if (verbosity > 1)
1400 16224 : info.bbox.union_with(*node);
1401 : }
1402 :
1403 : // Add nodeset info
1404 28 : if (!global || this->processor_id() == 0)
1405 : {
1406 16 : oss << "\n " << (global ? "" : "Local ") << "Mesh Nodesets:\n";
1407 16 : if (nodeset_ids.empty())
1408 0 : oss << " None\n";
1409 : }
1410 :
1411 8 : const auto & nodeset_name_map = this->get_boundary_info().get_nodeset_name_map();
1412 196 : for (const auto id : nodeset_ids)
1413 : {
1414 168 : NodesetInfo & info = nodeset_info_map[id];
1415 :
1416 : // Reduce the local information for this nodeset if required
1417 168 : if (global)
1418 : {
1419 168 : this->comm().sum(info.num_nodes);
1420 168 : if (verbosity > 1)
1421 : {
1422 216 : this->comm().min(info.bbox.min());
1423 216 : this->comm().max(info.bbox.max());
1424 : }
1425 : }
1426 :
1427 168 : const bool has_name = nodeset_name_map.count(id) && nodeset_name_map.at(id).size();
1428 264 : const std::string name = has_name ? nodeset_name_map.at(id) : "";
1429 96 : if (global)
1430 48 : libmesh_assert(this->comm().verify(name));
1431 :
1432 168 : if (global ? this->processor_id() == 0 : info.num_nodes > 0)
1433 : {
1434 96 : oss << " Nodeset " << id;
1435 96 : if (has_name)
1436 144 : oss << " (" << name << ")";
1437 120 : oss << ", " << info.num_nodes << " " << (global ? "" : "local ") << "nodes\n";
1438 :
1439 96 : if (verbosity > 1)
1440 : {
1441 96 : oss << " " << (global ? "Bounding" : "Local bounding") << " box minimum: "
1442 96 : << info.bbox.min() << "\n"
1443 96 : << " " << (global ? "Bounding" : "Local bounding") << " box maximum: "
1444 96 : << info.bbox.max() << "\n"
1445 96 : << " " << (global ? "Bounding" : "Local bounding") << " box delta: "
1446 96 : << (info.bbox.max() - info.bbox.min()) << "\n";
1447 : }
1448 : }
1449 : }
1450 :
1451 : // Reduce the sideset ids
1452 16 : auto sideset_ids = this->get_boundary_info().get_side_boundary_ids();
1453 28 : if (global)
1454 28 : this->comm().set_union(sideset_ids);
1455 :
1456 : // Accumulate local information for each sideset
1457 : struct SidesetInfo
1458 : {
1459 : std::size_t num_sides = 0;
1460 : Real volume = 0;
1461 : std::set<int> side_elem_types;
1462 : std::set<int> elem_types;
1463 : std::set<dof_id_type> elem_ids;
1464 : std::set<dof_id_type> node_ids;
1465 : BoundingBox bbox;
1466 : };
1467 16 : ElemSideBuilder side_builder;
1468 16 : std::map<boundary_id_type, SidesetInfo> sideset_info_map;
1469 23324 : for (const auto & pair : this->get_boundary_info().get_sideset_map())
1470 : {
1471 23296 : const Elem * elem = pair.first;
1472 16640 : if (!include_object(*elem))
1473 9984 : continue;
1474 :
1475 13312 : const auto id = pair.second.second;
1476 13312 : SidesetInfo & info = sideset_info_map[id];
1477 :
1478 13312 : const auto s = pair.second.first;
1479 13312 : const Elem & side = side_builder(*elem, s);
1480 :
1481 13312 : ++info.num_sides;
1482 13312 : info.side_elem_types.insert(side.type());
1483 13312 : info.elem_types.insert(elem->type());
1484 13312 : info.elem_ids.insert(elem->id());
1485 :
1486 66560 : for (const Node & node : side.node_ref_range())
1487 39936 : if (include_object(node))
1488 52672 : info.node_ids.insert(node.id());
1489 :
1490 13312 : if (verbosity > 1)
1491 : {
1492 13312 : info.volume += side.volume();
1493 13312 : info.bbox.union_with(side.loose_bounding_box());
1494 : }
1495 : }
1496 :
1497 : // Add sideset info
1498 28 : if (!global || this->processor_id() == 0)
1499 : {
1500 16 : oss << "\n " << (global ? "" : "Local ") << "Mesh Sidesets:\n";
1501 16 : if (sideset_ids.empty())
1502 0 : oss << " None\n";
1503 : }
1504 8 : const auto & sideset_name_map = this->get_boundary_info().get_sideset_name_map();
1505 196 : for (const auto id : sideset_ids)
1506 : {
1507 168 : SidesetInfo & info = sideset_info_map[id];
1508 :
1509 168 : auto num_elems = info.elem_ids.size();
1510 168 : auto num_nodes = info.node_ids.size();
1511 :
1512 : // Reduce the local information for this sideset if required
1513 168 : if (global)
1514 : {
1515 168 : this->comm().sum(info.num_sides);
1516 168 : this->comm().set_union(info.side_elem_types, 0);
1517 168 : this->comm().sum(num_elems);
1518 168 : this->comm().set_union(info.elem_types, 0);
1519 168 : this->comm().sum(num_nodes);
1520 168 : if (verbosity > 1)
1521 : {
1522 168 : this->comm().sum(info.volume);
1523 216 : this->comm().min(info.bbox.min());
1524 216 : this->comm().max(info.bbox.max());
1525 : }
1526 : }
1527 :
1528 168 : const bool has_name = sideset_name_map.count(id) && sideset_name_map.at(id).size();
1529 264 : const std::string name = has_name ? sideset_name_map.at(id) : "";
1530 96 : if (global)
1531 48 : libmesh_assert(this->comm().verify(name));
1532 :
1533 168 : if (global ? this->processor_id() == 0 : info.num_sides > 0)
1534 : {
1535 96 : oss << " Sideset " << id;
1536 96 : if (has_name)
1537 144 : oss << " (" << name << ")";
1538 168 : oss << ", " << info.num_sides << " sides (" << elem_type_helper(info.side_elem_types) << ")"
1539 288 : << ", " << num_elems << " " << (global ? "" : "local ") << "elems (" << elem_type_helper(info.elem_types) << ")"
1540 216 : << ", " << num_nodes << " " << (global ? "" : "local ") << "nodes\n";
1541 :
1542 96 : if (verbosity > 1)
1543 : {
1544 120 : oss << " " << (global ? "Side" : "Local side") << " volume: " << info.volume << "\n"
1545 96 : << " " << (global ? "Bounding" : "Local bounding") << " box minimum: "
1546 96 : << info.bbox.min() << "\n"
1547 96 : << " " << (global ? "Bounding" : "Local bounding") << " box maximum: "
1548 96 : << info.bbox.max() << "\n"
1549 96 : << " " << (global ? "Bounding" : "Local bounding") << " box delta: "
1550 96 : << (info.bbox.max() - info.bbox.min()) << "\n";
1551 : }
1552 : }
1553 : }
1554 :
1555 : // Reduce the edgeset ids
1556 16 : auto edgeset_ids = this->get_boundary_info().get_edge_boundary_ids();
1557 28 : if (global)
1558 28 : this->comm().set_union(edgeset_ids);
1559 :
1560 : // Accumulate local information for each edgeset
1561 0 : struct EdgesetInfo
1562 : {
1563 : std::size_t num_edges = 0;
1564 : std::set<int> edge_elem_types;
1565 : BoundingBox bbox;
1566 : };
1567 16 : std::map<boundary_id_type, EdgesetInfo> edgeset_info_map;
1568 28 : std::unique_ptr<const Elem> edge;
1569 :
1570 28 : for (const auto & pair : this->get_boundary_info().get_edgeset_map())
1571 : {
1572 0 : const Elem * elem = pair.first;
1573 0 : if (!include_object(*elem))
1574 0 : continue;
1575 :
1576 0 : const auto id = pair.second.second;
1577 0 : EdgesetInfo & info = edgeset_info_map[id];
1578 :
1579 0 : elem->build_edge_ptr(edge, pair.second.first);
1580 :
1581 0 : ++info.num_edges;
1582 0 : info.edge_elem_types.insert(edge->type());
1583 :
1584 0 : if (verbosity > 1)
1585 0 : info.bbox.union_with(edge->loose_bounding_box());
1586 : }
1587 :
1588 : // Add edgeset info
1589 28 : if (!global || this->processor_id() == 0)
1590 : {
1591 16 : oss << "\n " << (global ? "" : "Local ") << "Mesh Edgesets:\n";
1592 16 : if (edgeset_ids.empty())
1593 16 : oss << " None\n";
1594 : }
1595 :
1596 8 : const auto & edgeset_name_map = this->get_boundary_info().get_edgeset_name_map();
1597 28 : for (const auto id : edgeset_ids)
1598 : {
1599 0 : EdgesetInfo & info = edgeset_info_map[id];
1600 :
1601 : // Reduce the local information for this edgeset if required
1602 0 : if (global)
1603 : {
1604 0 : this->comm().sum(info.num_edges);
1605 0 : this->comm().set_union(info.edge_elem_types, 0);
1606 0 : if (verbosity > 1)
1607 : {
1608 0 : this->comm().min(info.bbox.min());
1609 0 : this->comm().min(info.bbox.max());
1610 : }
1611 : }
1612 :
1613 0 : const bool has_name = edgeset_name_map.count(id) && edgeset_name_map.at(id).size();
1614 0 : const std::string name = has_name ? edgeset_name_map.at(id) : "";
1615 0 : if (global)
1616 0 : libmesh_assert(this->comm().verify(name));
1617 :
1618 0 : if (global ? this->processor_id() == 0 : info.num_edges > 0)
1619 : {
1620 0 : oss << " Edgeset " << id;
1621 0 : if (has_name)
1622 0 : oss << " (" << name << ")";
1623 0 : oss << ", " << info.num_edges << " " << (global ? "" : "local ") << "edges ("
1624 0 : << elem_type_helper(info.edge_elem_types) << ")\n";
1625 :
1626 0 : if (verbosity > 1)
1627 : {
1628 0 : oss << " " << (global ? "Bounding" : "Local bounding") << " box minimum: "
1629 0 : << info.bbox.min() << "\n"
1630 0 : << " " << (global ? "Bounding" : "Local bounding") << " box maximum: "
1631 0 : << info.bbox.max() << "\n"
1632 0 : << " " << (global ? "Bounding" : "Local bounding") << " box delta: "
1633 0 : << (info.bbox.max() - info.bbox.min()) << "\n";
1634 : }
1635 : }
1636 : }
1637 :
1638 : // Reduce the block IDs and block names
1639 16 : std::set<subdomain_id_type> subdomains;
1640 49200 : for (const Elem * elem : this->active_element_ptr_range())
1641 20480 : if (include_object(*elem))
1642 16396 : subdomains.insert(elem->subdomain_id());
1643 28 : if (global)
1644 28 : this->comm().set_union(subdomains);
1645 :
1646 : // Accumulate local information for each subdomain
1647 0 : struct SubdomainInfo
1648 : {
1649 : std::size_t num_elems = 0;
1650 : Real volume = 0;
1651 : std::set<int> elem_types;
1652 : std::set<dof_id_type> active_node_ids;
1653 : #ifdef LIBMESH_ENABLE_AMR
1654 : std::size_t num_active_elems = 0;
1655 : #endif
1656 : BoundingBox bbox;
1657 : };
1658 16 : std::map<subdomain_id_type, SubdomainInfo> subdomain_info_map;
1659 49200 : for (const Elem * elem : this->element_ptr_range())
1660 20480 : if (include_object(*elem))
1661 : {
1662 16384 : SubdomainInfo & info = subdomain_info_map[elem->subdomain_id()];
1663 :
1664 16384 : ++info.num_elems;
1665 16384 : info.elem_types.insert(elem->type());
1666 :
1667 : #ifdef LIBMESH_ENABLE_AMR
1668 4096 : if (elem->active())
1669 16384 : ++info.num_active_elems;
1670 : #endif
1671 :
1672 147456 : for (const Node & node : elem->node_ref_range())
1673 130560 : if (include_object(node) && node.active())
1674 129536 : info.active_node_ids.insert(node.id());
1675 :
1676 16384 : if (verbosity > 1 && elem->active())
1677 : {
1678 16384 : info.volume += elem->volume();
1679 16384 : info.bbox.union_with(elem->loose_bounding_box());
1680 : }
1681 12 : }
1682 :
1683 : // Add subdomain info
1684 28 : oss << "\n " << (global ? "" : "Local ") << "Mesh Subdomains:\n";
1685 8 : const auto & subdomain_name_map = this->get_subdomain_name_map();
1686 56 : for (const auto id : subdomains)
1687 : {
1688 28 : SubdomainInfo & info = subdomain_info_map[id];
1689 :
1690 28 : auto num_active_nodes = info.active_node_ids.size();
1691 :
1692 : // Reduce the information for this subdomain if needed
1693 28 : if (global)
1694 : {
1695 28 : this->comm().sum(info.num_elems);
1696 : #ifdef LIBMESH_ENABLE_AMR
1697 28 : this->comm().sum(info.num_active_elems);
1698 : #endif
1699 28 : this->comm().sum(num_active_nodes);
1700 28 : this->comm().set_union(info.elem_types, 0);
1701 28 : if (verbosity > 1)
1702 : {
1703 36 : this->comm().min(info.bbox.min());
1704 36 : this->comm().max(info.bbox.max());
1705 28 : this->comm().sum(info.volume);
1706 : }
1707 : }
1708 28 : if (verbosity > 1)
1709 28 : volume += info.volume;
1710 :
1711 8 : const bool has_name = subdomain_name_map.count(id);
1712 36 : const std::string name = has_name ? subdomain_name_map.at(id) : "";
1713 16 : if (global)
1714 8 : libmesh_assert(this->comm().verify(name));
1715 :
1716 28 : if (!global || this->processor_id() == 0)
1717 : {
1718 16 : oss << " Subdomain " << id;
1719 16 : if (has_name)
1720 0 : oss << " (" << name << ")";
1721 16 : oss << ": " << info.num_elems << " " << (global ? "" : "local ") << "elems "
1722 24 : << "(" << elem_type_helper(info.elem_types);
1723 : #ifdef LIBMESH_ENABLE_AMR
1724 16 : oss << ", " << info.num_active_elems << " active";
1725 : #endif
1726 20 : oss << "), " << num_active_nodes << " " << (global ? "" : "local ") << "active nodes\n";
1727 16 : if (verbosity > 1)
1728 : {
1729 20 : oss << " " << (global ? "Volume" : "Local volume") << ": " << info.volume << "\n";
1730 16 : oss << " " << (global ? "Bounding" : "Local bounding") << " box minimum: "
1731 16 : << info.bbox.min() << "\n"
1732 16 : << " " << (global ? "Bounding" : "Local bounding") << " box maximum: "
1733 16 : << info.bbox.max() << "\n"
1734 16 : << " " << (global ? "Bounding" : "Local bounding") << " box delta: "
1735 16 : << (info.bbox.max() - info.bbox.min()) << "\n";
1736 : }
1737 : }
1738 : }
1739 :
1740 48 : oss << " " << (global ? "Global" : "Local") << " mesh volume = " << volume << "\n";
1741 :
1742 12 : }
1743 :
1744 1844 : return oss.str();
1745 584 : }
1746 :
1747 :
1748 1424 : void MeshBase::print_info(std::ostream & os, const unsigned int verbosity /* = 0 */, const bool global /* = true */) const
1749 : {
1750 1844 : os << this->get_info(verbosity, global)
1751 420 : << std::endl;
1752 1424 : }
1753 :
1754 :
1755 0 : std::ostream & operator << (std::ostream & os, const MeshBase & m)
1756 : {
1757 0 : m.print_info(os);
1758 0 : return os;
1759 : }
1760 :
1761 :
1762 44133 : void MeshBase::partition (const unsigned int n_parts)
1763 : {
1764 : // If we get here and we have unpartitioned elements, we need that
1765 : // fixed.
1766 44133 : if (this->n_unpartitioned_elem() > 0)
1767 : {
1768 9412 : libmesh_assert (partitioner().get());
1769 9412 : libmesh_assert (this->is_serial());
1770 32206 : partitioner()->partition (*this, n_parts);
1771 : }
1772 : // A nullptr partitioner or a skip_partitioning(true) call or a
1773 : // skip_noncritical_partitioning(true) call means don't repartition;
1774 : // skip_noncritical_partitioning() checks all these.
1775 3530 : else if (!skip_noncritical_partitioning())
1776 : {
1777 11342 : partitioner()->partition (*this, n_parts);
1778 : }
1779 : else
1780 : {
1781 : // Adaptive coarsening may have "orphaned" nodes on processors
1782 : // whose elements no longer share them. We need to check for
1783 : // and possibly fix that.
1784 585 : MeshTools::correct_node_proc_ids(*this);
1785 :
1786 : // Make sure locally cached partition count is correct
1787 585 : this->recalculate_n_partitions();
1788 :
1789 : // Make sure any other locally cached data is correct
1790 585 : this->update_post_partitioning();
1791 : }
1792 :
1793 44133 : _preparation.is_partitioned = true;
1794 44133 : }
1795 :
1796 1888 : void MeshBase::all_second_order (const bool full_ordered)
1797 : {
1798 1888 : this->all_second_order_range(this->element_ptr_range(), full_ordered);
1799 1888 : }
1800 :
1801 1474 : void MeshBase::all_complete_order ()
1802 : {
1803 1474 : this->all_complete_order_range(this->element_ptr_range());
1804 1474 : }
1805 :
1806 585 : unsigned int MeshBase::recalculate_n_partitions()
1807 : {
1808 : // This requires an inspection on every processor
1809 170 : parallel_object_only();
1810 :
1811 585 : unsigned int max_proc_id=0;
1812 :
1813 69322 : for (const auto & elem : this->active_local_element_ptr_range())
1814 68725 : max_proc_id = std::max(max_proc_id, static_cast<unsigned int>(elem->processor_id()));
1815 :
1816 : // The number of partitions is one more than the max processor ID.
1817 585 : _n_parts = max_proc_id+1;
1818 :
1819 585 : this->comm().max(_n_parts);
1820 :
1821 585 : return _n_parts;
1822 : }
1823 :
1824 :
1825 :
1826 179757 : std::unique_ptr<PointLocatorBase> MeshBase::sub_point_locator () const
1827 : {
1828 : // If there's no master point locator, then we need one.
1829 179757 : if (_point_locator.get() == nullptr)
1830 : {
1831 : // PointLocator construction may not be safe within threads
1832 760 : libmesh_assert(!Threads::in_threads);
1833 :
1834 : // And it may require parallel communication
1835 760 : parallel_object_only();
1836 :
1837 : #ifdef LIBMESH_ENABLE_NANOFLANN_POINTLOCATOR
1838 : _point_locator = PointLocatorBase::build(NANOFLANN, *this);
1839 : #else
1840 3360 : _point_locator = PointLocatorBase::build(TREE_ELEMENTS, *this);
1841 : #endif
1842 :
1843 2440 : if (_point_locator_close_to_point_tol > 0.)
1844 0 : _point_locator->set_close_to_point_tol(_point_locator_close_to_point_tol);
1845 : }
1846 :
1847 : // Otherwise there was a master point locator, and we can grab a
1848 : // sub-locator easily.
1849 : return
1850 : #ifdef LIBMESH_ENABLE_NANOFLANN_POINTLOCATOR
1851 : PointLocatorBase::build(NANOFLANN, *this, _point_locator.get());
1852 : #else
1853 179757 : PointLocatorBase::build(TREE_ELEMENTS, *this, _point_locator.get());
1854 : #endif
1855 : }
1856 :
1857 :
1858 :
1859 5164238 : void MeshBase::clear_point_locator ()
1860 : {
1861 1564851 : _point_locator.reset(nullptr);
1862 5164238 : }
1863 :
1864 :
1865 :
1866 0 : void MeshBase::set_count_lower_dim_elems_in_point_locator(bool count_lower_dim_elems)
1867 : {
1868 0 : _count_lower_dim_elems_in_point_locator = count_lower_dim_elems;
1869 0 : }
1870 :
1871 :
1872 :
1873 1564877 : bool MeshBase::get_count_lower_dim_elems_in_point_locator() const
1874 : {
1875 1564877 : return _count_lower_dim_elems_in_point_locator;
1876 : }
1877 :
1878 :
1879 :
1880 159173 : std::string & MeshBase::subdomain_name(subdomain_id_type id)
1881 : {
1882 159173 : return _block_id_to_name[id];
1883 : }
1884 :
1885 5487 : const std::string & MeshBase::subdomain_name(subdomain_id_type id) const
1886 : {
1887 : // An empty string to return when no matching subdomain name is found
1888 5487 : static const std::string empty;
1889 :
1890 5487 : if (const auto iter = _block_id_to_name.find(id);
1891 1623 : iter == _block_id_to_name.end())
1892 1417 : return empty;
1893 : else
1894 828 : return iter->second;
1895 : }
1896 :
1897 :
1898 :
1899 :
1900 0 : subdomain_id_type MeshBase::get_id_by_name(std::string_view name) const
1901 : {
1902 : // Linear search over the map values.
1903 0 : for (const auto & [sbd_id, sbd_name] : _block_id_to_name)
1904 0 : if (sbd_name == name)
1905 0 : return sbd_id;
1906 :
1907 : // If we made it here without returning, we don't have a subdomain
1908 : // with the requested name, so return Elem::invalid_subdomain_id.
1909 0 : return Elem::invalid_subdomain_id;
1910 : }
1911 :
1912 :
1913 158081 : const ElemRange & MeshBase::element_stored_range()
1914 : {
1915 158081 : if (!_element_stored_range)
1916 : {
1917 : // Range construction may not be safe within threads
1918 13984 : libmesh_assert(!Threads::in_threads);
1919 :
1920 : _element_stored_range =
1921 113666 : std::make_unique<ElemRange>(this->elements_begin(),
1922 108396 : this->elements_end());
1923 : }
1924 :
1925 158081 : return *_element_stored_range;
1926 : }
1927 :
1928 30850 : const ConstElemRange & MeshBase::active_local_element_stored_range() const
1929 : {
1930 30850 : if (!_const_active_local_element_stored_range)
1931 : {
1932 : // Range construction may not be safe within threads
1933 1754 : libmesh_assert(!Threads::in_threads);
1934 :
1935 : _const_active_local_element_stored_range =
1936 11624 : std::make_unique<ConstElemRange>(this->active_local_elements_begin(),
1937 11842 : this->active_local_elements_end());
1938 : }
1939 :
1940 30850 : return *_const_active_local_element_stored_range;
1941 : }
1942 :
1943 5148280 : void MeshBase::clear_stored_ranges()
1944 : {
1945 1561331 : _element_stored_range.reset(nullptr);
1946 1561331 : _const_active_local_element_stored_range.reset(nullptr);
1947 5148280 : }
1948 :
1949 :
1950 : #ifdef LIBMESH_ENABLE_DEPRECATED
1951 0 : void MeshBase::cache_elem_dims()
1952 : {
1953 : libmesh_deprecated();
1954 :
1955 0 : this->cache_elem_data();
1956 0 : }
1957 : #endif // LIBMESH_ENABLE_DEPRECATED
1958 :
1959 51238 : void MeshBase::cache_elem_data()
1960 : {
1961 : // This requires an inspection on every processor
1962 15108 : parallel_object_only();
1963 :
1964 : // Need to clear containers first in case all elements of a
1965 : // particular dimension/order/subdomain have been deleted.
1966 30200 : _elem_dims.clear();
1967 30200 : _elem_default_orders.clear();
1968 30200 : _mesh_subdomains.clear();
1969 51238 : _supported_nodal_order = MAXIMUM;
1970 :
1971 4409615 : for (const auto & elem : this->active_element_ptr_range())
1972 : {
1973 4322247 : _elem_dims.insert(cast_int<unsigned char>(elem->dim()));
1974 4322247 : _elem_default_orders.insert(elem->default_order());
1975 4322247 : _mesh_subdomains.insert(elem->subdomain_id());
1976 4322247 : _supported_nodal_order =
1977 4322247 : static_cast<Order>
1978 8644494 : (std::min(static_cast<int>(_supported_nodal_order),
1979 4358375 : static_cast<int>(elem->supported_nodal_order())));
1980 21038 : }
1981 :
1982 51238 : if (!this->is_serial())
1983 : {
1984 : // Some different dimension/order/subdomain elements may only live
1985 : // on other processors
1986 540 : this->comm().set_union(_elem_dims);
1987 540 : this->comm().set_union(_elem_default_orders);
1988 540 : this->comm().min(_supported_nodal_order);
1989 540 : this->comm().set_union(_mesh_subdomains);
1990 : }
1991 :
1992 : // If the largest element dimension found is larger than the current
1993 : // _spatial_dimension, increase _spatial_dimension.
1994 51238 : unsigned int max_dim = this->mesh_dimension();
1995 51238 : if (max_dim > _spatial_dimension)
1996 4875 : _spatial_dimension = cast_int<unsigned char>(max_dim);
1997 :
1998 : // _spatial_dimension may need to increase from 1->2 or 2->3 if the
1999 : // mesh is full of 1D elements but they are not x-aligned, or the
2000 : // mesh is full of 2D elements but they are not in the x-y plane.
2001 : // If the mesh is x-aligned or x-y planar, we will end up checking
2002 : // every node's coordinates and not breaking out of the loop
2003 : // early...
2004 51238 : if (_spatial_dimension < LIBMESH_DIM)
2005 : {
2006 5133874 : for (const auto & node : this->node_ptr_range())
2007 : {
2008 : // Note: the exact floating point comparison is intentional,
2009 : // we don't want to get tripped up by tolerances.
2010 3654025 : if ((*node)(0) != 0. && _spatial_dimension < 1)
2011 0 : _spatial_dimension = 1;
2012 :
2013 3654025 : if ((*node)(1) != 0. && _spatial_dimension < 2)
2014 : {
2015 292 : _spatial_dimension = 2;
2016 : #if LIBMESH_DIM == 2
2017 : // If libmesh is compiled in 2D mode, this is the
2018 : // largest spatial dimension possible so we can break
2019 : // out.
2020 : break;
2021 : #endif
2022 : }
2023 :
2024 : #if LIBMESH_DIM > 2
2025 3654025 : if ((*node)(2) != 0.)
2026 : {
2027 : // Spatial dimension can't get any higher than this, so
2028 : // we can break out.
2029 973 : _spatial_dimension = 3;
2030 973 : break;
2031 : }
2032 : #endif
2033 11280 : }
2034 : }
2035 :
2036 51238 : _preparation.has_cached_elem_data = true;
2037 51238 : }
2038 :
2039 :
2040 5460 : void MeshBase::sync_subdomain_name_map()
2041 : {
2042 : // This requires every processor
2043 1628 : parallel_object_only();
2044 :
2045 5460 : this->comm().set_union(_block_id_to_name);
2046 5460 : }
2047 :
2048 :
2049 44391 : void MeshBase::detect_interior_parents()
2050 : {
2051 13028 : LOG_SCOPE("detect_interior_parents()", "MeshBase");
2052 :
2053 : // This requires an inspection on every processor
2054 13028 : parallel_object_only();
2055 :
2056 : // This requires up-to-date mesh dimensions in cache
2057 13028 : libmesh_assert(_preparation.has_cached_elem_data);
2058 :
2059 : // Early return if the mesh is empty or has elements of a single spatial dimension.
2060 44391 : if (this->elem_dimensions().size() <= 1)
2061 : {
2062 43722 : _preparation.has_interior_parent_ptrs = true;
2063 43787 : return;
2064 : }
2065 :
2066 : // Convenient elem_dimensions iterators
2067 481 : const auto dim_start = this->elem_dimensions().begin();
2068 188 : const auto dim_end = this->elem_dimensions().end();
2069 :
2070 : // In this function we find only +1 dimensional interior parents,
2071 : // (so, for a given element el, the interior parent p must satisfy p.dim() == el.dim() + 1).
2072 : // Therefore, we can avoid checking the existence of interior parents
2073 : // for all those elements el such there there is no p with p.dim() == el.dim() + 1.
2074 : // We store whether to skip any given dimension in the construction of interior parents
2075 : // inside the vector in dimensions_to_skip_for_interior_parents.
2076 669 : std::vector<bool> skip_dimension_for_interior_parents(/*count=*/LIBMESH_DIM+1, /*value=*/false);
2077 481 : skip_dimension_for_interior_parents.back() = true;
2078 :
2079 : // Moreover, in the following, we will build a node-to-elem map.
2080 : // It is among the elems of this map that we will look for interior parents.
2081 : // Therefore, we can skip all elems p such that there is no el with el.dim() == p.dim() - 1.
2082 : // We store whether to skip any given dimension in the construction of the node-to-elem map
2083 : // in the vector skip_dimensions_for_node_to_el_map.
2084 669 : std::vector<bool> skip_dimensions_for_node_to_el_map(/*count=*/LIBMESH_DIM+1, /*value=*/false);
2085 669 : skip_dimensions_for_node_to_el_map[*dim_start] = true;
2086 :
2087 : // We also create a flag to know if all dimensions should be skipped,
2088 : // and if we should therefore return early.
2089 188 : bool skip_all_dimensions = true;
2090 :
2091 : // Fill dimensions_to_skip_for_interior_parents and dimensions_to_skip_for_node_to_el_map.
2092 861 : for (auto [it, next] = std::make_tuple(dim_start, std::next(dim_start));
2093 1544 : next != dim_end; ++it, ++next)
2094 : {
2095 683 : if (*it + 1 != *next) // if sequential dimensions differ by exactly 1
2096 : {
2097 80 : skip_dimension_for_interior_parents[*it] = true;
2098 151 : skip_dimensions_for_node_to_el_map[*next] = true;
2099 : }
2100 684 : else if (!skip_dimension_for_interior_parents[*it])
2101 152 : skip_all_dimensions = false;
2102 : }
2103 :
2104 : // There is nothing to do if all dimensions should be
2105 : // skipped. Before returning, we must also set the flag that says
2106 : // "interior parent pointers have been set up" even though we
2107 : // determined there was no work to be done.
2108 669 : if (skip_all_dimensions)
2109 : {
2110 137 : _preparation.has_interior_parent_ptrs = true;
2111 36 : return;
2112 : }
2113 :
2114 : // Do we have interior parent pointers going to a different mesh?
2115 : // If so then we'll still check to make sure that's the only place
2116 : // they go, so we can libmesh_not_implemented() if not.
2117 304 : const bool separate_interior_mesh = (&(this->interior_mesh()) != this);
2118 :
2119 : // This map will be used to set interior parents
2120 152 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> node_to_elem;
2121 :
2122 41642 : for (const auto & elem : this->element_ptr_range())
2123 : {
2124 : // Ignore element if it cannot be interior parent of any other elem.
2125 28511 : if (skip_dimensions_for_node_to_el_map[elem->dim()])
2126 5755 : continue;
2127 :
2128 : // Populating the node_to_elem map, same as MeshTools::build_nodes_to_elem_map
2129 107926 : for (auto n : make_range(elem->n_vertices()))
2130 : {
2131 24992 : libmesh_assert_less (elem->id(), this->max_elem_id());
2132 :
2133 112464 : node_to_elem[elem->node_id(n)].push_back(elem->id());
2134 : }
2135 228 : }
2136 :
2137 : // Automatically set interior parents
2138 41642 : for (const auto & element : this->element_ptr_range())
2139 : {
2140 : // Ignore elements with dimensions to skip
2141 : // or elements that already have an interior parent.
2142 28511 : if (skip_dimension_for_interior_parents[element->dim()] || element->interior_parent())
2143 4641 : continue;
2144 :
2145 : // Start by generating sets of dim+1 dimensional elements that
2146 : // touch each vertex of the current element. If we encounter a
2147 : // vertex not connected to _any_ dim+1 dimensional elements,
2148 : // then we can exit the loop without checking the remaining
2149 : // vertices since an interior parent (if it exists) will be
2150 : // connected to all vertices of the current element.
2151 37510 : std::vector<std::set<dof_id_type>> neighbors( element->n_vertices() );
2152 :
2153 6820 : bool found_interior_parents = true;
2154 :
2155 24185 : for (auto n : make_range(element->n_vertices()))
2156 : {
2157 30969 : auto it = node_to_elem.find(element->node_id(n));
2158 :
2159 : // Check at first that this node is not isolated.
2160 24087 : if (it == node_to_elem.end())
2161 : {
2162 1472 : found_interior_parents = false;
2163 6792 : break; // out of n-loop
2164 : }
2165 :
2166 105399 : for (const auto & vertex_neighbor_id : it->second)
2167 86464 : if (this->elem_ref(vertex_neighbor_id).dim() == element->dim()+1)
2168 861 : neighbors[n].insert(vertex_neighbor_id);
2169 :
2170 24345 : if (neighbors[n].empty())
2171 : {
2172 : // We have found an empty set for one vertex, no reason
2173 : // to continue.
2174 5320 : found_interior_parents = false;
2175 5320 : break; // out of n-loop
2176 : }
2177 : }
2178 :
2179 : // If we have generated a non-empty set of elements for each
2180 : // vertex, we will now look for a vertex_neighbor_id that
2181 : // appears in _all_ of those sets. If found, this is our interior
2182 : // parent id. If multiple such common ids are found, we will
2183 : // take the lowest such id to be the interior parent id.
2184 6820 : if (found_interior_parents)
2185 : {
2186 56 : std::set<dof_id_type> & neighbors_0 = neighbors[0];
2187 161 : for (const auto & interior_parent_id : neighbors_0)
2188 : {
2189 40 : found_interior_parents = false;
2190 259 : for (auto n : make_range(1u, element->n_vertices()))
2191 : {
2192 161 : if (neighbors[n].count(interior_parent_id))
2193 : {
2194 34 : found_interior_parents = true;
2195 : }
2196 : else
2197 : {
2198 12 : found_interior_parents = false;
2199 12 : break;
2200 : }
2201 : }
2202 :
2203 140 : if (found_interior_parents)
2204 : {
2205 77 : element->set_interior_parent(this->elem_ptr(interior_parent_id));
2206 22 : break;
2207 : }
2208 : }
2209 :
2210 : // Do we have a mixed dimensional mesh that contains some of
2211 : // its own interior parents, but we already expect to have
2212 : // interior parents on a different mesh? That's going to
2213 : // take some work to support if anyone needs it.
2214 98 : if (separate_interior_mesh)
2215 0 : libmesh_not_implemented_msg
2216 : ("interior_parent() values in multiple meshes are unsupported.");
2217 : }
2218 10458 : }
2219 :
2220 : // This flag doesn't necessarily mean any Elems actually have
2221 : // interior parent pointers, just that we did all the work to
2222 : // determine whether or not they do.
2223 532 : _preparation.has_interior_parent_ptrs = true;
2224 : }
2225 :
2226 :
2227 :
2228 : #ifdef LIBMESH_ENABLE_PERIODIC
2229 : /**
2230 : * Register a pair of boundaries as disjoint neighbor boundary pairs.
2231 : */
2232 98 : void MeshBase::add_disjoint_neighbor_boundary_pairs(const boundary_id_type b1,
2233 : const boundary_id_type b2,
2234 : const RealVectorValue & translation)
2235 : {
2236 : // Lazily allocate the container the first time it’s needed
2237 98 : if (!_disjoint_neighbor_boundary_pairs)
2238 88 : _disjoint_neighbor_boundary_pairs = std::make_unique<PeriodicBoundaries>();
2239 :
2240 28 : PeriodicBoundaries & db = *_disjoint_neighbor_boundary_pairs;
2241 :
2242 : // Create forward and inverse boundary mappings
2243 126 : PeriodicBoundary forward(translation);
2244 98 : PeriodicBoundary inverse(translation * -1.0);
2245 :
2246 98 : forward.myboundary = b1;
2247 98 : forward.pairedboundary = b2;
2248 98 : inverse.myboundary = b2;
2249 98 : inverse.pairedboundary = b1;
2250 :
2251 : // Add both directions into the container
2252 98 : db.emplace(b1, forward.clone());
2253 126 : db.emplace(b2, inverse.clone());
2254 98 : }
2255 :
2256 47808 : PeriodicBoundaries * MeshBase::get_disjoint_neighbor_boundary_pairs()
2257 : {
2258 47808 : return _disjoint_neighbor_boundary_pairs.get();
2259 : }
2260 :
2261 94436 : const PeriodicBoundaries * MeshBase::get_disjoint_neighbor_boundary_pairs() const
2262 : {
2263 94436 : return _disjoint_neighbor_boundary_pairs.get();
2264 : }
2265 :
2266 21 : void MeshBase::remove_disjoint_boundary_pair(const boundary_id_type b1,
2267 : const boundary_id_type b2)
2268 : {
2269 : // Nothing to remove if not allocated or empty
2270 21 : if (!_disjoint_neighbor_boundary_pairs || _disjoint_neighbor_boundary_pairs->empty())
2271 0 : return;
2272 :
2273 6 : auto & pairs = *_disjoint_neighbor_boundary_pairs;
2274 :
2275 : // Helper to check and erase both directions
2276 18 : auto erase_if_match = [](boundary_id_type key,
2277 : boundary_id_type pair,
2278 24 : PeriodicBoundaries & pb_map)
2279 : {
2280 12 : auto it = pb_map.find(key);
2281 42 : if (it != pb_map.end())
2282 : {
2283 12 : const auto & pb = *(it->second);
2284 : // Check both directions
2285 42 : if ((pb.myboundary == key && pb.pairedboundary == pair) ||
2286 0 : (pb.pairedboundary == key && pb.myboundary == pair))
2287 30 : pb_map.erase(it);
2288 : }
2289 42 : };
2290 :
2291 21 : erase_if_match(b1, b2, pairs);
2292 21 : erase_if_match(b2, b1, pairs);
2293 : }
2294 :
2295 :
2296 : #endif
2297 :
2298 :
2299 :
2300 0 : void MeshBase::set_point_locator_close_to_point_tol(Real val)
2301 : {
2302 0 : _point_locator_close_to_point_tol = val;
2303 0 : if (_point_locator)
2304 : {
2305 0 : if (val > 0.)
2306 0 : _point_locator->set_close_to_point_tol(val);
2307 : else
2308 0 : _point_locator->unset_close_to_point_tol();
2309 : }
2310 0 : }
2311 :
2312 :
2313 :
2314 42 : Real MeshBase::get_point_locator_close_to_point_tol() const
2315 : {
2316 42 : return _point_locator_close_to_point_tol;
2317 : }
2318 :
2319 :
2320 :
2321 535 : void MeshBase::size_elem_extra_integers()
2322 : {
2323 308 : const std::size_t new_size = _elem_integer_names.size();
2324 :
2325 : Threads::parallel_for
2326 535 : (this->element_stored_range(),
2327 1372 : [new_size, this](const ElemRange & range)
2328 : {
2329 1872 : for (Elem * elem : range)
2330 1337 : elem->add_extra_integers(new_size, this->_elem_integer_default_values);
2331 381 : });
2332 535 : }
2333 :
2334 :
2335 :
2336 1714 : void MeshBase::size_node_extra_integers()
2337 : {
2338 966 : const std::size_t new_size = _node_integer_names.size();
2339 45875 : for (auto node : this->node_ptr_range())
2340 30799 : node->add_extra_integers(new_size, _node_integer_default_values);
2341 1714 : }
2342 :
2343 :
2344 : std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
2345 2807 : MeshBase::merge_extra_integer_names(const MeshBase & other)
2346 : {
2347 854 : std::pair<std::vector<unsigned int>, std::vector<unsigned int>> returnval;
2348 2807 : returnval.first = this->add_elem_integers(other._elem_integer_names, true, &other._elem_integer_default_values);
2349 2807 : returnval.second = this->add_node_integers(other._node_integer_names, true, &other._node_integer_default_values);
2350 2807 : return returnval;
2351 : }
2352 :
2353 :
2354 :
2355 : void
2356 42 : MeshBase::post_dofobject_moves(MeshBase && other_mesh)
2357 : {
2358 : // Now that all the DofObject moving is done, we can move the GhostingFunctor objects
2359 : // which include the _default_ghosting,_ghosting_functors and _shared_functors. We also need
2360 : // to set the mesh object associated with these functors to the assignee mesh.
2361 :
2362 : // _default_ghosting
2363 12 : _default_ghosting = std::move(other_mesh._default_ghosting);
2364 42 : _default_ghosting->set_mesh(this);
2365 :
2366 : // _ghosting_functors
2367 42 : _ghosting_functors = std::move(other_mesh._ghosting_functors);
2368 :
2369 84 : for (const auto gf : _ghosting_functors )
2370 : {
2371 42 : gf->set_mesh(this);
2372 : }
2373 :
2374 : // _shared_functors
2375 12 : _shared_functors = std::move(other_mesh._shared_functors);
2376 :
2377 42 : for (const auto & sf : _shared_functors )
2378 : {
2379 0 : (sf.second)->set_mesh(this);
2380 : }
2381 :
2382 : // _constraint_rows
2383 12 : _constraint_rows = std::move(other_mesh._constraint_rows);
2384 :
2385 42 : if (other_mesh.partitioner())
2386 42 : _partitioner = std::move(other_mesh.partitioner());
2387 42 : }
2388 :
2389 :
2390 : void
2391 0 : MeshBase::copy_cached_data(const MeshBase & other_mesh)
2392 : {
2393 0 : this->_spatial_dimension = other_mesh._spatial_dimension;
2394 0 : this->_elem_dims = other_mesh._elem_dims;
2395 0 : this->_elem_default_orders = other_mesh._elem_default_orders;
2396 0 : this->_supported_nodal_order = other_mesh._supported_nodal_order;
2397 0 : this->_mesh_subdomains = other_mesh._mesh_subdomains;
2398 0 : }
2399 :
2400 :
2401 2053 : bool MeshBase::nodes_and_elements_equal(const MeshBase & other_mesh) const
2402 : {
2403 378586 : for (const auto & other_node : other_mesh.node_ptr_range())
2404 : {
2405 312721 : const Node * node = this->query_node_ptr(other_node->id());
2406 312721 : if (!node)
2407 0 : return false;
2408 312721 : if (*other_node != *node)
2409 0 : return false;
2410 607 : }
2411 378586 : for (const auto & node : this->node_ptr_range())
2412 312721 : if (!other_mesh.query_node_ptr(node->id()))
2413 607 : return false;
2414 :
2415 400146 : for (const auto & other_elem : other_mesh.element_ptr_range())
2416 : {
2417 375164 : const Elem * elem = this->query_elem_ptr(other_elem->id());
2418 375164 : if (!elem)
2419 0 : return false;
2420 375164 : if (!other_elem->topologically_equal(*elem))
2421 6 : return false;
2422 607 : }
2423 400110 : for (const auto & elem : this->element_ptr_range())
2424 375143 : if (!other_mesh.query_elem_ptr(elem->id()))
2425 598 : return false;
2426 :
2427 2032 : return true;
2428 : }
2429 :
2430 :
2431 764 : dof_id_type MeshBase::n_constraint_rows() const
2432 : {
2433 764 : dof_id_type n_local_rows=0, n_unpartitioned_rows=0;
2434 764 : for (const auto & [node, node_constraints] : _constraint_rows)
2435 : {
2436 : // Unpartitioned nodes
2437 0 : if (node->processor_id() == DofObject::invalid_processor_id)
2438 0 : n_unpartitioned_rows++;
2439 0 : else if (node->processor_id() == this->processor_id())
2440 0 : n_local_rows++;
2441 : }
2442 :
2443 764 : this->comm().sum(n_local_rows);
2444 :
2445 764 : return n_unpartitioned_rows + n_local_rows;
2446 : }
2447 :
2448 :
2449 : void
2450 2632 : MeshBase::copy_constraint_rows(const MeshBase & other_mesh)
2451 : {
2452 1608 : LOG_SCOPE("copy_constraint_rows(mesh)", "MeshBase");
2453 :
2454 1592 : _constraint_rows.clear();
2455 :
2456 804 : const auto & other_constraint_rows = other_mesh.get_constraint_rows();
2457 14679 : for (const auto & [other_node, other_node_constraints] : other_constraint_rows)
2458 : {
2459 12047 : const Node * const our_node = this->node_ptr(other_node->id());
2460 6416 : constraint_rows_mapped_type our_node_constraints;
2461 44238 : for (const auto & [other_inner_key_pair, constraint_value] : other_node_constraints)
2462 : {
2463 9048 : const auto & [other_elem, local_node_id] = other_inner_key_pair;
2464 32191 : const Elem * const our_elem = this->elem_ptr(other_elem->id());
2465 32191 : our_node_constraints.emplace_back(std::make_pair(our_elem, local_node_id), constraint_value);
2466 : }
2467 12047 : _constraint_rows[our_node] = std::move(our_node_constraints);
2468 : }
2469 2632 : }
2470 :
2471 :
2472 : template <typename T>
2473 : void
2474 53 : MeshBase::copy_constraint_rows(const SparseMatrix<T> & constraint_operator,
2475 : bool precondition_constraint_operator)
2476 : {
2477 32 : LOG_SCOPE("copy_constraint_rows(mat)", "MeshBase");
2478 :
2479 16 : this->_constraint_rows.clear();
2480 :
2481 : // We're not going to support doing this distributed yet; it'd be
2482 : // pointless unless we temporarily had a linear partitioning to
2483 : // better match the constraint operator.
2484 85 : MeshSerializer serialize(*this);
2485 :
2486 : // Our current mesh should already reflect the desired assembly space
2487 53 : libmesh_error_msg_if(this->n_nodes() != constraint_operator.m(),
2488 : "Constraint operator matrix with " <<
2489 : constraint_operator.m() <<
2490 : "rows does not match this mesh with " <<
2491 : this->n_nodes() << " nodes");
2492 :
2493 : // First, find what new unconstrained DoFs we need to add. We can't
2494 : // iterate over columns in a SparseMatrix, so we'll iterate over
2495 : // rows and keep track of columns.
2496 :
2497 : // If we have nodes that will work unconstrained, keep track of
2498 : // their node ids and corresponding column indices.
2499 : // existing_unconstrained_nodes[column_id] = node_id
2500 32 : std::map<dof_id_type, dof_id_type> existing_unconstrained_columns;
2501 32 : std::set<dof_id_type> existing_unconstrained_nodes;
2502 :
2503 : // In case we need new nodes, keep track of their columns.
2504 : // columns[j][k] will be the kth row index and value of column j
2505 : typedef
2506 : std::unordered_map<dof_id_type,
2507 : std::vector<std::pair<dof_id_type, Real>>>
2508 : columns_type;
2509 90 : columns_type columns(constraint_operator.n());
2510 :
2511 : // If we need to precondition the constraint operator (e.g. it's an
2512 : // unpreconditioned extraction operator for a Flex IGA matrix),
2513 : // we'll want to keep track of the sum of each column, because we'll
2514 : // be dividing each column by that sum (Jacobi preconditioning on
2515 : // the right, which then leads to symmetric preconditioning on a
2516 : // physics Jacobian).
2517 32 : std::unordered_map<dof_id_type, Real> column_sums;
2518 :
2519 : // Work in parallel, though we'll have to sync shortly
2520 1044 : for (auto i : make_range(constraint_operator.row_start(),
2521 21 : constraint_operator.row_stop()))
2522 : {
2523 618 : std::vector<numeric_index_type> indices;
2524 618 : std::vector<T> values;
2525 :
2526 991 : constraint_operator.get_row(i, indices, values);
2527 309 : libmesh_assert_equal_to(indices.size(), values.size());
2528 :
2529 1375 : if (indices.size() == 1 &&
2530 225 : values[0] == T(1))
2531 : {
2532 : // If we have multiple simple Ui=Uj constraints, let the
2533 : // first one be our "unconstrained" node and let the others
2534 : // be constrained to it.
2535 258 : if (existing_unconstrained_columns.find(indices[0]) !=
2536 138 : existing_unconstrained_columns.end())
2537 : {
2538 12 : const auto j = indices[0];
2539 12 : columns[j].emplace_back(i, 1);
2540 : }
2541 : else
2542 : {
2543 180 : existing_unconstrained_nodes.insert(i);
2544 246 : existing_unconstrained_columns.emplace(indices[0],i);
2545 : }
2546 : }
2547 : else
2548 6036 : for (auto jj : index_range(indices))
2549 : {
2550 5303 : const auto j = indices[jj];
2551 7062 : const Real coef = libmesh_real(values[jj]);
2552 1759 : libmesh_assert_equal_to(coef, values[jj]);
2553 5303 : columns[j].emplace_back(i, coef);
2554 : }
2555 : }
2556 :
2557 : // Merge data from different processors' slabs of the matrix
2558 53 : this->comm().set_union(existing_unconstrained_nodes);
2559 53 : this->comm().set_union(existing_unconstrained_columns);
2560 :
2561 48 : std::vector<columns_type> all_columns;
2562 53 : this->comm().allgather(columns, all_columns);
2563 :
2564 16 : columns.clear();
2565 154 : for (auto p : index_range(all_columns))
2566 1975 : for (auto & [j, subcol] : all_columns[p])
2567 12475 : for (auto [i, v] : subcol)
2568 10601 : columns[j].emplace_back(i,v);
2569 :
2570 : // Keep track of elements on which unconstrained nodes exist, and
2571 : // their local node indices.
2572 : // node_to_elem_ptrs[node] = [elem_id, local_node_num]
2573 32 : std::unordered_map<const Node *, std::pair<dof_id_type, unsigned int>> node_to_elem_ptrs;
2574 :
2575 : // Find elements attached to any existing nodes that will stay
2576 : // unconstrained. We'll also build a subdomain set here so we don't
2577 : // have to assert that the mesh is already prepared before we pick a
2578 : // new subdomain for any NodeElems we need to add.
2579 32 : std::set<subdomain_id_type> subdomain_ids;
2580 7325 : for (const Elem * elem : this->element_ptr_range())
2581 : {
2582 4298 : subdomain_ids.insert(elem->subdomain_id());
2583 21084 : for (auto n : make_range(elem->n_nodes()))
2584 : {
2585 16786 : const Node * node = elem->node_ptr(n);
2586 16786 : if (existing_unconstrained_nodes.count(node->id()))
2587 1103 : node_to_elem_ptrs.emplace(node, std::make_pair(elem->id(), n));
2588 : }
2589 : }
2590 :
2591 53 : const subdomain_id_type new_sbd_id = *subdomain_ids.rbegin() + 1;
2592 :
2593 1557 : for (auto j : make_range(constraint_operator.n()))
2594 : {
2595 : // If we already have a good node for this then we're done
2596 576 : if (existing_unconstrained_columns.count(j))
2597 444 : continue;
2598 :
2599 : // Get a half-decent spot to place a new NodeElem for
2600 : // unconstrained DoF(s) here. Getting a *fully*-decent spot
2601 : // would require finding a Moore-Penrose pseudoinverse, and I'm
2602 : // not going to do that, but scaling a transpose will at least
2603 : // get us a little uniqueness to make visualization reasonable.
2604 264 : Point newpt;
2605 264 : Real total_scaling = 0;
2606 264 : unsigned int total_entries = 0;
2607 :
2608 : // We'll get a decent initial pid choice here too, if only to
2609 : // aid in later repartitioning.
2610 528 : std::map<processor_id_type, int> pids;
2611 :
2612 264 : auto & column = columns[j];
2613 10486 : for (auto [i, r] : column)
2614 : {
2615 9692 : Node & constrained_node = this->node_ref(i);
2616 9692 : const Point constrained_pt = constrained_node;
2617 3228 : newpt += r*constrained_pt;
2618 9692 : total_scaling += r;
2619 9692 : ++total_entries;
2620 9692 : ++pids[constrained_node.processor_id()];
2621 : }
2622 :
2623 794 : if (precondition_constraint_operator)
2624 0 : column_sums[j] = total_scaling;
2625 :
2626 794 : libmesh_error_msg_if
2627 : (!total_entries,
2628 : "Empty column " << j <<
2629 : " found in constraint operator matrix");
2630 :
2631 : // If we have *cancellation* here then we can end up dividing by
2632 : // zero; try just evenly scaling across all constrained node
2633 : // points instead.
2634 794 : if (total_scaling > TOLERANCE)
2635 264 : newpt /= total_scaling;
2636 : else
2637 0 : newpt /= total_entries;
2638 :
2639 794 : Node *n = this->add_point(newpt);
2640 1058 : std::unique_ptr<Elem> elem = Elem::build(NODEELEM);
2641 794 : elem->set_node(0, n);
2642 794 : elem->subdomain_id() = new_sbd_id;
2643 :
2644 1058 : Elem * added_elem = this->add_elem(std::move(elem));
2645 794 : this->_elem_dims.insert(0);
2646 794 : this->_elem_default_orders.insert(added_elem->default_order());
2647 794 : this->_supported_nodal_order =
2648 266 : static_cast<Order>
2649 1588 : (std::min(static_cast<int>(this->_supported_nodal_order),
2650 794 : static_cast<int>(added_elem->supported_nodal_order())));
2651 530 : this->_mesh_subdomains.insert(new_sbd_id);
2652 794 : node_to_elem_ptrs.emplace(n, std::make_pair(added_elem->id(), 0));
2653 794 : existing_unconstrained_columns.emplace(j,n->id());
2654 :
2655 : // Repartition the new objects *after* adding them, so a
2656 : // DistributedMesh doesn't get confused and think you're not
2657 : // adding them on all processors at once.
2658 264 : int n_pids = 0;
2659 2098 : for (auto [pid, count] : pids)
2660 1304 : if (count >= n_pids)
2661 : {
2662 324 : n_pids = count;
2663 974 : added_elem->processor_id() = pid;
2664 974 : n->processor_id() = pid;
2665 : }
2666 : }
2667 :
2668 : // Calculate constraint rows in an indexed form that's easy for us
2669 : // to allgather
2670 : std::unordered_map<dof_id_type,
2671 : std::vector<std::pair<std::pair<dof_id_type, unsigned int>,Real>>>
2672 32 : indexed_constraint_rows;
2673 :
2674 1044 : for (auto i : make_range(constraint_operator.row_start(),
2675 21 : constraint_operator.row_stop()))
2676 : {
2677 423 : if (existing_unconstrained_nodes.count(i))
2678 246 : continue;
2679 :
2680 486 : std::vector<numeric_index_type> indices;
2681 486 : std::vector<T> values;
2682 :
2683 745 : constraint_operator.get_row(i, indices, values);
2684 :
2685 486 : std::vector<std::pair<std::pair<dof_id_type, unsigned int>, Real>> constraint_row;
2686 :
2687 6060 : for (auto jj : index_range(indices))
2688 : {
2689 7077 : const dof_id_type node_id =
2690 : existing_unconstrained_columns[indices[jj]];
2691 :
2692 5315 : Node & constraining_node = this->node_ref(node_id);
2693 :
2694 1762 : libmesh_assert(node_to_elem_ptrs.count(&constraining_node));
2695 :
2696 5315 : auto p = node_to_elem_ptrs[&constraining_node];
2697 :
2698 5286 : Real coef = libmesh_real(values[jj]);
2699 1762 : libmesh_assert_equal_to(coef, values[jj]);
2700 :
2701 : // If we're preconditioning and we created a nodeelem then
2702 : // we can scale the meaning of that nodeelem's value to give
2703 : // us a better-conditioned matrix after the constraints are
2704 : // applied.
2705 5315 : if (precondition_constraint_operator)
2706 0 : if (auto sum_it = column_sums.find(indices[jj]);
2707 0 : sum_it != column_sums.end())
2708 : {
2709 0 : const Real scaling = sum_it->second;
2710 :
2711 0 : if (scaling > TOLERANCE)
2712 0 : coef /= scaling;
2713 : }
2714 :
2715 7077 : constraint_row.emplace_back(std::make_pair(p, coef));
2716 : }
2717 :
2718 243 : indexed_constraint_rows.emplace(i, std::move(constraint_row));
2719 : }
2720 :
2721 53 : this->comm().set_union(indexed_constraint_rows);
2722 :
2723 : // Add constraint rows as mesh constraint rows
2724 1527 : for (auto & [node_id, indexed_row] : indexed_constraint_rows)
2725 : {
2726 1474 : Node * constrained_node = this->node_ptr(node_id);
2727 :
2728 972 : constraint_rows_mapped_type constraint_row;
2729 :
2730 12075 : for (auto [p, coef] : indexed_row)
2731 : {
2732 10601 : const Elem * elem = this->elem_ptr(p.first);
2733 7048 : constraint_row.emplace_back
2734 14125 : (std::make_pair(std::make_pair(elem, p.second), coef));
2735 : }
2736 :
2737 988 : this->_constraint_rows.emplace(constrained_node,
2738 486 : std::move(constraint_row));
2739 : }
2740 74 : }
2741 :
2742 :
2743 0 : void MeshBase::print_constraint_rows(std::ostream & os,
2744 : bool print_nonlocal) const
2745 : {
2746 0 : parallel_object_only();
2747 :
2748 : std::string local_constraints =
2749 0 : this->get_local_constraints(print_nonlocal);
2750 :
2751 0 : if (this->processor_id())
2752 : {
2753 0 : this->comm().send(0, local_constraints);
2754 : }
2755 : else
2756 : {
2757 0 : os << "Processor 0:\n";
2758 0 : os << local_constraints;
2759 :
2760 0 : for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
2761 : {
2762 0 : this->comm().receive(p, local_constraints);
2763 0 : os << "Processor " << p << ":\n";
2764 0 : os << local_constraints;
2765 : }
2766 : }
2767 0 : }
2768 :
2769 :
2770 :
2771 0 : std::string MeshBase::get_local_constraints(bool print_nonlocal) const
2772 : {
2773 0 : std::ostringstream os;
2774 :
2775 0 : if (print_nonlocal)
2776 0 : os << "All ";
2777 : else
2778 0 : os << "Local ";
2779 :
2780 0 : os << "Mesh Constraint Rows:"
2781 0 : << std::endl;
2782 :
2783 0 : for (const auto & [node, row] : _constraint_rows)
2784 : {
2785 0 : const bool local = (node->processor_id() == this->processor_id());
2786 :
2787 : // Skip non-local dofs if requested
2788 0 : if (!print_nonlocal && !local)
2789 0 : continue;
2790 :
2791 0 : os << "Constraints for " << (local ? "Local" : "Ghost") << " Node " << node->id()
2792 0 : << ": \t";
2793 :
2794 0 : for (const auto & [elem_and_node, coef] : row)
2795 0 : os << " ((" << elem_and_node.first->id() << ',' << elem_and_node.second << "), " << coef << ")\t";
2796 :
2797 0 : os << std::endl;
2798 : }
2799 :
2800 0 : return os.str();
2801 0 : }
2802 :
2803 33396 : MeshBase::Preparation::Preparation() :
2804 14080 : is_partitioned(false),
2805 14080 : has_synched_id_counts(false),
2806 14080 : has_neighbor_ptrs(false),
2807 14080 : has_cached_elem_data(false),
2808 14080 : has_interior_parent_ptrs(false),
2809 14080 : has_removed_remote_elements(false),
2810 14080 : has_removed_orphaned_nodes(false),
2811 14080 : has_boundary_id_sets(false),
2812 33396 : has_reinit_ghosting_functors(false)
2813 33396 : {}
2814 :
2815 46355 : MeshBase::Preparation::operator bool() const
2816 : {
2817 89770 : return is_partitioned &&
2818 43415 : has_synched_id_counts &&
2819 43415 : has_neighbor_ptrs &&
2820 43415 : has_cached_elem_data &&
2821 43415 : has_interior_parent_ptrs &&
2822 43415 : has_removed_remote_elements &&
2823 43415 : has_removed_orphaned_nodes &&
2824 131010 : has_reinit_ghosting_functors &&
2825 85533 : has_boundary_id_sets;
2826 : }
2827 :
2828 : MeshBase::Preparation &
2829 148143 : MeshBase::Preparation::operator= (bool set_all)
2830 : {
2831 148143 : is_partitioned = set_all;
2832 148143 : has_synched_id_counts = set_all;
2833 148143 : has_neighbor_ptrs = set_all;
2834 148143 : has_cached_elem_data = set_all;
2835 148143 : has_interior_parent_ptrs = set_all;
2836 148143 : has_removed_remote_elements = set_all;
2837 148143 : has_removed_orphaned_nodes = set_all;
2838 148143 : has_reinit_ghosting_functors = set_all;
2839 148143 : has_boundary_id_sets = set_all;
2840 :
2841 148143 : return *this;
2842 : }
2843 :
2844 : bool
2845 2095 : MeshBase::Preparation::operator== (const Preparation & other) const
2846 : {
2847 2095 : if (is_partitioned != other.is_partitioned)
2848 0 : return false;
2849 2095 : if (has_synched_id_counts != other.has_synched_id_counts)
2850 0 : return false;
2851 2095 : if (has_neighbor_ptrs != other.has_neighbor_ptrs)
2852 0 : return false;
2853 2095 : if (has_cached_elem_data != other.has_cached_elem_data)
2854 0 : return false;
2855 2095 : if (has_interior_parent_ptrs != other.has_interior_parent_ptrs)
2856 0 : return false;
2857 2095 : if (has_removed_remote_elements != other.has_removed_remote_elements)
2858 0 : return false;
2859 2095 : if (has_removed_orphaned_nodes != other.has_removed_orphaned_nodes)
2860 0 : return false;
2861 2095 : if (has_reinit_ghosting_functors != other.has_reinit_ghosting_functors)
2862 0 : return false;
2863 2095 : if (has_boundary_id_sets != other.has_boundary_id_sets)
2864 0 : return false;
2865 :
2866 1062 : return true;
2867 : }
2868 :
2869 : bool
2870 2095 : MeshBase::Preparation::operator!= (const Preparation & other) const
2871 : {
2872 2095 : return !(*this == other);
2873 : }
2874 :
2875 :
2876 : // Explicit instantiations for our template function
2877 : template LIBMESH_EXPORT void
2878 : MeshBase::copy_constraint_rows(const SparseMatrix<Real> & constraint_operator,
2879 : bool precondition_constraint_operator);
2880 :
2881 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
2882 : template LIBMESH_EXPORT void
2883 : MeshBase::copy_constraint_rows(const SparseMatrix<Complex> & constraint_operator,
2884 : bool precondition_constraint_operator);
2885 : #endif
2886 :
2887 :
2888 : } // namespace libMesh
|
