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