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