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