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