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