Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 :
45 : // C++ includes
46 : #include <algorithm> // for std::min
47 : #include <map> // for std::multimap
48 : #include <memory>
49 : #include <sstream> // for std::ostringstream
50 : #include <unordered_map>
51 :
52 : namespace libMesh
53 : {
54 :
55 :
56 :
57 : // ------------------------------------------------------------
58 : // MeshBase class member functions
59 302284 : MeshBase::MeshBase (const Parallel::Communicator & comm_in,
60 302284 : unsigned char d) :
61 : ParallelObject (comm_in),
62 302284 : boundary_info (new BoundaryInfo(*this)), // BoundaryInfo has protected ctor, can't use std::make_unique
63 284640 : _n_parts (1),
64 284640 : _default_mapping_type(LAGRANGE_MAP),
65 284640 : _default_mapping_data(0),
66 284640 : _is_prepared (false),
67 284640 : _point_locator (),
68 284640 : _count_lower_dim_elems_in_point_locator(true),
69 284640 : _partitioner (),
70 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
71 284640 : _next_unique_id(DofObject::invalid_unique_id),
72 : #endif
73 284640 : _interior_mesh(this),
74 284640 : _skip_noncritical_partitioning(false),
75 311106 : _skip_all_partitioning(libMesh::on_command_line("--skip-partitioning")),
76 284640 : _skip_renumber_nodes_and_elements(false),
77 284640 : _skip_find_neighbors(false),
78 284640 : _allow_remote_element_removal(true),
79 284640 : _spatial_dimension(d),
80 302284 : _default_ghosting(std::make_unique<GhostPointNeighbors>(*this)),
81 968606 : _point_locator_close_to_point_tol(0.)
82 : {
83 293462 : _elem_dims.insert(d);
84 302284 : _ghosting_functors.push_back(_default_ghosting.get());
85 8822 : libmesh_assert_less_equal (LIBMESH_DIM, 3);
86 8822 : libmesh_assert_greater_equal (LIBMESH_DIM, d);
87 8822 : libmesh_assert (libMesh::initialized());
88 302284 : }
89 :
90 :
91 :
92 21765 : MeshBase::MeshBase (const MeshBase & other_mesh) :
93 : ParallelObject (other_mesh),
94 21765 : boundary_info (new BoundaryInfo(*this)), // BoundaryInfo has protected ctor, can't use std::make_unique
95 21765 : _n_parts (other_mesh._n_parts),
96 21765 : _default_mapping_type(other_mesh._default_mapping_type),
97 21765 : _default_mapping_data(other_mesh._default_mapping_data),
98 21765 : _is_prepared (other_mesh._is_prepared),
99 20321 : _point_locator (),
100 21765 : _count_lower_dim_elems_in_point_locator(other_mesh._count_lower_dim_elems_in_point_locator),
101 20321 : _partitioner (),
102 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
103 21765 : _next_unique_id(other_mesh._next_unique_id),
104 : #endif
105 : // If the other mesh interior_parent pointers just go back to
106 : // itself, so should we
107 21765 : _interior_mesh((other_mesh._interior_mesh == &other_mesh) ?
108 : this : other_mesh._interior_mesh),
109 21765 : _skip_noncritical_partitioning(other_mesh._skip_noncritical_partitioning),
110 21765 : _skip_all_partitioning(other_mesh._skip_all_partitioning),
111 21765 : _skip_renumber_nodes_and_elements(other_mesh._skip_renumber_nodes_and_elements),
112 21765 : _skip_find_neighbors(other_mesh._skip_find_neighbors),
113 21765 : _allow_remote_element_removal(other_mesh._allow_remote_element_removal),
114 730 : _elem_dims(other_mesh._elem_dims),
115 730 : _elem_default_orders(other_mesh._elem_default_orders),
116 21765 : _supported_nodal_order(other_mesh._supported_nodal_order),
117 730 : _elemset_codes_inverse_map(other_mesh._elemset_codes_inverse_map),
118 730 : _all_elemset_ids(other_mesh._all_elemset_ids),
119 21765 : _spatial_dimension(other_mesh._spatial_dimension),
120 21765 : _default_ghosting(std::make_unique<GhostPointNeighbors>(*this)),
121 89202 : _point_locator_close_to_point_tol(other_mesh._point_locator_close_to_point_tol)
122 : {
123 730 : const GhostingFunctor * const other_default_ghosting = other_mesh._default_ghosting.get();
124 :
125 74214 : for (GhostingFunctor * const gf : other_mesh._ghosting_functors)
126 : {
127 : // If the other mesh is using default ghosting, then we will use our own
128 : // default ghosting
129 52449 : if (gf == other_default_ghosting)
130 : {
131 21765 : _ghosting_functors.push_back(_default_ghosting.get());
132 21765 : continue;
133 : }
134 :
135 60352 : std::shared_ptr<GhostingFunctor> clone_gf = gf->clone();
136 : // Some subclasses of GhostingFunctor might not override the
137 : // clone function yet. If this is the case, GhostingFunctor will
138 : // return nullptr by default. The clone function should be overridden
139 : // in all derived classes. This following code ("else") is written
140 : // for API upgrade. That will allow users gradually to update their code.
141 : // Once the API upgrade is done, we will come back and delete "else."
142 30684 : if (clone_gf)
143 : {
144 30684 : clone_gf->set_mesh(this);
145 60352 : add_ghosting_functor(clone_gf);
146 : }
147 : else
148 : {
149 : libmesh_deprecated();
150 0 : add_ghosting_functor(*gf);
151 : }
152 : }
153 :
154 21765 : if (other_mesh._partitioner.get())
155 42800 : _partitioner = other_mesh._partitioner->clone();
156 :
157 : // _elemset_codes stores pointers to entries in _elemset_codes_inverse_map,
158 : // so it is not possible to simply copy it directly from other_mesh
159 21765 : for (const auto & [set, code] : _elemset_codes_inverse_map)
160 0 : _elemset_codes.emplace(code, &set);
161 21765 : }
162 :
163 426 : MeshBase& MeshBase::operator= (MeshBase && other_mesh)
164 : {
165 12 : LOG_SCOPE("operator=(&&)", "MeshBase");
166 :
167 : // Move assign as a ParallelObject.
168 12 : this->ParallelObject::operator=(other_mesh);
169 :
170 426 : _n_parts = other_mesh.n_partitions();
171 426 : _default_mapping_type = other_mesh.default_mapping_type();
172 426 : _default_mapping_data = other_mesh.default_mapping_data();
173 426 : _is_prepared = other_mesh.is_prepared();
174 12 : _point_locator = std::move(other_mesh._point_locator);
175 426 : _count_lower_dim_elems_in_point_locator = other_mesh.get_count_lower_dim_elems_in_point_locator();
176 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
177 426 : _next_unique_id = other_mesh.next_unique_id();
178 : #endif
179 : // If the other mesh interior_parent pointers just go back to
180 : // itself, so should we
181 426 : _interior_mesh = (other_mesh._interior_mesh == &other_mesh) ?
182 : this : other_mesh._interior_mesh;
183 426 : _skip_noncritical_partitioning = other_mesh.skip_noncritical_partitioning();
184 426 : _skip_all_partitioning = other_mesh.skip_partitioning();
185 426 : _skip_renumber_nodes_and_elements = !(other_mesh.allow_renumbering());
186 426 : _skip_find_neighbors = !(other_mesh.allow_find_neighbors());
187 426 : _allow_remote_element_removal = other_mesh.allow_remote_element_removal();
188 12 : _block_id_to_name = std::move(other_mesh._block_id_to_name);
189 12 : _elem_dims = std::move(other_mesh.elem_dimensions());
190 12 : _elem_default_orders = std::move(other_mesh.elem_default_orders());
191 426 : _supported_nodal_order = other_mesh.supported_nodal_order();
192 12 : _elemset_codes = std::move(other_mesh._elemset_codes);
193 12 : _elemset_codes_inverse_map = std::move(other_mesh._elemset_codes_inverse_map);
194 12 : _all_elemset_ids = std::move(other_mesh._all_elemset_ids);
195 426 : _spatial_dimension = other_mesh.spatial_dimension();
196 426 : _elem_integer_names = std::move(other_mesh._elem_integer_names);
197 426 : _elem_integer_default_values = std::move(other_mesh._elem_integer_default_values);
198 426 : _node_integer_names = std::move(other_mesh._node_integer_names);
199 426 : _node_integer_default_values = std::move(other_mesh._node_integer_default_values);
200 426 : _point_locator_close_to_point_tol = other_mesh.get_point_locator_close_to_point_tol();
201 :
202 : // This relies on our subclasses *not* invalidating pointers when we
203 : // do their portion of the move assignment later!
204 12 : boundary_info = std::move(other_mesh.boundary_info);
205 12 : boundary_info->set_mesh(*this);
206 :
207 : #ifdef DEBUG
208 : // Make sure that move assignment worked for pointers
209 12 : for (const auto & [set, code] : _elemset_codes_inverse_map)
210 : {
211 0 : auto it = _elemset_codes.find(code);
212 0 : libmesh_assert_msg(it != _elemset_codes.end(),
213 : "Elemset code " << code << " not found in _elmset_codes container.");
214 0 : libmesh_assert_equal_to(it->second, &set);
215 : }
216 : #endif
217 :
218 : // We're *not* really done at this point, but we have the problem
219 : // that some of our data movement might be expecting subclasses data
220 : // movement to happen first. We'll let subclasses handle that by
221 : // calling our post_dofobject_moves()
222 438 : return *this;
223 : }
224 :
225 :
226 6250 : bool MeshBase::operator== (const MeshBase & other_mesh) const
227 : {
228 786 : LOG_SCOPE("operator==()", "MeshBase");
229 :
230 6250 : bool is_equal = this->locally_equals(other_mesh);
231 6250 : this->comm().min(is_equal);
232 7036 : return is_equal;
233 : }
234 :
235 :
236 6250 : bool MeshBase::locally_equals (const MeshBase & other_mesh) const
237 : {
238 : // Check whether (almost) everything in the base is equal
239 : //
240 : // We don't check _next_unique_id here, because it's expected to
241 : // change in a DistributedMesh prepare_for_use(); it's conceptually
242 : // "mutable".
243 : //
244 : // We use separate if statements instead of logical operators here,
245 : // to make it easy to see the failing condition when using a
246 : // debugger to figure out why a MeshTools::valid_is_prepared(mesh)
247 : // is failing.
248 6250 : if (_n_parts != other_mesh._n_parts)
249 0 : return false;
250 6250 : if (_default_mapping_type != other_mesh._default_mapping_type)
251 0 : return false;
252 6250 : if (_default_mapping_data != other_mesh._default_mapping_data)
253 0 : return false;
254 6250 : if (_is_prepared != other_mesh._is_prepared)
255 0 : return false;
256 7194 : if (_count_lower_dim_elems_in_point_locator !=
257 6250 : other_mesh._count_lower_dim_elems_in_point_locator)
258 0 : return false;
259 :
260 : // We should either both have our own interior parents or both not;
261 : // but if we both don't then we can't really assert anything else
262 : // because pointing at the same interior mesh is fair but so is
263 : // pointing at two different copies of "the same" interior mesh.
264 7194 : if ((_interior_mesh == this) !=
265 6250 : (other_mesh._interior_mesh == &other_mesh))
266 0 : return false;
267 :
268 6250 : if (_skip_noncritical_partitioning != other_mesh._skip_noncritical_partitioning)
269 0 : return false;
270 6250 : if (_skip_all_partitioning != other_mesh._skip_all_partitioning)
271 0 : return false;
272 6250 : if (_skip_renumber_nodes_and_elements != other_mesh._skip_renumber_nodes_and_elements)
273 0 : return false;
274 6250 : if (_skip_find_neighbors != other_mesh._skip_find_neighbors)
275 0 : return false;
276 6250 : if (_allow_remote_element_removal != other_mesh._allow_remote_element_removal)
277 0 : return false;
278 6250 : if (_spatial_dimension != other_mesh._spatial_dimension)
279 0 : return false;
280 6250 : if (_point_locator_close_to_point_tol != other_mesh._point_locator_close_to_point_tol)
281 0 : return false;
282 6250 : if (_block_id_to_name != other_mesh._block_id_to_name)
283 0 : return false;
284 6250 : if (_elem_dims != other_mesh._elem_dims)
285 0 : return false;
286 6250 : if (_elem_default_orders != other_mesh._elem_default_orders)
287 0 : return false;
288 6250 : if (_supported_nodal_order != other_mesh._supported_nodal_order)
289 0 : return false;
290 6250 : if (_mesh_subdomains != other_mesh._mesh_subdomains)
291 0 : return false;
292 6250 : if (_all_elemset_ids != other_mesh._all_elemset_ids)
293 0 : return false;
294 6250 : if (_elem_integer_names != other_mesh._elem_integer_names)
295 0 : return false;
296 6408 : if (_elem_integer_default_values != other_mesh._elem_integer_default_values)
297 0 : return false;
298 6250 : if (_node_integer_names != other_mesh._node_integer_names)
299 0 : return false;
300 6408 : if (_node_integer_default_values != other_mesh._node_integer_default_values)
301 0 : return false;
302 6250 : if (bool(_default_ghosting) != bool(other_mesh._default_ghosting))
303 0 : return false;
304 6250 : if (bool(_partitioner) != bool(other_mesh._partitioner))
305 0 : return false;
306 6250 : if (*boundary_info != *other_mesh.boundary_info)
307 0 : return false;
308 :
309 : const constraint_rows_type & other_rows =
310 786 : other_mesh.get_constraint_rows();
311 57886 : for (const auto & [node, row] : this->_constraint_rows)
312 : {
313 51636 : const dof_id_type node_id = node->id();
314 51636 : const Node * other_node = other_mesh.query_node_ptr(node_id);
315 51636 : if (!other_node)
316 0 : return false;
317 :
318 6416 : auto it = other_rows.find(other_node);
319 51636 : if (it == other_rows.end())
320 0 : return false;
321 :
322 6416 : const auto & other_row = it->second;
323 58052 : if (row.size() != other_row.size())
324 0 : return false;
325 :
326 231893 : for (auto i : index_range(row))
327 : {
328 18096 : const auto & [elem_pair, coef] = row[i];
329 18096 : const auto & [other_elem_pair, other_coef] = other_row[i];
330 18096 : libmesh_assert(elem_pair.first);
331 18096 : libmesh_assert(other_elem_pair.first);
332 180257 : if (elem_pair.first->id() !=
333 342418 : other_elem_pair.first->id() ||
334 180257 : elem_pair.second !=
335 198353 : other_elem_pair.second ||
336 180257 : coef != other_coef)
337 0 : return false;
338 : }
339 : }
340 :
341 6250 : for (const auto & [elemset_code, elemset_ptr] : this->_elemset_codes)
342 0 : if (const auto it = other_mesh._elemset_codes.find(elemset_code);
343 0 : it == other_mesh._elemset_codes.end() || *elemset_ptr != *it->second)
344 0 : return false;
345 :
346 : // FIXME: we have no good way to compare ghosting functors, since
347 : // they're in a vector of pointers, and we have no way *at all*
348 : // to compare ghosting functors, since they don't have operator==
349 : // defined and we encourage users to subclass them. We can check if
350 : // we have the same number, is all.
351 7194 : if (_ghosting_functors.size() !=
352 944 : other_mesh._ghosting_functors.size())
353 0 : return false;
354 :
355 : // Same deal for partitioners. We tested that we both have one or
356 : // both don't, but are they equivalent? Let's guess "yes".
357 :
358 : // Now let the subclasses decide whether everything else is equal
359 6250 : return this->subclass_locally_equals(other_mesh);
360 : }
361 :
362 :
363 333585 : MeshBase::~MeshBase()
364 : {
365 324049 : this->MeshBase::clear();
366 :
367 9552 : libmesh_exceptionless_assert (!libMesh::closed());
368 1238932 : }
369 :
370 :
371 :
372 21990178 : unsigned int MeshBase::mesh_dimension() const
373 : {
374 21990178 : if (!_elem_dims.empty())
375 21990163 : return cast_int<unsigned int>(*_elem_dims.rbegin());
376 0 : return 0;
377 : }
378 :
379 :
380 :
381 0 : void MeshBase::set_elem_dimensions(std::set<unsigned char> elem_dims)
382 : {
383 : #ifdef DEBUG
384 : // In debug mode, we call cache_elem_data() and then make sure
385 : // the result actually agrees with what the user specified.
386 0 : parallel_object_only();
387 :
388 0 : this->cache_elem_data();
389 0 : libmesh_assert_msg(_elem_dims == elem_dims, \
390 : "Specified element dimensions does not match true element dimensions!");
391 : #endif
392 :
393 0 : _elem_dims = std::move(elem_dims);
394 0 : }
395 :
396 :
397 :
398 1775 : void MeshBase::add_elemset_code(dof_id_type code, MeshBase::elemset_type id_set)
399 : {
400 : // Populate inverse map, stealing id_set's resources
401 1725 : auto [it1, inserted1] = _elemset_codes_inverse_map.emplace(std::move(id_set), code);
402 :
403 : // Reference to the newly inserted (or previously existing) id_set
404 1775 : const auto & inserted_id_set = it1->first;
405 :
406 : // Keep track of all elemset ids ever added for O(1) n_elemsets()
407 : // performance. Only need to do this if we didn't know about this
408 : // id_set before...
409 1775 : if (inserted1)
410 1775 : _all_elemset_ids.insert(inserted_id_set.begin(), inserted_id_set.end());
411 :
412 : // Take the address of the newly emplaced set to use in
413 : // _elemset_codes, avoid duplicating std::set storage
414 1775 : auto [it2, inserted2] = _elemset_codes.emplace(code, &inserted_id_set);
415 :
416 : // Throw an error if this code already exists with a pointer to a
417 : // different set of ids.
418 1775 : libmesh_error_msg_if(!inserted2 && it2->second != &inserted_id_set,
419 : "The elemset code " << code << " already exists with a different id_set.");
420 1775 : }
421 :
422 :
423 :
424 14859 : unsigned int MeshBase::n_elemsets() const
425 : {
426 14859 : return _all_elemset_ids.size();
427 : }
428 :
429 2963 : void MeshBase::get_elemsets(dof_id_type elemset_code, MeshBase::elemset_type & id_set_to_fill) const
430 : {
431 : // If we don't recognize this elemset_code, hand back an empty set
432 102 : id_set_to_fill.clear();
433 :
434 2963 : if (const auto it = _elemset_codes.find(elemset_code);
435 102 : it != _elemset_codes.end())
436 1386 : id_set_to_fill.insert(it->second->begin(), it->second->end());
437 2963 : }
438 :
439 852 : dof_id_type MeshBase::get_elemset_code(const MeshBase::elemset_type & id_set) const
440 : {
441 24 : auto it = _elemset_codes_inverse_map.find(id_set);
442 852 : return (it == _elemset_codes_inverse_map.end()) ? DofObject::invalid_id : it->second;
443 : }
444 :
445 4123 : std::vector<dof_id_type> MeshBase::get_elemset_codes() const
446 : {
447 118 : std::vector<dof_id_type> ret;
448 4123 : ret.reserve(_elemset_codes.size());
449 4904 : for (const auto & pr : _elemset_codes)
450 781 : ret.push_back(pr.first);
451 4123 : return ret;
452 : }
453 :
454 284 : void MeshBase::change_elemset_code(dof_id_type old_code, dof_id_type new_code)
455 : {
456 : // Look up elemset ids for old_code
457 8 : auto it = _elemset_codes.find(old_code);
458 :
459 : // If we don't have the old_code, then do nothing. Alternatively, we
460 : // could throw an error since trying to change an elemset code you
461 : // don't have could indicate there's a problem...
462 284 : if (it == _elemset_codes.end())
463 0 : return;
464 :
465 : // Make copy of the set of elemset ids. We are not changing these,
466 : // only updating the elemset code it corresponds to.
467 284 : elemset_type id_set_copy = *(it->second);
468 :
469 : // Look up the corresponding entry in the inverse map. Note: we want
470 : // the iterator because we are going to remove it.
471 8 : auto inverse_it = _elemset_codes_inverse_map.find(id_set_copy);
472 284 : libmesh_error_msg_if(inverse_it == _elemset_codes_inverse_map.end(),
473 : "Expected _elemset_codes_inverse_map entry for elemset code " << old_code);
474 :
475 : // Erase entry from inverse map
476 276 : _elemset_codes_inverse_map.erase(inverse_it);
477 :
478 : // Erase entry from forward map
479 276 : _elemset_codes.erase(it);
480 :
481 : // Add new code with original set of ids.
482 560 : this->add_elemset_code(new_code, id_set_copy);
483 :
484 : // We can't update any actual elemset codes if there is no extra integer defined for it.
485 284 : if (!this->has_elem_integer("elemset_code"))
486 0 : return;
487 :
488 : // Get index of elemset_code extra integer
489 284 : unsigned int elemset_index = this->get_elem_integer_index("elemset_code");
490 :
491 : // Loop over all elems and update code
492 10220 : for (auto & elem : this->element_ptr_range())
493 : {
494 : dof_id_type elemset_code =
495 4970 : elem->get_extra_integer(elemset_index);
496 :
497 4970 : if (elemset_code == old_code)
498 2415 : elem->set_extra_integer(elemset_index, new_code);
499 268 : }
500 : }
501 :
502 284 : void MeshBase::change_elemset_id(elemset_id_type old_id, elemset_id_type new_id)
503 : {
504 : // Early return if we don't have old_id
505 8 : if (!_all_elemset_ids.count(old_id))
506 0 : return;
507 :
508 : // Throw an error if the new_id is already used
509 8 : libmesh_error_msg_if(_all_elemset_ids.count(new_id),
510 : "Cannot change elemset id " << old_id <<
511 : " to " << new_id << ", " << new_id << " already exists.");
512 :
513 : // We will build up a new version of the inverse map so we can iterate over
514 : // the current one without invalidating anything.
515 16 : std::map<MeshBase::elemset_type, dof_id_type> new_elemset_codes_inverse_map;
516 852 : for (const auto & [id_set, elemset_code] : _elemset_codes_inverse_map)
517 : {
518 32 : auto id_set_copy = id_set;
519 16 : if (id_set_copy.count(old_id))
520 : {
521 : // Remove old_id, insert new_id
522 8 : id_set_copy.erase(old_id);
523 276 : id_set_copy.insert(new_id);
524 : }
525 :
526 : // Store in new version of map
527 552 : new_elemset_codes_inverse_map.emplace(id_set_copy, elemset_code);
528 : }
529 :
530 : // Swap existing map with newly-built one
531 8 : _elemset_codes_inverse_map.swap(new_elemset_codes_inverse_map);
532 :
533 : // Reconstruct _elemset_codes map
534 8 : _elemset_codes.clear();
535 852 : for (const auto & [id_set, elemset_code] : _elemset_codes_inverse_map)
536 568 : _elemset_codes.emplace(elemset_code, &id_set);
537 :
538 : // Update _all_elemset_ids
539 8 : _all_elemset_ids.erase(old_id);
540 276 : _all_elemset_ids.insert(new_id);
541 : }
542 :
543 259771 : unsigned int MeshBase::spatial_dimension () const
544 : {
545 268551 : return cast_int<unsigned int>(_spatial_dimension);
546 : }
547 :
548 :
549 :
550 353402 : void MeshBase::set_spatial_dimension(unsigned char d)
551 : {
552 : // The user can set the _spatial_dimension however they wish,
553 : // libMesh will only *increase* the spatial dimension, however,
554 : // never decrease it.
555 353402 : _spatial_dimension = d;
556 353402 : }
557 :
558 :
559 :
560 4077 : unsigned int MeshBase::add_elem_integer(std::string name,
561 : bool allocate_data,
562 : dof_id_type default_value)
563 : {
564 5273 : for (auto i : index_range(_elem_integer_names))
565 1416 : if (_elem_integer_names[i] == name)
566 : {
567 8 : libmesh_assert_less(i, _elem_integer_default_values.size());
568 284 : _elem_integer_default_values[i] = default_value;
569 284 : return i;
570 : }
571 :
572 118 : libmesh_assert_equal_to(_elem_integer_names.size(),
573 : _elem_integer_default_values.size());
574 3793 : _elem_integer_names.push_back(std::move(name));
575 3793 : _elem_integer_default_values.push_back(default_value);
576 3793 : if (allocate_data)
577 3205 : this->size_elem_extra_integers();
578 3911 : return _elem_integer_names.size()-1;
579 : }
580 :
581 :
582 :
583 27194 : std::vector<unsigned int> MeshBase::add_elem_integers(const std::vector<std::string> & names,
584 : bool allocate_data,
585 : const std::vector<dof_id_type> * default_values)
586 : {
587 898 : libmesh_assert(!default_values || default_values->size() == names.size());
588 898 : libmesh_assert_equal_to(_elem_integer_names.size(), _elem_integer_default_values.size());
589 :
590 1796 : std::unordered_map<std::string, std::size_t> name_indices;
591 28046 : for (auto i : index_range(_elem_integer_names))
592 852 : name_indices[_elem_integer_names[i]] = i;
593 :
594 28076 : std::vector<unsigned int> returnval(names.size());
595 :
596 898 : bool added_an_integer = false;
597 29699 : for (auto i : index_range(names))
598 : {
599 156 : const std::string & name = names[i];
600 2505 : if (const auto it = name_indices.find(name);
601 78 : it != name_indices.end())
602 : {
603 426 : returnval[i] = it->second;
604 438 : _elem_integer_default_values[it->second] =
605 426 : default_values ? (*default_values)[i] : DofObject::invalid_id;
606 : }
607 : else
608 : {
609 2145 : returnval[i] = _elem_integer_names.size();
610 2079 : name_indices[name] = returnval[i];
611 2079 : _elem_integer_names.push_back(name);
612 : _elem_integer_default_values.push_back
613 2883 : (default_values ? (*default_values)[i] : DofObject::invalid_id);
614 66 : added_an_integer = true;
615 : }
616 : }
617 :
618 27194 : if (allocate_data && added_an_integer)
619 1065 : this->size_elem_extra_integers();
620 :
621 28092 : return returnval;
622 : }
623 :
624 :
625 :
626 3278 : unsigned int MeshBase::get_elem_integer_index(std::string_view name) const
627 : {
628 4130 : for (auto i : index_range(_elem_integer_names))
629 4037 : if (_elem_integer_names[i] == name)
630 3278 : return i;
631 :
632 0 : libmesh_error_msg("Unknown elem integer " << name);
633 : return libMesh::invalid_uint;
634 : }
635 :
636 :
637 :
638 8219 : bool MeshBase::has_elem_integer(std::string_view name) const
639 : {
640 9213 : for (auto & entry : _elem_integer_names)
641 2367 : if (entry == name)
642 40 : return true;
643 :
644 622 : return false;
645 : }
646 :
647 :
648 :
649 14636 : unsigned int MeshBase::add_node_integer(std::string name,
650 : bool allocate_data,
651 : dof_id_type default_value)
652 : {
653 24899 : for (auto i : index_range(_node_integer_names))
654 10118 : if (_node_integer_names[i] == name)
655 : {
656 16 : libmesh_assert_less(i, _node_integer_default_values.size());
657 568 : _node_integer_default_values[i] = default_value;
658 568 : return i;
659 : }
660 :
661 582 : libmesh_assert_equal_to(_node_integer_names.size(),
662 : _node_integer_default_values.size());
663 14068 : _node_integer_names.push_back(std::move(name));
664 14068 : _node_integer_default_values.push_back(default_value);
665 14068 : if (allocate_data)
666 7120 : this->size_node_extra_integers();
667 14650 : return _node_integer_names.size()-1;
668 : }
669 :
670 :
671 :
672 27052 : std::vector<unsigned int> MeshBase::add_node_integers(const std::vector<std::string> & names,
673 : bool allocate_data,
674 : const std::vector<dof_id_type> * default_values)
675 : {
676 894 : libmesh_assert(!default_values || default_values->size() == names.size());
677 894 : libmesh_assert_equal_to(_node_integer_names.size(), _node_integer_default_values.size());
678 :
679 1788 : std::unordered_map<std::string, std::size_t> name_indices;
680 27620 : for (auto i : index_range(_node_integer_names))
681 568 : name_indices[_node_integer_names[i]] = i;
682 :
683 27930 : std::vector<unsigned int> returnval(names.size());
684 :
685 894 : bool added_an_integer = false;
686 29684 : for (auto i : index_range(names))
687 : {
688 192 : const std::string & name = names[i];
689 2632 : if (const auto it = name_indices.find(name);
690 96 : it != name_indices.end())
691 : {
692 0 : returnval[i] = it->second;
693 0 : _node_integer_default_values[it->second] =
694 0 : default_values ? (*default_values)[i] : DofObject::invalid_id;
695 : }
696 : else
697 : {
698 2728 : returnval[i] = _node_integer_names.size();
699 2632 : name_indices[name] = returnval[i];
700 2632 : _node_integer_names.push_back(name);
701 : _node_integer_default_values.push_back
702 3724 : (default_values ? (*default_values)[i] : DofObject::invalid_id);
703 96 : added_an_integer = true;
704 : }
705 : }
706 :
707 27052 : if (allocate_data && added_an_integer)
708 1010 : this->size_node_extra_integers();
709 :
710 27946 : return returnval;
711 : }
712 :
713 :
714 :
715 710 : unsigned int MeshBase::get_node_integer_index(std::string_view name) const
716 : {
717 1988 : for (auto i : index_range(_node_integer_names))
718 1968 : if (_node_integer_names[i] == name)
719 710 : return i;
720 :
721 0 : libmesh_error_msg("Unknown node integer " << name);
722 : return libMesh::invalid_uint;
723 : }
724 :
725 :
726 :
727 1136 : bool MeshBase::has_node_integer(std::string_view name) const
728 : {
729 3134 : for (auto & entry : _node_integer_names)
730 3134 : if (entry == name)
731 32 : return true;
732 :
733 0 : return false;
734 : }
735 :
736 :
737 :
738 58231 : void MeshBase::remove_orphaned_nodes ()
739 : {
740 3700 : LOG_SCOPE("remove_orphaned_nodes()", "MeshBase");
741 :
742 : // Will hold the set of nodes that are currently connected to elements
743 3700 : std::unordered_set<Node *> connected_nodes;
744 :
745 : // Loop over the elements. Find which nodes are connected to at
746 : // least one of them.
747 22985286 : for (const auto & element : this->element_ptr_range())
748 70791017 : for (auto & n : element->node_ref_range())
749 58533043 : connected_nodes.insert(&n);
750 :
751 18276420 : for (const auto & node : this->node_ptr_range())
752 1044710 : if (!connected_nodes.count(node))
753 135899 : this->delete_node(node);
754 58231 : }
755 :
756 :
757 :
758 : #ifdef LIBMESH_ENABLE_DEPRECATED
759 0 : void MeshBase::prepare_for_use (const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
760 : {
761 : libmesh_deprecated();
762 :
763 : // We only respect the users wish if they tell us to skip renumbering. If they tell us not to
764 : // skip renumbering but someone previously called allow_renumbering(false), then the latter takes
765 : // precedence
766 0 : if (skip_renumber_nodes_and_elements)
767 0 : this->allow_renumbering(false);
768 :
769 : // We always accept the user's value for skip_find_neighbors, in contrast to skip_renumber
770 0 : const bool old_allow_find_neighbors = this->allow_find_neighbors();
771 0 : this->allow_find_neighbors(!skip_find_neighbors);
772 :
773 0 : this->prepare_for_use();
774 :
775 0 : this->allow_find_neighbors(old_allow_find_neighbors);
776 0 : }
777 :
778 0 : void MeshBase::prepare_for_use (const bool skip_renumber_nodes_and_elements)
779 : {
780 : libmesh_deprecated();
781 :
782 : // We only respect the users wish if they tell us to skip renumbering. If they tell us not to
783 : // skip renumbering but someone previously called allow_renumbering(false), then the latter takes
784 : // precedence
785 0 : if (skip_renumber_nodes_and_elements)
786 0 : this->allow_renumbering(false);
787 :
788 0 : this->prepare_for_use();
789 0 : }
790 : #endif // LIBMESH_ENABLE_DEPRECATED
791 :
792 :
793 :
794 412566 : void MeshBase::prepare_for_use ()
795 : {
796 25172 : LOG_SCOPE("prepare_for_use()", "MeshBase");
797 :
798 12586 : parallel_object_only();
799 :
800 12586 : libmesh_assert(this->comm().verify(this->is_serial()));
801 :
802 : // If we don't go into this method with valid constraint rows, we're
803 : // only going to be able to make that worse.
804 : #ifdef DEBUG
805 12586 : MeshTools::libmesh_assert_valid_constraint_rows(*this);
806 : #endif
807 :
808 : // A distributed mesh may have processors with no elements (or
809 : // processors with no elements of higher dimension, if we ever
810 : // support mixed-dimension meshes), but we want consistent
811 : // mesh_dimension anyways.
812 : //
813 : // cache_elem_data() should get the elem_dimensions() and
814 : // mesh_dimension() correct later, and we don't need it earlier.
815 :
816 :
817 : // Renumber the nodes and elements so that they in contiguous
818 : // blocks. By default, _skip_renumber_nodes_and_elements is false.
819 : //
820 : // Instances where you if prepare_for_use() should not renumber the nodes
821 : // and elements include reading in e.g. an xda/r or gmv file. In
822 : // this case, the ordering of the nodes may depend on an accompanying
823 : // solution, and the node ordering cannot be changed.
824 :
825 :
826 : // Mesh modification operations might not leave us with consistent
827 : // id counts, or might leave us with orphaned nodes we're no longer
828 : // using, but our partitioner might need that consistency and/or
829 : // might be confused by orphaned nodes.
830 412566 : if (!_skip_renumber_nodes_and_elements)
831 354335 : this->renumber_nodes_and_elements();
832 : else
833 : {
834 58231 : this->remove_orphaned_nodes();
835 58231 : this->update_parallel_id_counts();
836 : }
837 :
838 : // Let all the elements find their neighbors
839 412566 : if (!_skip_find_neighbors)
840 380826 : this->find_neighbors();
841 :
842 : // The user may have set boundary conditions. We require that the
843 : // boundary conditions were set consistently. Because we examine
844 : // neighbors when evaluating non-raw boundary condition IDs, this
845 : // assert is only valid when our neighbor links are in place.
846 : #ifdef DEBUG
847 12586 : MeshTools::libmesh_assert_valid_boundary_ids(*this);
848 : #endif
849 :
850 : // Search the mesh for all the dimensions of the elements
851 : // and cache them.
852 412566 : this->cache_elem_data();
853 :
854 : // Search the mesh for elements that have a neighboring element
855 : // of dim+1 and set that element as the interior parent
856 412566 : this->detect_interior_parents();
857 :
858 : // Fix up node unique ids in case mesh generation code didn't take
859 : // exceptional care to do so.
860 : // MeshCommunication().make_node_unique_ids_parallel_consistent(*this);
861 :
862 : // We're going to still require that mesh generation code gets
863 : // element unique ids consistent.
864 : #if defined(DEBUG) && defined(LIBMESH_ENABLE_UNIQUE_ID)
865 12586 : MeshTools::libmesh_assert_valid_unique_ids(*this);
866 : #endif
867 :
868 : // Reset our PointLocator. Any old locator is invalidated any time
869 : // the elements in the underlying elements in the mesh have changed,
870 : // so we clear it here.
871 412566 : this->clear_point_locator();
872 :
873 : // Allow our GhostingFunctor objects to reinit if necessary.
874 : // Do this before partitioning and redistributing, and before
875 : // deleting remote elements.
876 412566 : this->reinit_ghosting_functors();
877 :
878 : // Partition the mesh unless *all* partitioning is to be skipped.
879 : // If only noncritical partitioning is to be skipped, the
880 : // partition() call will still check for orphaned nodes.
881 412566 : if (!skip_partitioning())
882 11768 : this->partition();
883 :
884 : // If we're using DistributedMesh, we'll probably want it
885 : // parallelized.
886 412566 : if (this->_allow_remote_element_removal)
887 380628 : this->delete_remote_elements();
888 :
889 : // Much of our boundary info may have been for now-remote parts of the mesh,
890 : // in which case we don't want to keep local copies of data meant to be
891 : // local. On the other hand we may have deleted, or the user may have added in
892 : // a distributed fashion, boundary data that is meant to be global. So we
893 : // handle both of those scenarios here
894 412566 : this->get_boundary_info().regenerate_id_sets();
895 :
896 412566 : if (!_skip_renumber_nodes_and_elements)
897 354335 : this->renumber_nodes_and_elements();
898 :
899 : // The mesh is now prepared for use.
900 412566 : _is_prepared = true;
901 :
902 : #ifdef DEBUG
903 12586 : MeshTools::libmesh_assert_valid_boundary_ids(*this);
904 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
905 12586 : MeshTools::libmesh_assert_valid_unique_ids(*this);
906 : #endif
907 : #endif
908 412566 : }
909 :
910 : void
911 412566 : MeshBase::reinit_ghosting_functors()
912 : {
913 951666 : for (auto & gf : _ghosting_functors)
914 : {
915 16738 : libmesh_assert(gf);
916 539100 : gf->mesh_reinit();
917 : }
918 412566 : }
919 :
920 930658 : void MeshBase::clear ()
921 : {
922 : // Reset the number of partitions
923 930658 : _n_parts = 1;
924 :
925 : // Reset the _is_prepared flag
926 930658 : _is_prepared = false;
927 :
928 : // Clear boundary information
929 930658 : if (boundary_info)
930 929380 : boundary_info->clear();
931 :
932 : // Clear cached element data
933 28060 : _elem_dims.clear();
934 28060 : _elem_default_orders.clear();
935 930658 : _supported_nodal_order = MAXIMUM;
936 :
937 28060 : _elemset_codes.clear();
938 28060 : _elemset_codes_inverse_map.clear();
939 :
940 28060 : _constraint_rows.clear();
941 :
942 : // Clear our point locator.
943 930658 : this->clear_point_locator();
944 930658 : }
945 :
946 :
947 :
948 983713 : void MeshBase::add_ghosting_functor(GhostingFunctor & ghosting_functor)
949 : {
950 : // We used to implicitly support duplicate inserts to std::set
951 : #ifdef LIBMESH_ENABLE_DEPRECATED
952 : _ghosting_functors.erase
953 927277 : (std::remove(_ghosting_functors.begin(),
954 : _ghosting_functors.end(),
955 983713 : &ghosting_functor),
956 84654 : _ghosting_functors.end());
957 : #endif
958 :
959 : // We shouldn't have two copies of the same functor
960 28218 : libmesh_assert(std::find(_ghosting_functors.begin(),
961 : _ghosting_functors.end(),
962 : &ghosting_functor) ==
963 : _ghosting_functors.end());
964 :
965 983713 : _ghosting_functors.push_back(&ghosting_functor);
966 983713 : }
967 :
968 :
969 :
970 953171 : void MeshBase::remove_ghosting_functor(GhostingFunctor & ghosting_functor)
971 : {
972 898759 : auto raw_it = std::find(_ghosting_functors.begin(),
973 980377 : _ghosting_functors.end(), &ghosting_functor);
974 :
975 : // The DofMap has a "to_mesh" parameter that tells it to avoid
976 : // registering a new functor with the mesh, but it doesn't keep
977 : // track of which functors weren't added, so we'll support "remove a
978 : // functor that isn't there" just like we did with set::erase
979 : // before.
980 953171 : if (raw_it != _ghosting_functors.end())
981 952248 : _ghosting_functors.erase(raw_it);
982 :
983 : // We shouldn't have had two copies of the same functor
984 27206 : libmesh_assert(std::find(_ghosting_functors.begin(),
985 : _ghosting_functors.end(),
986 : &ghosting_functor) ==
987 : _ghosting_functors.end());
988 :
989 953171 : if (const auto it = _shared_functors.find(&ghosting_functor);
990 27206 : it != _shared_functors.end())
991 69 : _shared_functors.erase(it);
992 953171 : }
993 :
994 :
995 :
996 23798 : void MeshBase::subdomain_ids (std::set<subdomain_id_type> & ids, const bool global /* = true */) const
997 : {
998 : // This requires an inspection on every processor
999 1202 : if (global)
1000 1202 : parallel_object_only();
1001 :
1002 1202 : ids.clear();
1003 :
1004 3983786 : for (const auto & elem : this->active_local_element_ptr_range())
1005 2196172 : ids.insert(elem->subdomain_id());
1006 :
1007 23798 : if (global)
1008 : {
1009 : // Only include the unpartitioned elements if the user requests the global IDs.
1010 : // In the case of the local subdomain IDs, it doesn't make sense to include the
1011 : // unpartitioned elements because said elements do not have a sense of locality.
1012 49868 : for (const auto & elem : this->active_unpartitioned_element_ptr_range())
1013 24868 : ids.insert(elem->subdomain_id());
1014 :
1015 : // Some subdomains may only live on other processors
1016 23798 : this->comm().set_union(ids);
1017 : }
1018 23798 : }
1019 :
1020 :
1021 :
1022 567528 : void MeshBase::redistribute()
1023 : {
1024 : // We now have all elements and nodes redistributed; our ghosting
1025 : // functors should be ready to redistribute and/or recompute any
1026 : // cached data they use too.
1027 570304 : for (auto & gf : as_range(this->ghosting_functors_begin(),
1028 1284517 : this->ghosting_functors_end()))
1029 690793 : gf->redistribute();
1030 567528 : }
1031 :
1032 :
1033 :
1034 16973 : subdomain_id_type MeshBase::n_subdomains() const
1035 : {
1036 : // This requires an inspection on every processor
1037 500 : parallel_object_only();
1038 :
1039 1000 : std::set<subdomain_id_type> ids;
1040 :
1041 16973 : this->subdomain_ids (ids);
1042 :
1043 17473 : return cast_int<subdomain_id_type>(ids.size());
1044 : }
1045 :
1046 :
1047 :
1048 0 : subdomain_id_type MeshBase::n_local_subdomains() const
1049 : {
1050 0 : std::set<subdomain_id_type> ids;
1051 :
1052 0 : this->subdomain_ids (ids, /* global = */ false);
1053 :
1054 0 : return cast_int<subdomain_id_type>(ids.size());
1055 : }
1056 :
1057 :
1058 :
1059 :
1060 3204007 : dof_id_type MeshBase::n_nodes_on_proc (const processor_id_type proc_id) const
1061 : {
1062 : // We're either counting a processor's nodes or unpartitioned
1063 : // nodes
1064 5294 : libmesh_assert (proc_id < this->n_processors() ||
1065 : proc_id == DofObject::invalid_processor_id);
1066 :
1067 6402720 : return static_cast<dof_id_type>(std::distance (this->pid_nodes_begin(proc_id),
1068 9606727 : this->pid_nodes_end (proc_id)));
1069 : }
1070 :
1071 :
1072 :
1073 3592607 : dof_id_type MeshBase::n_elem_on_proc (const processor_id_type proc_id) const
1074 : {
1075 : // We're either counting a processor's elements or unpartitioned
1076 : // elements
1077 17072 : libmesh_assert (proc_id < this->n_processors() ||
1078 : proc_id == DofObject::invalid_processor_id);
1079 :
1080 7168142 : return static_cast<dof_id_type>(std::distance (this->pid_elements_begin(proc_id),
1081 10760749 : this->pid_elements_end (proc_id)));
1082 : }
1083 :
1084 :
1085 :
1086 1151467 : dof_id_type MeshBase::n_active_elem_on_proc (const processor_id_type proc_id) const
1087 : {
1088 2004 : libmesh_assert_less (proc_id, this->n_processors());
1089 2300930 : return static_cast<dof_id_type>(std::distance (this->active_pid_elements_begin(proc_id),
1090 3452397 : this->active_pid_elements_end (proc_id)));
1091 : }
1092 :
1093 :
1094 :
1095 0 : dof_id_type MeshBase::n_sub_elem () const
1096 : {
1097 0 : dof_id_type ne=0;
1098 :
1099 0 : for (const auto & elem : this->element_ptr_range())
1100 0 : ne += elem->n_sub_elem();
1101 :
1102 0 : return ne;
1103 : }
1104 :
1105 :
1106 :
1107 32834 : dof_id_type MeshBase::n_active_sub_elem () const
1108 : {
1109 1591 : dof_id_type ne=0;
1110 :
1111 26592389 : for (const auto & elem : this->active_element_ptr_range())
1112 26557964 : ne += elem->n_sub_elem();
1113 :
1114 32834 : return ne;
1115 : }
1116 :
1117 :
1118 :
1119 14844 : std::string MeshBase::get_info(const unsigned int verbosity /* = 0 */, const bool global /* = true */) const
1120 : {
1121 15720 : std::ostringstream oss;
1122 :
1123 14844 : oss << " Mesh Information:" << '\n';
1124 :
1125 14844 : if (!_elem_dims.empty())
1126 : {
1127 14844 : oss << " elem_dimensions()={";
1128 14844 : std::copy(_elem_dims.begin(),
1129 438 : --_elem_dims.end(), // --end() is valid if the set is non-empty
1130 14406 : std::ostream_iterator<unsigned int>(oss, ", "));
1131 14844 : oss << cast_int<unsigned int>(*_elem_dims.rbegin());
1132 14844 : oss << "}\n";
1133 : }
1134 :
1135 14844 : if (!_elem_default_orders.empty())
1136 : {
1137 14844 : oss << " elem_default_orders()={";
1138 14844 : std::transform(_elem_default_orders.begin(),
1139 438 : --_elem_default_orders.end(),
1140 14844 : std::ostream_iterator<std::string>(oss, ", "),
1141 8 : [](Order o)
1142 276 : { return Utility::enum_to_string<Order>(o); });
1143 14844 : oss << Utility::enum_to_string<Order>(*_elem_default_orders.rbegin());
1144 14844 : oss << "}\n";
1145 : }
1146 :
1147 14844 : oss << " supported_nodal_order()=" << this->supported_nodal_order() << '\n'
1148 15720 : << " spatial_dimension()=" << this->spatial_dimension() << '\n'
1149 15282 : << " n_nodes()=" << this->n_nodes() << '\n'
1150 14844 : << " n_local_nodes()=" << this->n_local_nodes() << '\n'
1151 15282 : << " n_elem()=" << this->n_elem() << '\n'
1152 29250 : << " n_local_elem()=" << this->n_local_elem() << '\n';
1153 : #ifdef LIBMESH_ENABLE_AMR
1154 28812 : oss << " n_active_elem()=" << this->n_active_elem() << '\n';
1155 : #endif
1156 14844 : if (global)
1157 15282 : oss << " n_subdomains()=" << static_cast<std::size_t>(this->n_subdomains()) << '\n';
1158 : else
1159 0 : oss << " n_local_subdomains()= " << static_cast<std::size_t>(this->n_local_subdomains()) << '\n';
1160 15282 : oss << " n_elemsets()=" << static_cast<std::size_t>(this->n_elemsets()) << '\n';
1161 14844 : if (!_elemset_codes.empty())
1162 0 : oss << " n_elemset_codes=" << _elemset_codes.size() << '\n';
1163 14844 : oss << " n_partitions()=" << static_cast<std::size_t>(this->n_partitions()) << '\n'
1164 15720 : << " n_processors()=" << static_cast<std::size_t>(this->n_processors()) << '\n'
1165 15282 : << " n_threads()=" << static_cast<std::size_t>(libMesh::n_threads()) << '\n'
1166 15720 : << " processor_id()=" << static_cast<std::size_t>(this->processor_id()) << '\n'
1167 876 : << " is_prepared()=" << (this->is_prepared() ? "true" : "false") << '\n'
1168 28773 : << " is_replicated()=" << (this->is_replicated() ? "true" : "false") << '\n';
1169 :
1170 14844 : if (verbosity > 0)
1171 : {
1172 284 : if (global)
1173 : {
1174 8 : libmesh_parallel_only(this->comm());
1175 292 : if (this->processor_id() != 0)
1176 236 : oss << "\n Detailed global get_info() (verbosity > 0) is reduced and output to only rank 0.";
1177 : }
1178 :
1179 : // Helper for printing element types
1180 672 : const auto elem_type_helper = [](const std::set<int> & elem_types) {
1181 784 : std::stringstream ss;
1182 1344 : for (auto it = elem_types.begin(); it != elem_types.end();)
1183 : {
1184 1288 : ss << Utility::enum_to_string((ElemType)*it);
1185 672 : if (++it != elem_types.end())
1186 0 : ss << ", ";
1187 : }
1188 728 : return ss.str();
1189 560 : };
1190 :
1191 : // Helper for whether or not the given DofObject is to be included. If we're doing
1192 : // a global reduction, we also count unpartitioned objects on rank 0.
1193 871156 : const auto include_object = [this, &global](const DofObject & dof_object) {
1194 978188 : return this->processor_id() == dof_object.processor_id() ||
1195 250756 : (global &&
1196 42324 : this->processor_id() == 0 &&
1197 846328 : dof_object.processor_id() == DofObject::invalid_processor_id);
1198 276 : };
1199 :
1200 8 : Real volume = 0;
1201 :
1202 : // Add bounding box information
1203 284 : const auto bbox = global ? MeshTools::create_bounding_box(*this) : MeshTools::create_local_bounding_box(*this);
1204 284 : if (!global || this->processor_id() == 0)
1205 : oss << "\n " << (global ? "" : "Local ") << "Mesh Bounding Box:\n"
1206 48 : << " Minimum: " << bbox.min() << "\n"
1207 44 : << " Maximum: " << bbox.max() << "\n"
1208 88 : << " Delta: " << (bbox.max() - bbox.min()) << "\n";
1209 :
1210 : // Obtain the global or local element types
1211 16 : std::set<int> elem_types;
1212 49712 : for (const Elem * elem : this->active_local_element_ptr_range())
1213 49420 : elem_types.insert(elem->type());
1214 284 : if (global)
1215 : {
1216 : // Pick up unpartitioned elems on rank 0
1217 292 : if (this->processor_id() == 0)
1218 92 : for (const Elem * elem : this->active_unpartitioned_element_ptr_range())
1219 40 : elem_types.insert(elem->type());
1220 :
1221 284 : this->comm().set_union(elem_types);
1222 : }
1223 :
1224 : // Add element types
1225 284 : if (!global || this->processor_id() == 0)
1226 : oss << "\n " << (global ? "" : "Local ") << "Mesh Element Type(s):\n "
1227 140 : << elem_type_helper(elem_types) << "\n";
1228 :
1229 : // Reduce the nodeset ids
1230 16 : auto nodeset_ids = this->get_boundary_info().get_node_boundary_ids();
1231 284 : if (global)
1232 284 : this->comm().set_union(nodeset_ids);
1233 :
1234 : // Accumulate local information for each nodeset
1235 1608 : struct NodesetInfo
1236 : {
1237 : std::size_t num_nodes = 0;
1238 : BoundingBox bbox;
1239 : };
1240 16 : std::map<boundary_id_type, NodesetInfo> nodeset_info_map;
1241 91376 : for (const auto & [node, id] : this->get_boundary_info().get_nodeset_map())
1242 : {
1243 91092 : if (!include_object(*node))
1244 38364 : continue;
1245 :
1246 48672 : NodesetInfo & info = nodeset_info_map[id];
1247 :
1248 48672 : ++info.num_nodes;
1249 :
1250 48672 : if (verbosity > 1)
1251 48672 : info.bbox.union_with(*node);
1252 : }
1253 :
1254 : // Add nodeset info
1255 284 : if (!global || this->processor_id() == 0)
1256 : {
1257 48 : oss << "\n " << (global ? "" : "Local ") << "Mesh Nodesets:\n";
1258 48 : if (nodeset_ids.empty())
1259 0 : oss << " None\n";
1260 : }
1261 :
1262 8 : const auto & nodeset_name_map = this->get_boundary_info().get_nodeset_name_map();
1263 1988 : for (const auto id : nodeset_ids)
1264 : {
1265 1704 : NodesetInfo & info = nodeset_info_map[id];
1266 :
1267 : // Reduce the local information for this nodeset if required
1268 1704 : if (global)
1269 : {
1270 1704 : this->comm().sum(info.num_nodes);
1271 1704 : if (verbosity > 1)
1272 : {
1273 1752 : this->comm().min(info.bbox.min());
1274 1752 : this->comm().max(info.bbox.max());
1275 : }
1276 : }
1277 :
1278 1704 : const bool has_name = nodeset_name_map.count(id) && nodeset_name_map.at(id).size();
1279 1800 : const std::string name = has_name ? nodeset_name_map.at(id) : "";
1280 96 : if (global)
1281 48 : libmesh_assert(this->comm().verify(name));
1282 :
1283 1704 : if (global ? this->processor_id() == 0 : info.num_nodes > 0)
1284 : {
1285 288 : oss << " Nodeset " << id;
1286 288 : if (has_name)
1287 528 : oss << " (" << name << ")";
1288 312 : oss << ", " << info.num_nodes << " " << (global ? "" : "local ") << "nodes\n";
1289 :
1290 288 : if (verbosity > 1)
1291 : {
1292 288 : oss << " " << (global ? "Bounding" : "Local bounding") << " box minimum: "
1293 288 : << info.bbox.min() << "\n"
1294 288 : << " " << (global ? "Bounding" : "Local bounding") << " box maximum: "
1295 288 : << info.bbox.max() << "\n"
1296 288 : << " " << (global ? "Bounding" : "Local bounding") << " box delta: "
1297 288 : << (info.bbox.max() - info.bbox.min()) << "\n";
1298 : }
1299 : }
1300 : }
1301 :
1302 : // Reduce the sideset ids
1303 16 : auto sideset_ids = this->get_boundary_info().get_side_boundary_ids();
1304 284 : if (global)
1305 284 : this->comm().set_union(sideset_ids);
1306 :
1307 : // Accumulate local information for each sideset
1308 : struct SidesetInfo
1309 : {
1310 : std::size_t num_sides = 0;
1311 : Real volume = 0;
1312 : std::set<int> side_elem_types;
1313 : std::set<int> elem_types;
1314 : std::set<dof_id_type> elem_ids;
1315 : std::set<dof_id_type> node_ids;
1316 : BoundingBox bbox;
1317 : };
1318 16 : ElemSideBuilder side_builder;
1319 16 : std::map<boundary_id_type, SidesetInfo> sideset_info_map;
1320 70396 : for (const auto & pair : this->get_boundary_info().get_sideset_map())
1321 : {
1322 70112 : const Elem * elem = pair.first;
1323 63456 : if (!include_object(*elem))
1324 30176 : continue;
1325 :
1326 39936 : const auto id = pair.second.second;
1327 39936 : SidesetInfo & info = sideset_info_map[id];
1328 :
1329 39936 : const auto s = pair.second.first;
1330 39936 : const Elem & side = side_builder(*elem, s);
1331 :
1332 39936 : ++info.num_sides;
1333 39936 : info.side_elem_types.insert(side.type());
1334 39936 : info.elem_types.insert(elem->type());
1335 39936 : info.elem_ids.insert(elem->id());
1336 :
1337 199680 : for (const Node & node : side.node_ref_range())
1338 146432 : if (include_object(node))
1339 147552 : info.node_ids.insert(node.id());
1340 :
1341 39936 : if (verbosity > 1)
1342 : {
1343 39936 : info.volume += side.volume();
1344 39936 : info.bbox.union_with(side.loose_bounding_box());
1345 : }
1346 : }
1347 :
1348 : // Add sideset info
1349 284 : if (!global || this->processor_id() == 0)
1350 : {
1351 48 : oss << "\n " << (global ? "" : "Local ") << "Mesh Sidesets:\n";
1352 48 : if (sideset_ids.empty())
1353 0 : oss << " None\n";
1354 : }
1355 8 : const auto & sideset_name_map = this->get_boundary_info().get_sideset_name_map();
1356 1988 : for (const auto id : sideset_ids)
1357 : {
1358 1704 : SidesetInfo & info = sideset_info_map[id];
1359 :
1360 1704 : auto num_elems = info.elem_ids.size();
1361 1704 : auto num_nodes = info.node_ids.size();
1362 :
1363 : // Reduce the local information for this sideset if required
1364 1704 : if (global)
1365 : {
1366 1704 : this->comm().sum(info.num_sides);
1367 1704 : this->comm().set_union(info.side_elem_types, 0);
1368 1704 : this->comm().sum(num_elems);
1369 1704 : this->comm().set_union(info.elem_types, 0);
1370 1704 : this->comm().sum(num_nodes);
1371 1704 : if (verbosity > 1)
1372 : {
1373 1704 : this->comm().sum(info.volume);
1374 1752 : this->comm().min(info.bbox.min());
1375 1752 : this->comm().max(info.bbox.max());
1376 : }
1377 : }
1378 :
1379 1704 : const bool has_name = sideset_name_map.count(id) && sideset_name_map.at(id).size();
1380 1800 : const std::string name = has_name ? sideset_name_map.at(id) : "";
1381 96 : if (global)
1382 48 : libmesh_assert(this->comm().verify(name));
1383 :
1384 1704 : if (global ? this->processor_id() == 0 : info.num_sides > 0)
1385 : {
1386 288 : oss << " Sideset " << id;
1387 288 : if (has_name)
1388 528 : oss << " (" << name << ")";
1389 552 : oss << ", " << info.num_sides << " sides (" << elem_type_helper(info.side_elem_types) << ")"
1390 1056 : << ", " << num_elems << " " << (global ? "" : "local ") << "elems (" << elem_type_helper(info.elem_types) << ")"
1391 600 : << ", " << num_nodes << " " << (global ? "" : "local ") << "nodes\n";
1392 :
1393 288 : if (verbosity > 1)
1394 : {
1395 312 : oss << " " << (global ? "Side" : "Local side") << " volume: " << info.volume << "\n"
1396 288 : << " " << (global ? "Bounding" : "Local bounding") << " box minimum: "
1397 288 : << info.bbox.min() << "\n"
1398 288 : << " " << (global ? "Bounding" : "Local bounding") << " box maximum: "
1399 288 : << info.bbox.max() << "\n"
1400 288 : << " " << (global ? "Bounding" : "Local bounding") << " box delta: "
1401 288 : << (info.bbox.max() - info.bbox.min()) << "\n";
1402 : }
1403 : }
1404 : }
1405 :
1406 : // Reduce the edgeset ids
1407 16 : auto edgeset_ids = this->get_boundary_info().get_edge_boundary_ids();
1408 284 : if (global)
1409 284 : this->comm().set_union(edgeset_ids);
1410 :
1411 : // Accumulate local information for each edgeset
1412 0 : struct EdgesetInfo
1413 : {
1414 : std::size_t num_edges = 0;
1415 : std::set<int> edge_elem_types;
1416 : BoundingBox bbox;
1417 : };
1418 16 : std::map<boundary_id_type, EdgesetInfo> edgeset_info_map;
1419 284 : std::unique_ptr<const Elem> edge;
1420 :
1421 284 : for (const auto & pair : this->get_boundary_info().get_edgeset_map())
1422 : {
1423 0 : const Elem * elem = pair.first;
1424 0 : if (!include_object(*elem))
1425 0 : continue;
1426 :
1427 0 : const auto id = pair.second.second;
1428 0 : EdgesetInfo & info = edgeset_info_map[id];
1429 :
1430 0 : elem->build_edge_ptr(edge, pair.second.first);
1431 :
1432 0 : ++info.num_edges;
1433 0 : info.edge_elem_types.insert(edge->type());
1434 :
1435 0 : if (verbosity > 1)
1436 0 : info.bbox.union_with(edge->loose_bounding_box());
1437 : }
1438 :
1439 : // Add edgeset info
1440 284 : if (!global || this->processor_id() == 0)
1441 : {
1442 48 : oss << "\n " << (global ? "" : "Local ") << "Mesh Edgesets:\n";
1443 48 : if (edgeset_ids.empty())
1444 48 : oss << " None\n";
1445 : }
1446 :
1447 8 : const auto & edgeset_name_map = this->get_boundary_info().get_edgeset_name_map();
1448 284 : for (const auto id : edgeset_ids)
1449 : {
1450 0 : EdgesetInfo & info = edgeset_info_map[id];
1451 :
1452 : // Reduce the local information for this edgeset if required
1453 0 : if (global)
1454 : {
1455 0 : this->comm().sum(info.num_edges);
1456 0 : this->comm().set_union(info.edge_elem_types, 0);
1457 0 : if (verbosity > 1)
1458 : {
1459 0 : this->comm().min(info.bbox.min());
1460 0 : this->comm().min(info.bbox.max());
1461 : }
1462 : }
1463 :
1464 0 : const bool has_name = edgeset_name_map.count(id) && edgeset_name_map.at(id).size();
1465 0 : const std::string name = has_name ? edgeset_name_map.at(id) : "";
1466 0 : if (global)
1467 0 : libmesh_assert(this->comm().verify(name));
1468 :
1469 0 : if (global ? this->processor_id() == 0 : info.num_edges > 0)
1470 : {
1471 0 : oss << " Edgeset " << id;
1472 0 : if (has_name)
1473 0 : oss << " (" << name << ")";
1474 0 : oss << ", " << info.num_edges << " " << (global ? "" : "local ") << "edges ("
1475 0 : << elem_type_helper(info.edge_elem_types) << ")\n";
1476 :
1477 0 : if (verbosity > 1)
1478 : {
1479 0 : oss << " " << (global ? "Bounding" : "Local bounding") << " box minimum: "
1480 0 : << info.bbox.min() << "\n"
1481 0 : << " " << (global ? "Bounding" : "Local bounding") << " box maximum: "
1482 0 : << info.bbox.max() << "\n"
1483 0 : << " " << (global ? "Bounding" : "Local bounding") << " box delta: "
1484 0 : << (info.bbox.max() - info.bbox.min()) << "\n";
1485 : }
1486 : }
1487 : }
1488 :
1489 : // Reduce the block IDs and block names
1490 16 : std::set<subdomain_id_type> subdomains;
1491 194032 : for (const Elem * elem : this->active_element_ptr_range())
1492 92640 : if (include_object(*elem))
1493 49420 : subdomains.insert(elem->subdomain_id());
1494 284 : if (global)
1495 284 : this->comm().set_union(subdomains);
1496 :
1497 : // Accumulate local information for each subdomain
1498 0 : struct SubdomainInfo
1499 : {
1500 : std::size_t num_elems = 0;
1501 : Real volume = 0;
1502 : std::set<int> elem_types;
1503 : std::set<dof_id_type> active_node_ids;
1504 : #ifdef LIBMESH_ENABLE_AMR
1505 : std::size_t num_active_elems = 0;
1506 : #endif
1507 : BoundingBox bbox;
1508 : };
1509 16 : std::map<subdomain_id_type, SubdomainInfo> subdomain_info_map;
1510 194032 : for (const Elem * elem : this->element_ptr_range())
1511 92640 : if (include_object(*elem))
1512 : {
1513 49152 : SubdomainInfo & info = subdomain_info_map[elem->subdomain_id()];
1514 :
1515 49152 : ++info.num_elems;
1516 49152 : info.elem_types.insert(elem->type());
1517 :
1518 : #ifdef LIBMESH_ENABLE_AMR
1519 4096 : if (elem->active())
1520 49152 : ++info.num_active_elems;
1521 : #endif
1522 :
1523 442368 : for (const Node & node : elem->node_ref_range())
1524 392704 : if (include_object(node) && node.active())
1525 330608 : info.active_node_ids.insert(node.id());
1526 :
1527 49152 : if (verbosity > 1 && elem->active())
1528 : {
1529 49152 : info.volume += elem->volume();
1530 49152 : info.bbox.union_with(elem->loose_bounding_box());
1531 : }
1532 268 : }
1533 :
1534 : // Add subdomain info
1535 284 : oss << "\n " << (global ? "" : "Local ") << "Mesh Subdomains:\n";
1536 8 : const auto & subdomain_name_map = this->get_subdomain_name_map();
1537 568 : for (const auto id : subdomains)
1538 : {
1539 284 : SubdomainInfo & info = subdomain_info_map[id];
1540 :
1541 284 : auto num_active_nodes = info.active_node_ids.size();
1542 :
1543 : // Reduce the information for this subdomain if needed
1544 284 : if (global)
1545 : {
1546 284 : this->comm().sum(info.num_elems);
1547 : #ifdef LIBMESH_ENABLE_AMR
1548 284 : this->comm().sum(info.num_active_elems);
1549 : #endif
1550 284 : this->comm().sum(num_active_nodes);
1551 284 : this->comm().set_union(info.elem_types, 0);
1552 284 : if (verbosity > 1)
1553 : {
1554 292 : this->comm().min(info.bbox.min());
1555 292 : this->comm().max(info.bbox.max());
1556 284 : this->comm().sum(info.volume);
1557 : }
1558 : }
1559 284 : if (verbosity > 1)
1560 284 : volume += info.volume;
1561 :
1562 8 : const bool has_name = subdomain_name_map.count(id);
1563 292 : const std::string name = has_name ? subdomain_name_map.at(id) : "";
1564 16 : if (global)
1565 8 : libmesh_assert(this->comm().verify(name));
1566 :
1567 284 : if (!global || this->processor_id() == 0)
1568 : {
1569 48 : oss << " Subdomain " << id;
1570 48 : if (has_name)
1571 0 : oss << " (" << name << ")";
1572 48 : oss << ": " << info.num_elems << " " << (global ? "" : "local ") << "elems "
1573 56 : << "(" << elem_type_helper(info.elem_types);
1574 : #ifdef LIBMESH_ENABLE_AMR
1575 48 : oss << ", " << info.num_active_elems << " active";
1576 : #endif
1577 52 : oss << "), " << num_active_nodes << " " << (global ? "" : "local ") << "active nodes\n";
1578 48 : if (verbosity > 1)
1579 : {
1580 52 : oss << " " << (global ? "Volume" : "Local volume") << ": " << info.volume << "\n";
1581 48 : oss << " " << (global ? "Bounding" : "Local bounding") << " box minimum: "
1582 48 : << info.bbox.min() << "\n"
1583 48 : << " " << (global ? "Bounding" : "Local bounding") << " box maximum: "
1584 48 : << info.bbox.max() << "\n"
1585 48 : << " " << (global ? "Bounding" : "Local bounding") << " box delta: "
1586 48 : << (info.bbox.max() - info.bbox.min()) << "\n";
1587 : }
1588 : }
1589 : }
1590 :
1591 560 : oss << " " << (global ? "Global" : "Local") << " mesh volume = " << volume << "\n";
1592 :
1593 268 : }
1594 :
1595 15282 : return oss.str();
1596 13968 : }
1597 :
1598 :
1599 14844 : void MeshBase::print_info(std::ostream & os, const unsigned int verbosity /* = 0 */, const bool global /* = true */) const
1600 : {
1601 15282 : os << this->get_info(verbosity, global)
1602 438 : << std::endl;
1603 14844 : }
1604 :
1605 :
1606 0 : std::ostream & operator << (std::ostream & os, const MeshBase & m)
1607 : {
1608 0 : m.print_info(os);
1609 0 : return os;
1610 : }
1611 :
1612 :
1613 389084 : void MeshBase::partition (const unsigned int n_parts)
1614 : {
1615 : // If we get here and we have unpartitioned elements, we need that
1616 : // fixed.
1617 389084 : if (this->n_unpartitioned_elem() > 0)
1618 : {
1619 8824 : libmesh_assert (partitioner().get());
1620 8824 : libmesh_assert (this->is_serial());
1621 294155 : partitioner()->partition (*this, n_parts);
1622 : }
1623 : // A nullptr partitioner or a skip_partitioning(true) call or a
1624 : // skip_noncritical_partitioning(true) call means don't repartition;
1625 : // skip_noncritical_partitioning() checks all these.
1626 2968 : else if (!skip_noncritical_partitioning())
1627 : {
1628 89052 : partitioner()->partition (*this, n_parts);
1629 : }
1630 : else
1631 : {
1632 : // Adaptive coarsening may have "orphaned" nodes on processors
1633 : // whose elements no longer share them. We need to check for
1634 : // and possibly fix that.
1635 5877 : MeshTools::correct_node_proc_ids(*this);
1636 :
1637 : // Make sure locally cached partition count is correct
1638 5877 : this->recalculate_n_partitions();
1639 :
1640 : // Make sure any other locally cached data is correct
1641 5877 : this->update_post_partitioning();
1642 : }
1643 389084 : }
1644 :
1645 15948 : void MeshBase::all_second_order (const bool full_ordered)
1646 : {
1647 15948 : this->all_second_order_range(this->element_ptr_range(), full_ordered);
1648 15948 : }
1649 :
1650 13916 : void MeshBase::all_complete_order ()
1651 : {
1652 13916 : this->all_complete_order_range(this->element_ptr_range());
1653 13916 : }
1654 :
1655 5877 : unsigned int MeshBase::recalculate_n_partitions()
1656 : {
1657 : // This requires an inspection on every processor
1658 146 : parallel_object_only();
1659 :
1660 5877 : unsigned int max_proc_id=0;
1661 :
1662 217698 : for (const auto & elem : this->active_local_element_ptr_range())
1663 215361 : max_proc_id = std::max(max_proc_id, static_cast<unsigned int>(elem->processor_id()));
1664 :
1665 : // The number of partitions is one more than the max processor ID.
1666 5877 : _n_parts = max_proc_id+1;
1667 :
1668 5877 : this->comm().max(_n_parts);
1669 :
1670 5877 : return _n_parts;
1671 : }
1672 :
1673 :
1674 :
1675 487532 : std::unique_ptr<PointLocatorBase> MeshBase::sub_point_locator () const
1676 : {
1677 : // If there's no master point locator, then we need one.
1678 487532 : if (_point_locator.get() == nullptr)
1679 : {
1680 : // PointLocator construction may not be safe within threads
1681 504 : libmesh_assert(!Threads::in_threads);
1682 :
1683 : // And it may require parallel communication
1684 504 : parallel_object_only();
1685 :
1686 : #ifdef LIBMESH_ENABLE_NANOFLANN_POINTLOCATOR
1687 : _point_locator = PointLocatorBase::build(NANOFLANN, *this);
1688 : #else
1689 24074 : _point_locator = PointLocatorBase::build(TREE_ELEMENTS, *this);
1690 : #endif
1691 :
1692 12541 : if (_point_locator_close_to_point_tol > 0.)
1693 0 : _point_locator->set_close_to_point_tol(_point_locator_close_to_point_tol);
1694 : }
1695 :
1696 : // Otherwise there was a master point locator, and we can grab a
1697 : // sub-locator easily.
1698 : return
1699 : #ifdef LIBMESH_ENABLE_NANOFLANN_POINTLOCATOR
1700 : PointLocatorBase::build(NANOFLANN, *this, _point_locator.get());
1701 : #else
1702 487532 : PointLocatorBase::build(TREE_ELEMENTS, *this, _point_locator.get());
1703 : #endif
1704 : }
1705 :
1706 :
1707 :
1708 1921238 : void MeshBase::clear_point_locator ()
1709 : {
1710 42274 : _point_locator.reset(nullptr);
1711 1921238 : }
1712 :
1713 :
1714 :
1715 0 : void MeshBase::set_count_lower_dim_elems_in_point_locator(bool count_lower_dim_elems)
1716 : {
1717 0 : _count_lower_dim_elems_in_point_locator = count_lower_dim_elems;
1718 0 : }
1719 :
1720 :
1721 :
1722 6385477 : bool MeshBase::get_count_lower_dim_elems_in_point_locator() const
1723 : {
1724 6385477 : return _count_lower_dim_elems_in_point_locator;
1725 : }
1726 :
1727 :
1728 :
1729 882727 : std::string & MeshBase::subdomain_name(subdomain_id_type id)
1730 : {
1731 882727 : return _block_id_to_name[id];
1732 : }
1733 :
1734 43761 : const std::string & MeshBase::subdomain_name(subdomain_id_type id) const
1735 : {
1736 : // An empty string to return when no matching subdomain name is found
1737 43761 : static const std::string empty;
1738 :
1739 43761 : if (const auto iter = _block_id_to_name.find(id);
1740 1602 : iter == _block_id_to_name.end())
1741 1396 : return empty;
1742 : else
1743 2476 : return iter->second;
1744 : }
1745 :
1746 :
1747 :
1748 :
1749 0 : subdomain_id_type MeshBase::get_id_by_name(std::string_view name) const
1750 : {
1751 : // Linear search over the map values.
1752 0 : for (const auto & [sbd_id, sbd_name] : _block_id_to_name)
1753 0 : if (sbd_name == name)
1754 0 : return sbd_id;
1755 :
1756 : // If we made it here without returning, we don't have a subdomain
1757 : // with the requested name, so return Elem::invalid_subdomain_id.
1758 0 : return Elem::invalid_subdomain_id;
1759 : }
1760 :
1761 : #ifdef LIBMESH_ENABLE_DEPRECATED
1762 0 : void MeshBase::cache_elem_dims()
1763 : {
1764 : libmesh_deprecated();
1765 :
1766 0 : this->cache_elem_data();
1767 0 : }
1768 : #endif // LIBMESH_ENABLE_DEPRECATED
1769 :
1770 464662 : void MeshBase::cache_elem_data()
1771 : {
1772 : // This requires an inspection on every processor
1773 14288 : parallel_object_only();
1774 :
1775 : // Need to clear containers first in case all elements of a
1776 : // particular dimension/order/subdomain have been deleted.
1777 28544 : _elem_dims.clear();
1778 28544 : _elem_default_orders.clear();
1779 28544 : _mesh_subdomains.clear();
1780 464662 : _supported_nodal_order = MAXIMUM;
1781 :
1782 41305250 : for (const auto & elem : this->active_element_ptr_range())
1783 : {
1784 40390214 : _elem_dims.insert(cast_int<unsigned char>(elem->dim()));
1785 40390214 : _elem_default_orders.insert(elem->default_order());
1786 40390214 : _mesh_subdomains.insert(elem->subdomain_id());
1787 40390214 : _supported_nodal_order =
1788 40390214 : static_cast<Order>
1789 80780428 : (std::min(static_cast<int>(_supported_nodal_order),
1790 40810066 : static_cast<int>(elem->supported_nodal_order())));
1791 436118 : }
1792 :
1793 464662 : if (!this->is_serial())
1794 : {
1795 : // Some different dimension/order/subdomain elements may only live
1796 : // on other processors
1797 122928 : this->comm().set_union(_elem_dims);
1798 122928 : this->comm().set_union(_elem_default_orders);
1799 122928 : this->comm().min(_supported_nodal_order);
1800 122928 : this->comm().set_union(_mesh_subdomains);
1801 : }
1802 :
1803 : // If the largest element dimension found is larger than the current
1804 : // _spatial_dimension, increase _spatial_dimension.
1805 464662 : unsigned int max_dim = this->mesh_dimension();
1806 464662 : if (max_dim > _spatial_dimension)
1807 29661 : _spatial_dimension = cast_int<unsigned char>(max_dim);
1808 :
1809 : // _spatial_dimension may need to increase from 1->2 or 2->3 if the
1810 : // mesh is full of 1D elements but they are not x-aligned, or the
1811 : // mesh is full of 2D elements but they are not in the x-y plane.
1812 : // If the mesh is x-aligned or x-y planar, we will end up checking
1813 : // every node's coordinates and not breaking out of the loop
1814 : // early...
1815 464662 : if (_spatial_dimension < LIBMESH_DIM)
1816 : {
1817 41540458 : for (const auto & node : this->node_ptr_range())
1818 : {
1819 : // Note: the exact floating point comparison is intentional,
1820 : // we don't want to get tripped up by tolerances.
1821 21609135 : if ((*node)(0) != 0. && _spatial_dimension < 1)
1822 0 : _spatial_dimension = 1;
1823 :
1824 21609135 : if ((*node)(1) != 0. && _spatial_dimension < 2)
1825 : {
1826 2195 : _spatial_dimension = 2;
1827 : #if LIBMESH_DIM == 2
1828 : // If libmesh is compiled in 2D mode, this is the
1829 : // largest spatial dimension possible so we can break
1830 : // out.
1831 : break;
1832 : #endif
1833 : }
1834 :
1835 : #if LIBMESH_DIM > 2
1836 21609135 : if ((*node)(2) != 0.)
1837 : {
1838 : // Spatial dimension can't get any higher than this, so
1839 : // we can break out.
1840 5387 : _spatial_dimension = 3;
1841 5387 : break;
1842 : }
1843 : #endif
1844 240624 : }
1845 : }
1846 464662 : }
1847 :
1848 :
1849 41988 : void MeshBase::sync_subdomain_name_map()
1850 : {
1851 : // This requires every processor
1852 1404 : parallel_object_only();
1853 :
1854 41988 : this->comm().set_union(_block_id_to_name);
1855 41988 : }
1856 :
1857 :
1858 412566 : void MeshBase::detect_interior_parents()
1859 : {
1860 : // This requires an inspection on every processor
1861 12586 : parallel_object_only();
1862 :
1863 : // Check if the mesh contains mixed dimensions. If so, then we may
1864 : // have interior parents to set. Otherwise return.
1865 412566 : if (this->elem_dimensions().size() == 1)
1866 403994 : return;
1867 :
1868 : // Do we have interior parent pointers going to a different mesh?
1869 : // If so then we'll still check to make sure that's the only place
1870 : // they go, so we can libmesh_not_implemented() if not.
1871 460 : const bool separate_interior_mesh = (&(this->interior_mesh()) != this);
1872 :
1873 : // This map will be used to set interior parents
1874 460 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> node_to_elem;
1875 :
1876 20428190 : for (const auto & elem : this->element_ptr_range())
1877 : {
1878 : // Populating the node_to_elem map, same as MeshTools::build_nodes_to_elem_map
1879 51952965 : for (auto n : make_range(elem->n_vertices()))
1880 : {
1881 1206440 : libmesh_assert_less (elem->id(), this->max_elem_id());
1882 :
1883 42647269 : node_to_elem[elem->node_id(n)].push_back(elem->id());
1884 : }
1885 8112 : }
1886 :
1887 : // Automatically set interior parents
1888 20428190 : for (const auto & element : this->element_ptr_range())
1889 : {
1890 : // Ignore an 3D element or an element that already has an interior parent
1891 10512136 : if (element->dim()>=LIBMESH_DIM || element->interior_parent())
1892 10069663 : continue;
1893 :
1894 : // Start by generating a SET of elements that are dim+1 to the current
1895 : // element at each vertex of the current element, thus ignoring interior nodes.
1896 : // If one of the SET of elements is empty, then we will not have an interior parent
1897 : // since an interior parent must be connected to all vertices of the current element
1898 475757 : std::vector<std::set<dof_id_type>> neighbors( element->n_vertices() );
1899 :
1900 16642 : bool found_interior_parents = false;
1901 :
1902 445668 : for (auto n : make_range(element->n_vertices()))
1903 : {
1904 461378 : std::vector<dof_id_type> & element_ids = node_to_elem[element->node_id(n)];
1905 1918983 : for (const auto & eid : element_ids)
1906 1474309 : if (this->elem_ref(eid).dim() == element->dim()+1)
1907 8733 : neighbors[n].insert(eid);
1908 :
1909 461378 : if (neighbors[n].size()>0)
1910 : {
1911 90 : found_interior_parents = true;
1912 : }
1913 : else
1914 : {
1915 : // We have found an empty set, no reason to continue
1916 : // Ensure we set this flag to false before the break since it could have
1917 : // been set to true for previous vertex
1918 16614 : found_interior_parents = false;
1919 16614 : break;
1920 : }
1921 : }
1922 :
1923 : // If we have successfully generated a set of elements for each vertex, we will compare
1924 : // the set for vertex 0 will the sets for the vertices until we find a id that exists in
1925 : // all sets. If found, this is our an interior parent id. The interior parent id found
1926 : // will be the lowest element id if there is potential for multiple interior parents.
1927 442473 : if (found_interior_parents)
1928 : {
1929 56 : std::set<dof_id_type> & neighbors_0 = neighbors[0];
1930 1633 : for (const auto & interior_parent_id : neighbors_0)
1931 : {
1932 40 : found_interior_parents = false;
1933 2627 : for (auto n : make_range(1u, element->n_vertices()))
1934 : {
1935 1633 : if (neighbors[n].count(interior_parent_id))
1936 : {
1937 34 : found_interior_parents = true;
1938 : }
1939 : else
1940 : {
1941 12 : found_interior_parents = false;
1942 12 : break;
1943 : }
1944 : }
1945 :
1946 1420 : if (found_interior_parents)
1947 : {
1948 781 : element->set_interior_parent(this->elem_ptr(interior_parent_id));
1949 22 : break;
1950 : }
1951 : }
1952 :
1953 : // Do we have a mixed dimensional mesh that contains some of
1954 : // its own interior parents, but we already expect to have
1955 : // interior parents on a different mesh? That's going to
1956 : // take some work to support if anyone needs it.
1957 994 : if (separate_interior_mesh)
1958 0 : libmesh_not_implemented_msg
1959 : ("interior_parent() values in multiple meshes are unsupported.");
1960 : }
1961 417301 : }
1962 : }
1963 :
1964 :
1965 :
1966 0 : void MeshBase::set_point_locator_close_to_point_tol(Real val)
1967 : {
1968 0 : _point_locator_close_to_point_tol = val;
1969 0 : if (_point_locator)
1970 : {
1971 0 : if (val > 0.)
1972 0 : _point_locator->set_close_to_point_tol(val);
1973 : else
1974 0 : _point_locator->unset_close_to_point_tol();
1975 : }
1976 0 : }
1977 :
1978 :
1979 :
1980 426 : Real MeshBase::get_point_locator_close_to_point_tol() const
1981 : {
1982 426 : return _point_locator_close_to_point_tol;
1983 : }
1984 :
1985 :
1986 :
1987 4696 : void MeshBase::size_elem_extra_integers()
1988 : {
1989 272 : const std::size_t new_size = _elem_integer_names.size();
1990 35102 : for (auto elem : this->element_ptr_range())
1991 17729 : elem->add_extra_integers(new_size, _elem_integer_default_values);
1992 4696 : }
1993 :
1994 :
1995 :
1996 14472 : void MeshBase::size_node_extra_integers()
1997 : {
1998 846 : const std::size_t new_size = _node_integer_names.size();
1999 216621 : for (auto node : this->node_ptr_range())
2000 110576 : node->add_extra_integers(new_size, _node_integer_default_values);
2001 14472 : }
2002 :
2003 :
2004 : std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
2005 23398 : MeshBase::merge_extra_integer_names(const MeshBase & other)
2006 : {
2007 776 : std::pair<std::vector<unsigned int>, std::vector<unsigned int>> returnval;
2008 23398 : returnval.first = this->add_elem_integers(other._elem_integer_names, true, &other._elem_integer_default_values);
2009 23398 : returnval.second = this->add_node_integers(other._node_integer_names, true, &other._node_integer_default_values);
2010 23398 : return returnval;
2011 : }
2012 :
2013 :
2014 :
2015 : void
2016 426 : MeshBase::post_dofobject_moves(MeshBase && other_mesh)
2017 : {
2018 : // Now that all the DofObject moving is done, we can move the GhostingFunctor objects
2019 : // which include the _default_ghosting,_ghosting_functors and _shared_functors. We also need
2020 : // to set the mesh object associated with these functors to the assignee mesh.
2021 :
2022 : // _default_ghosting
2023 12 : _default_ghosting = std::move(other_mesh._default_ghosting);
2024 426 : _default_ghosting->set_mesh(this);
2025 :
2026 : // _ghosting_functors
2027 426 : _ghosting_functors = std::move(other_mesh._ghosting_functors);
2028 :
2029 852 : for (const auto gf : _ghosting_functors )
2030 : {
2031 426 : gf->set_mesh(this);
2032 : }
2033 :
2034 : // _shared_functors
2035 12 : _shared_functors = std::move(other_mesh._shared_functors);
2036 :
2037 426 : for (const auto & sf : _shared_functors )
2038 : {
2039 0 : (sf.second)->set_mesh(this);
2040 : }
2041 :
2042 : // _constraint_rows
2043 12 : _constraint_rows = std::move(other_mesh._constraint_rows);
2044 :
2045 426 : if (other_mesh.partitioner())
2046 426 : _partitioner = std::move(other_mesh.partitioner());
2047 426 : }
2048 :
2049 :
2050 : void
2051 21765 : MeshBase::copy_cached_data(const MeshBase & other_mesh)
2052 : {
2053 21765 : this->_spatial_dimension = other_mesh._spatial_dimension;
2054 730 : this->_elem_dims = other_mesh._elem_dims;
2055 730 : this->_elem_default_orders = other_mesh._elem_default_orders;
2056 21765 : this->_supported_nodal_order = other_mesh._supported_nodal_order;
2057 730 : this->_mesh_subdomains = other_mesh._mesh_subdomains;
2058 21765 : }
2059 :
2060 :
2061 6250 : bool MeshBase::nodes_and_elements_equal(const MeshBase & other_mesh) const
2062 : {
2063 1188084 : for (const auto & other_node : other_mesh.node_ptr_range())
2064 : {
2065 689575 : const Node * node = this->query_node_ptr(other_node->id());
2066 689575 : if (!node)
2067 0 : return false;
2068 689575 : if (*other_node != *node)
2069 0 : return false;
2070 5306 : }
2071 1188084 : for (const auto & node : this->node_ptr_range())
2072 689575 : if (!other_mesh.query_node_ptr(node->id()))
2073 5306 : return false;
2074 :
2075 719574 : for (const auto & other_elem : other_mesh.element_ptr_range())
2076 : {
2077 523191 : const Elem * elem = this->query_elem_ptr(other_elem->id());
2078 523191 : if (!elem)
2079 0 : return false;
2080 523191 : if (!other_elem->topologically_equal(*elem))
2081 6 : return false;
2082 5306 : }
2083 719298 : for (const auto & elem : this->element_ptr_range())
2084 523050 : if (!other_mesh.query_elem_ptr(elem->id()))
2085 5177 : return false;
2086 :
2087 6109 : return true;
2088 : }
2089 :
2090 :
2091 736 : dof_id_type MeshBase::n_constraint_rows() const
2092 : {
2093 736 : dof_id_type n_local_rows=0, n_unpartitioned_rows=0;
2094 736 : for (const auto & [node, node_constraints] : _constraint_rows)
2095 : {
2096 : // Unpartitioned nodes
2097 0 : if (node->processor_id() == DofObject::invalid_processor_id)
2098 0 : n_unpartitioned_rows++;
2099 0 : else if (node->processor_id() == this->processor_id())
2100 0 : n_local_rows++;
2101 : }
2102 :
2103 736 : this->comm().sum(n_local_rows);
2104 :
2105 736 : return n_unpartitioned_rows + n_local_rows;
2106 : }
2107 :
2108 :
2109 : void
2110 21765 : MeshBase::copy_constraint_rows(const MeshBase & other_mesh)
2111 : {
2112 1460 : LOG_SCOPE("copy_constraint_rows(mesh)", "MeshBase");
2113 :
2114 1444 : _constraint_rows.clear();
2115 :
2116 730 : const auto & other_constraint_rows = other_mesh.get_constraint_rows();
2117 70193 : for (const auto & [other_node, other_node_constraints] : other_constraint_rows)
2118 : {
2119 48428 : const Node * const our_node = this->node_ptr(other_node->id());
2120 6416 : constraint_rows_mapped_type our_node_constraints;
2121 219637 : for (const auto & [other_inner_key_pair, constraint_value] : other_node_constraints)
2122 : {
2123 9048 : const auto & [other_elem, local_node_id] = other_inner_key_pair;
2124 171209 : const Elem * const our_elem = this->elem_ptr(other_elem->id());
2125 171209 : our_node_constraints.emplace_back(std::make_pair(our_elem, local_node_id), constraint_value);
2126 : }
2127 48428 : _constraint_rows[our_node] = std::move(our_node_constraints);
2128 : }
2129 21765 : }
2130 :
2131 :
2132 : template <typename T>
2133 : void
2134 565 : MeshBase::copy_constraint_rows(const SparseMatrix<T> & constraint_operator,
2135 : bool precondition_constraint_operator)
2136 : {
2137 32 : LOG_SCOPE("copy_constraint_rows(mat)", "MeshBase");
2138 :
2139 16 : this->_constraint_rows.clear();
2140 :
2141 : // We're not going to support doing this distributed yet; it'd be
2142 : // pointless unless we temporarily had a linear partitioning to
2143 : // better match the constraint operator.
2144 597 : MeshSerializer serialize(*this);
2145 :
2146 : // Our current mesh should already reflect the desired assembly space
2147 565 : libmesh_error_msg_if(this->n_nodes() != constraint_operator.m(),
2148 : "Constraint operator matrix with " <<
2149 : constraint_operator.m() <<
2150 : "rows does not match this mesh with " <<
2151 : this->n_nodes() << " nodes");
2152 :
2153 : // First, find what new unconstrained DoFs we need to add. We can't
2154 : // iterate over columns in a SparseMatrix, so we'll iterate over
2155 : // rows and keep track of columns.
2156 :
2157 : // If we have nodes that will work unconstrained, keep track of
2158 : // their node ids and corresponding column indices.
2159 : // existing_unconstrained_nodes[column_id] = node_id
2160 32 : std::map<dof_id_type, dof_id_type> existing_unconstrained_columns;
2161 32 : std::set<dof_id_type> existing_unconstrained_nodes;
2162 :
2163 : // In case we need new nodes, keep track of their columns.
2164 : // columns[j][k] will be the kth row index and value of column j
2165 : typedef
2166 : std::unordered_map<dof_id_type,
2167 : std::vector<std::pair<dof_id_type, Real>>>
2168 : columns_type;
2169 1114 : columns_type columns(constraint_operator.n());
2170 :
2171 : // If we need to precondition the constraint operator (e.g. it's an
2172 : // unpreconditioned extraction operator for a Flex IGA matrix),
2173 : // we'll want to keep track of the sum of each column, because we'll
2174 : // be dividing each column by that sum (Jacobi preconditioning on
2175 : // the right, which then leads to symmetric preconditioning on a
2176 : // physics Jacobian).
2177 32 : std::unordered_map<dof_id_type, Real> column_sums;
2178 :
2179 : // Work in parallel, though we'll have to sync shortly
2180 4028 : for (auto i : make_range(constraint_operator.row_start(),
2181 533 : constraint_operator.row_stop()))
2182 : {
2183 618 : std::vector<numeric_index_type> indices;
2184 618 : std::vector<T> values;
2185 :
2186 3463 : constraint_operator.get_row(i, indices, values);
2187 309 : libmesh_assert_equal_to(indices.size(), values.size());
2188 :
2189 3847 : if (indices.size() == 1 &&
2190 825 : values[0] == T(1))
2191 : {
2192 : // If we have multiple simple Ui=Uj constraints, let the
2193 : // first one be our "unconstrained" node and let the others
2194 : // be constrained to it.
2195 810 : if (existing_unconstrained_columns.find(indices[0]) !=
2196 138 : existing_unconstrained_columns.end())
2197 : {
2198 36 : const auto j = indices[0];
2199 36 : columns[j].emplace_back(i, 1);
2200 : }
2201 : else
2202 : {
2203 708 : existing_unconstrained_nodes.insert(i);
2204 774 : existing_unconstrained_columns.emplace(indices[0],i);
2205 : }
2206 : }
2207 : else
2208 22028 : for (auto jj : index_range(indices))
2209 : {
2210 19375 : const auto j = indices[jj];
2211 21134 : const Real coef = libmesh_real(values[jj]);
2212 1759 : libmesh_assert_equal_to(coef, values[jj]);
2213 19375 : columns[j].emplace_back(i, coef);
2214 : }
2215 : }
2216 :
2217 : // Merge data from different processors' slabs of the matrix
2218 565 : this->comm().set_union(existing_unconstrained_nodes);
2219 565 : this->comm().set_union(existing_unconstrained_columns);
2220 :
2221 48 : std::vector<columns_type> all_columns;
2222 565 : this->comm().allgather(columns, all_columns);
2223 :
2224 16 : columns.clear();
2225 6106 : for (auto p : index_range(all_columns))
2226 65601 : for (auto & [j, subcol] : all_columns[p])
2227 183429 : for (auto [i, v] : subcol)
2228 123369 : columns[j].emplace_back(i,v);
2229 :
2230 : // Keep track of elements on which unconstrained nodes exist, and
2231 : // their local node indices.
2232 : // node_to_elem_ptrs[node] = [elem_id, local_node_num]
2233 32 : std::unordered_map<const Node *, std::pair<dof_id_type, unsigned int>> node_to_elem_ptrs;
2234 :
2235 : // Find elements attached to any existing nodes that will stay
2236 : // unconstrained. We'll also build a subdomain set here so we don't
2237 : // have to assert that the mesh is already prepared before we pick a
2238 : // new subdomain for any NodeElems we need to add.
2239 32 : std::set<subdomain_id_type> subdomain_ids;
2240 145181 : for (const Elem * elem : this->element_ptr_range())
2241 : {
2242 49738 : subdomain_ids.insert(elem->subdomain_id());
2243 244572 : for (auto n : make_range(elem->n_nodes()))
2244 : {
2245 194834 : const Node * node = elem->node_ptr(n);
2246 194834 : if (existing_unconstrained_nodes.count(node->id()))
2247 12175 : node_to_elem_ptrs.emplace(node, std::make_pair(elem->id(), n));
2248 : }
2249 : }
2250 :
2251 565 : const subdomain_id_type new_sbd_id = *subdomain_ids.rbegin() + 1;
2252 :
2253 23189 : for (auto j : make_range(constraint_operator.n()))
2254 : {
2255 : // If we already have a good node for this then we're done
2256 4800 : if (existing_unconstrained_columns.count(j))
2257 4668 : continue;
2258 :
2259 : // Get a half-decent spot to place a new NodeElem for
2260 : // unconstrained DoF(s) here. Getting a *fully*-decent spot
2261 : // would require finding a Moore-Penrose pseudoinverse, and I'm
2262 : // not going to do that, but scaling a transpose will at least
2263 : // get us a little uniqueness to make visualization reasonable.
2264 264 : Point newpt;
2265 264 : Real total_scaling = 0;
2266 264 : unsigned int total_entries = 0;
2267 :
2268 : // We'll get a decent initial pid choice here too, if only to
2269 : // aid in later repartitioning.
2270 528 : std::map<processor_id_type, int> pids;
2271 :
2272 264 : auto & column = columns[j];
2273 122230 : for (auto [i, r] : column)
2274 : {
2275 112988 : Node & constrained_node = this->node_ref(i);
2276 112988 : const Point constrained_pt = constrained_node;
2277 3228 : newpt += r*constrained_pt;
2278 112988 : total_scaling += r;
2279 112988 : ++total_entries;
2280 112988 : ++pids[constrained_node.processor_id()];
2281 : }
2282 :
2283 9242 : if (precondition_constraint_operator)
2284 0 : column_sums[j] = total_scaling;
2285 :
2286 9242 : libmesh_error_msg_if
2287 : (!total_entries,
2288 : "Empty column " << j <<
2289 : " found in constraint operator matrix");
2290 :
2291 : // If we have *cancellation* here then we can end up dividing by
2292 : // zero; try just evenly scaling across all constrained node
2293 : // points instead.
2294 9242 : if (total_scaling > TOLERANCE)
2295 264 : newpt /= total_scaling;
2296 : else
2297 0 : newpt /= total_entries;
2298 :
2299 9242 : Node *n = this->add_point(newpt);
2300 9506 : std::unique_ptr<Elem> elem = Elem::build(NODEELEM);
2301 9242 : elem->set_node(0, n);
2302 9242 : elem->subdomain_id() = new_sbd_id;
2303 :
2304 9506 : Elem * added_elem = this->add_elem(std::move(elem));
2305 9242 : this->_elem_dims.insert(0);
2306 9242 : this->_elem_default_orders.insert(added_elem->default_order());
2307 9242 : this->_supported_nodal_order =
2308 8714 : static_cast<Order>
2309 18484 : (std::min(static_cast<int>(this->_supported_nodal_order),
2310 9242 : static_cast<int>(added_elem->supported_nodal_order())));
2311 8978 : this->_mesh_subdomains.insert(new_sbd_id);
2312 9242 : node_to_elem_ptrs.emplace(n, std::make_pair(added_elem->id(), 0));
2313 9242 : existing_unconstrained_columns.emplace(j,n->id());
2314 :
2315 : // Repartition the new objects *after* adding them, so a
2316 : // DistributedMesh doesn't get confused and think you're not
2317 : // adding them on all processors at once.
2318 264 : int n_pids = 0;
2319 43912 : for (auto [pid, count] : pids)
2320 34670 : if (count >= n_pids)
2321 : {
2322 324 : n_pids = count;
2323 15469 : added_elem->processor_id() = pid;
2324 15469 : n->processor_id() = pid;
2325 : }
2326 : }
2327 :
2328 : // Calculate constraint rows in an indexed form that's easy for us
2329 : // to allgather
2330 : std::unordered_map<dof_id_type,
2331 : std::vector<std::pair<std::pair<dof_id_type, unsigned int>,Real>>>
2332 32 : indexed_constraint_rows;
2333 :
2334 4028 : for (auto i : make_range(constraint_operator.row_start(),
2335 533 : constraint_operator.row_stop()))
2336 : {
2337 951 : if (existing_unconstrained_nodes.count(i))
2338 774 : continue;
2339 :
2340 486 : std::vector<numeric_index_type> indices;
2341 486 : std::vector<T> values;
2342 :
2343 2689 : constraint_operator.get_row(i, indices, values);
2344 :
2345 486 : std::vector<std::pair<std::pair<dof_id_type, unsigned int>, Real>> constraint_row;
2346 :
2347 22100 : for (auto jj : index_range(indices))
2348 : {
2349 21173 : const dof_id_type node_id =
2350 : existing_unconstrained_columns[indices[jj]];
2351 :
2352 19411 : Node & constraining_node = this->node_ref(node_id);
2353 :
2354 1762 : libmesh_assert(node_to_elem_ptrs.count(&constraining_node));
2355 :
2356 19411 : auto p = node_to_elem_ptrs[&constraining_node];
2357 :
2358 19382 : Real coef = libmesh_real(values[jj]);
2359 1762 : libmesh_assert_equal_to(coef, values[jj]);
2360 :
2361 : // If we're preconditioning and we created a nodeelem then
2362 : // we can scale the meaning of that nodeelem's value to give
2363 : // us a better-conditioned matrix after the constraints are
2364 : // applied.
2365 19411 : if (precondition_constraint_operator)
2366 0 : if (auto sum_it = column_sums.find(indices[jj]);
2367 0 : sum_it != column_sums.end())
2368 : {
2369 0 : const Real scaling = sum_it->second;
2370 :
2371 0 : if (scaling > TOLERANCE)
2372 0 : coef /= scaling;
2373 : }
2374 :
2375 21173 : constraint_row.emplace_back(std::make_pair(p, coef));
2376 : }
2377 :
2378 243 : indexed_constraint_rows.emplace(i, std::move(constraint_row));
2379 : }
2380 :
2381 565 : this->comm().set_union(indexed_constraint_rows);
2382 :
2383 : // Add constraint rows as mesh constraint rows
2384 17591 : for (auto & [node_id, indexed_row] : indexed_constraint_rows)
2385 : {
2386 17026 : Node * constrained_node = this->node_ptr(node_id);
2387 :
2388 972 : constraint_rows_mapped_type constraint_row;
2389 :
2390 140395 : for (auto [p, coef] : indexed_row)
2391 : {
2392 123369 : const Elem * elem = this->elem_ptr(p.first);
2393 7048 : constraint_row.emplace_back
2394 126893 : (std::make_pair(std::make_pair(elem, p.second), coef));
2395 : }
2396 :
2397 16540 : this->_constraint_rows.emplace(constrained_node,
2398 486 : std::move(constraint_row));
2399 : }
2400 1098 : }
2401 :
2402 :
2403 0 : void MeshBase::print_constraint_rows(std::ostream & os,
2404 : bool print_nonlocal) const
2405 : {
2406 0 : parallel_object_only();
2407 :
2408 : std::string local_constraints =
2409 0 : this->get_local_constraints(print_nonlocal);
2410 :
2411 0 : if (this->processor_id())
2412 : {
2413 0 : this->comm().send(0, local_constraints);
2414 : }
2415 : else
2416 : {
2417 0 : os << "Processor 0:\n";
2418 0 : os << local_constraints;
2419 :
2420 0 : for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
2421 : {
2422 0 : this->comm().receive(p, local_constraints);
2423 0 : os << "Processor " << p << ":\n";
2424 0 : os << local_constraints;
2425 : }
2426 : }
2427 0 : }
2428 :
2429 :
2430 :
2431 0 : std::string MeshBase::get_local_constraints(bool print_nonlocal) const
2432 : {
2433 0 : std::ostringstream os;
2434 :
2435 0 : if (print_nonlocal)
2436 0 : os << "All ";
2437 : else
2438 0 : os << "Local ";
2439 :
2440 0 : os << "Mesh Constraint Rows:"
2441 0 : << std::endl;
2442 :
2443 0 : for (const auto & [node, row] : _constraint_rows)
2444 : {
2445 0 : const bool local = (node->processor_id() == this->processor_id());
2446 :
2447 : // Skip non-local dofs if requested
2448 0 : if (!print_nonlocal && !local)
2449 0 : continue;
2450 :
2451 0 : os << "Constraints for " << (local ? "Local" : "Ghost") << " Node " << node->id()
2452 0 : << ": \t";
2453 :
2454 0 : for (const auto & [elem_and_node, coef] : row)
2455 0 : os << " ((" << elem_and_node.first->id() << ',' << elem_and_node.second << "), " << coef << ")\t";
2456 :
2457 0 : os << std::endl;
2458 : }
2459 :
2460 0 : return os.str();
2461 0 : }
2462 :
2463 :
2464 :
2465 :
2466 : // Explicit instantiations for our template function
2467 : template LIBMESH_EXPORT void
2468 : MeshBase::copy_constraint_rows(const SparseMatrix<Real> & constraint_operator,
2469 : bool precondition_constraint_operator);
2470 :
2471 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
2472 : template LIBMESH_EXPORT void
2473 : MeshBase::copy_constraint_rows(const SparseMatrix<Complex> & constraint_operator,
2474 : bool precondition_constraint_operator);
2475 : #endif
2476 :
2477 :
2478 : } // namespace libMesh
|
