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