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