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);
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  {
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 set sorted by pointer, 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  _ghosting_functors.erase(&ghosting_functor);
951 
952  if (const auto it = _shared_functors.find(&ghosting_functor);
953  it != _shared_functors.end())
954  _shared_functors.erase(it);
955 }
956 
957 
958 
959 void MeshBase::subdomain_ids (std::set<subdomain_id_type> & ids, const bool global /* = true */) const
960 {
961  // This requires an inspection on every processor
962  if (global)
963  parallel_object_only();
964 
965  ids.clear();
966 
967  for (const auto & elem : this->active_local_element_ptr_range())
968  ids.insert(elem->subdomain_id());
969 
970  if (global)
971  {
972  // Only include the unpartitioned elements if the user requests the global IDs.
973  // In the case of the local subdomain IDs, it doesn't make sense to include the
974  // unpartitioned elements because said elements do not have a sense of locality.
975  for (const auto & elem : this->active_unpartitioned_element_ptr_range())
976  ids.insert(elem->subdomain_id());
977 
978  // Some subdomains may only live on other processors
979  this->comm().set_union(ids);
980  }
981 }
982 
983 
984 
986 {
987  // We now have all elements and nodes redistributed; our ghosting
988  // functors should be ready to redistribute and/or recompute any
989  // cached data they use too.
990  for (auto & gf : as_range(this->ghosting_functors_begin(),
991  this->ghosting_functors_end()))
992  gf->redistribute();
993 }
994 
995 
996 
998 {
999  // This requires an inspection on every processor
1000  parallel_object_only();
1001 
1002  std::set<subdomain_id_type> ids;
1003 
1004  this->subdomain_ids (ids);
1005 
1006  return cast_int<subdomain_id_type>(ids.size());
1007 }
1008 
1009 
1010 
1012 {
1013  std::set<subdomain_id_type> ids;
1014 
1015  this->subdomain_ids (ids, /* global = */ false);
1016 
1017  return cast_int<subdomain_id_type>(ids.size());
1018 }
1019 
1020 
1021 
1022 
1024 {
1025  // We're either counting a processor's nodes or unpartitioned
1026  // nodes
1027  libmesh_assert (proc_id < this->n_processors() ||
1028  proc_id == DofObject::invalid_processor_id);
1029 
1030  return static_cast<dof_id_type>(std::distance (this->pid_nodes_begin(proc_id),
1031  this->pid_nodes_end (proc_id)));
1032 }
1033 
1034 
1035 
1037 {
1038  // We're either counting a processor's elements or unpartitioned
1039  // elements
1040  libmesh_assert (proc_id < this->n_processors() ||
1041  proc_id == DofObject::invalid_processor_id);
1042 
1043  return static_cast<dof_id_type>(std::distance (this->pid_elements_begin(proc_id),
1044  this->pid_elements_end (proc_id)));
1045 }
1046 
1047 
1048 
1050 {
1051  libmesh_assert_less (proc_id, this->n_processors());
1052  return static_cast<dof_id_type>(std::distance (this->active_pid_elements_begin(proc_id),
1053  this->active_pid_elements_end (proc_id)));
1054 }
1055 
1056 
1057 
1059 {
1060  dof_id_type ne=0;
1061 
1062  for (const auto & elem : this->element_ptr_range())
1063  ne += elem->n_sub_elem();
1064 
1065  return ne;
1066 }
1067 
1068 
1069 
1071 {
1072  dof_id_type ne=0;
1073 
1074  for (const auto & elem : this->active_element_ptr_range())
1075  ne += elem->n_sub_elem();
1076 
1077  return ne;
1078 }
1079 
1080 
1081 
1082 std::string MeshBase::get_info(const unsigned int verbosity /* = 0 */, const bool global /* = true */) const
1083 {
1084  std::ostringstream oss;
1085 
1086  oss << " Mesh Information:" << '\n';
1087 
1088  if (!_elem_dims.empty())
1089  {
1090  oss << " elem_dimensions()={";
1091  std::copy(_elem_dims.begin(),
1092  --_elem_dims.end(), // --end() is valid if the set is non-empty
1093  std::ostream_iterator<unsigned int>(oss, ", "));
1094  oss << cast_int<unsigned int>(*_elem_dims.rbegin());
1095  oss << "}\n";
1096  }
1097 
1098  if (!_elem_default_orders.empty())
1099  {
1100  oss << " elem_default_orders()={";
1101  std::transform(_elem_default_orders.begin(),
1102  --_elem_default_orders.end(),
1103  std::ostream_iterator<std::string>(oss, ", "),
1104  [](Order o)
1105  { return Utility::enum_to_string<Order>(o); });
1106  oss << Utility::enum_to_string<Order>(*_elem_default_orders.rbegin());
1107  oss << "}\n";
1108  }
1109 
1110  oss << " supported_nodal_order()=" << this->supported_nodal_order() << '\n'
1111  << " spatial_dimension()=" << this->spatial_dimension() << '\n'
1112  << " n_nodes()=" << this->n_nodes() << '\n'
1113  << " n_local_nodes()=" << this->n_local_nodes() << '\n'
1114  << " n_elem()=" << this->n_elem() << '\n'
1115  << " n_local_elem()=" << this->n_local_elem() << '\n';
1116 #ifdef LIBMESH_ENABLE_AMR
1117  oss << " n_active_elem()=" << this->n_active_elem() << '\n';
1118 #endif
1119  if (global)
1120  oss << " n_subdomains()=" << static_cast<std::size_t>(this->n_subdomains()) << '\n';
1121  else
1122  oss << " n_local_subdomains()= " << static_cast<std::size_t>(this->n_local_subdomains()) << '\n';
1123  oss << " n_elemsets()=" << static_cast<std::size_t>(this->n_elemsets()) << '\n';
1124  if (!_elemset_codes.empty())
1125  oss << " n_elemset_codes=" << _elemset_codes.size() << '\n';
1126  oss << " n_partitions()=" << static_cast<std::size_t>(this->n_partitions()) << '\n'
1127  << " n_processors()=" << static_cast<std::size_t>(this->n_processors()) << '\n'
1128  << " n_threads()=" << static_cast<std::size_t>(libMesh::n_threads()) << '\n'
1129  << " processor_id()=" << static_cast<std::size_t>(this->processor_id()) << '\n'
1130  << " is_prepared()=" << (this->is_prepared() ? "true" : "false") << '\n'
1131  << " is_replicated()=" << (this->is_replicated() ? "true" : "false") << '\n';
1132 
1133  if (verbosity > 0)
1134  {
1135  if (global)
1136  {
1137  libmesh_parallel_only(this->comm());
1138  if (this->processor_id() != 0)
1139  oss << "\n Detailed global get_info() (verbosity > 0) is reduced and output to only rank 0.";
1140  }
1141 
1142  // Helper for printing element types
1143  const auto elem_type_helper = [](const std::set<int> & elem_types) {
1144  std::stringstream ss;
1145  for (auto it = elem_types.begin(); it != elem_types.end();)
1146  {
1147  ss << Utility::enum_to_string((ElemType)*it);
1148  if (++it != elem_types.end())
1149  ss << ", ";
1150  }
1151  return ss.str();
1152  };
1153 
1154  // Helper for whether or not the given DofObject is to be included. If we're doing
1155  // a global reduction, we also count unpartitioned objects on rank 0.
1156  const auto include_object = [this, &global](const DofObject & dof_object) {
1157  return this->processor_id() == dof_object.processor_id() ||
1158  (global &&
1159  this->processor_id() == 0 &&
1160  dof_object.processor_id() == DofObject::invalid_processor_id);
1161  };
1162 
1163  Real volume = 0;
1164 
1165  // Add bounding box information
1166  const auto bbox = global ? MeshTools::create_bounding_box(*this) : MeshTools::create_local_bounding_box(*this);
1167  if (!global || this->processor_id() == 0)
1168  oss << "\n " << (global ? "" : "Local ") << "Mesh Bounding Box:\n"
1169  << " Minimum: " << bbox.min() << "\n"
1170  << " Maximum: " << bbox.max() << "\n"
1171  << " Delta: " << (bbox.max() - bbox.min()) << "\n";
1172 
1173  // Obtain the global or local element types
1174  std::set<int> elem_types;
1175  for (const Elem * elem : this->active_local_element_ptr_range())
1176  elem_types.insert(elem->type());
1177  if (global)
1178  {
1179  // Pick up unpartitioned elems on rank 0
1180  if (this->processor_id() == 0)
1181  for (const Elem * elem : this->active_unpartitioned_element_ptr_range())
1182  elem_types.insert(elem->type());
1183 
1184  this->comm().set_union(elem_types);
1185  }
1186 
1187  // Add element types
1188  if (!global || this->processor_id() == 0)
1189  oss << "\n " << (global ? "" : "Local ") << "Mesh Element Type(s):\n "
1190  << elem_type_helper(elem_types) << "\n";
1191 
1192  // Reduce the nodeset ids
1193  auto nodeset_ids = this->get_boundary_info().get_node_boundary_ids();
1194  if (global)
1195  this->comm().set_union(nodeset_ids);
1196 
1197  // Accumulate local information for each nodeset
1198  struct NodesetInfo
1199  {
1200  std::size_t num_nodes = 0;
1201  BoundingBox bbox;
1202  };
1203  std::map<boundary_id_type, NodesetInfo> nodeset_info_map;
1204  for (const auto & [node, id] : this->get_boundary_info().get_nodeset_map())
1205  {
1206  if (!include_object(*node))
1207  continue;
1208 
1209  NodesetInfo & info = nodeset_info_map[id];
1210 
1211  ++info.num_nodes;
1212 
1213  if (verbosity > 1)
1214  info.bbox.union_with(*node);
1215  }
1216 
1217  // Add nodeset info
1218  if (!global || this->processor_id() == 0)
1219  {
1220  oss << "\n " << (global ? "" : "Local ") << "Mesh Nodesets:\n";
1221  if (nodeset_ids.empty())
1222  oss << " None\n";
1223  }
1224 
1225  const auto & nodeset_name_map = this->get_boundary_info().get_nodeset_name_map();
1226  for (const auto id : nodeset_ids)
1227  {
1228  NodesetInfo & info = nodeset_info_map[id];
1229 
1230  // Reduce the local information for this nodeset if required
1231  if (global)
1232  {
1233  this->comm().sum(info.num_nodes);
1234  if (verbosity > 1)
1235  {
1236  this->comm().min(info.bbox.min());
1237  this->comm().max(info.bbox.max());
1238  }
1239  }
1240 
1241  const bool has_name = nodeset_name_map.count(id) && nodeset_name_map.at(id).size();
1242  const std::string name = has_name ? nodeset_name_map.at(id) : "";
1243  if (global)
1244  libmesh_assert(this->comm().verify(name));
1245 
1246  if (global ? this->processor_id() == 0 : info.num_nodes > 0)
1247  {
1248  oss << " Nodeset " << id;
1249  if (has_name)
1250  oss << " (" << name << ")";
1251  oss << ", " << info.num_nodes << " " << (global ? "" : "local ") << "nodes\n";
1252 
1253  if (verbosity > 1)
1254  {
1255  oss << " " << (global ? "Bounding" : "Local bounding") << " box minimum: "
1256  << info.bbox.min() << "\n"
1257  << " " << (global ? "Bounding" : "Local bounding") << " box maximum: "
1258  << info.bbox.max() << "\n"
1259  << " " << (global ? "Bounding" : "Local bounding") << " box delta: "
1260  << (info.bbox.max() - info.bbox.min()) << "\n";
1261  }
1262  }
1263  }
1264 
1265  // Reduce the sideset ids
1266  auto sideset_ids = this->get_boundary_info().get_side_boundary_ids();
1267  if (global)
1268  this->comm().set_union(sideset_ids);
1269 
1270  // Accumulate local information for each sideset
1271  struct SidesetInfo
1272  {
1273  std::size_t num_sides = 0;
1274  Real volume = 0;
1275  std::set<int> side_elem_types;
1276  std::set<int> elem_types;
1277  std::set<dof_id_type> elem_ids;
1278  std::set<dof_id_type> node_ids;
1279  BoundingBox bbox;
1280  };
1281  ElemSideBuilder side_builder;
1282  std::map<boundary_id_type, SidesetInfo> sideset_info_map;
1283  for (const auto & pair : this->get_boundary_info().get_sideset_map())
1284  {
1285  const Elem * elem = pair.first;
1286  if (!include_object(*elem))
1287  continue;
1288 
1289  const auto id = pair.second.second;
1290  SidesetInfo & info = sideset_info_map[id];
1291 
1292  const auto s = pair.second.first;
1293  const Elem & side = side_builder(*elem, s);
1294 
1295  ++info.num_sides;
1296  info.side_elem_types.insert(side.type());
1297  info.elem_types.insert(elem->type());
1298  info.elem_ids.insert(elem->id());
1299 
1300  for (const Node & node : side.node_ref_range())
1301  if (include_object(node))
1302  info.node_ids.insert(node.id());
1303 
1304  if (verbosity > 1)
1305  {
1306  info.volume += side.volume();
1307  info.bbox.union_with(side.loose_bounding_box());
1308  }
1309  }
1310 
1311  // Add sideset info
1312  if (!global || this->processor_id() == 0)
1313  {
1314  oss << "\n " << (global ? "" : "Local ") << "Mesh Sidesets:\n";
1315  if (sideset_ids.empty())
1316  oss << " None\n";
1317  }
1318  const auto & sideset_name_map = this->get_boundary_info().get_sideset_name_map();
1319  for (const auto id : sideset_ids)
1320  {
1321  SidesetInfo & info = sideset_info_map[id];
1322 
1323  auto num_elems = info.elem_ids.size();
1324  auto num_nodes = info.node_ids.size();
1325 
1326  // Reduce the local information for this sideset if required
1327  if (global)
1328  {
1329  this->comm().sum(info.num_sides);
1330  this->comm().set_union(info.side_elem_types, 0);
1331  this->comm().sum(num_elems);
1332  this->comm().set_union(info.elem_types, 0);
1333  this->comm().sum(num_nodes);
1334  if (verbosity > 1)
1335  {
1336  this->comm().sum(info.volume);
1337  this->comm().min(info.bbox.min());
1338  this->comm().max(info.bbox.max());
1339  }
1340  }
1341 
1342  const bool has_name = sideset_name_map.count(id) && sideset_name_map.at(id).size();
1343  const std::string name = has_name ? sideset_name_map.at(id) : "";
1344  if (global)
1345  libmesh_assert(this->comm().verify(name));
1346 
1347  if (global ? this->processor_id() == 0 : info.num_sides > 0)
1348  {
1349  oss << " Sideset " << id;
1350  if (has_name)
1351  oss << " (" << name << ")";
1352  oss << ", " << info.num_sides << " sides (" << elem_type_helper(info.side_elem_types) << ")"
1353  << ", " << num_elems << " " << (global ? "" : "local ") << "elems (" << elem_type_helper(info.elem_types) << ")"
1354  << ", " << num_nodes << " " << (global ? "" : "local ") << "nodes\n";
1355 
1356  if (verbosity > 1)
1357  {
1358  oss << " " << (global ? "Side" : "Local side") << " volume: " << info.volume << "\n"
1359  << " " << (global ? "Bounding" : "Local bounding") << " box minimum: "
1360  << info.bbox.min() << "\n"
1361  << " " << (global ? "Bounding" : "Local bounding") << " box maximum: "
1362  << info.bbox.max() << "\n"
1363  << " " << (global ? "Bounding" : "Local bounding") << " box delta: "
1364  << (info.bbox.max() - info.bbox.min()) << "\n";
1365  }
1366  }
1367  }
1368 
1369  // Reduce the edgeset ids
1370  auto edgeset_ids = this->get_boundary_info().get_edge_boundary_ids();
1371  if (global)
1372  this->comm().set_union(edgeset_ids);
1373 
1374  // Accumulate local information for each edgeset
1375  struct EdgesetInfo
1376  {
1377  std::size_t num_edges = 0;
1378  std::set<int> edge_elem_types;
1379  BoundingBox bbox;
1380  };
1381  std::map<boundary_id_type, EdgesetInfo> edgeset_info_map;
1382  std::unique_ptr<const Elem> edge;
1383 
1384  for (const auto & pair : this->get_boundary_info().get_edgeset_map())
1385  {
1386  const Elem * elem = pair.first;
1387  if (!include_object(*elem))
1388  continue;
1389 
1390  const auto id = pair.second.second;
1391  EdgesetInfo & info = edgeset_info_map[id];
1392 
1393  elem->build_edge_ptr(edge, pair.second.first);
1394 
1395  ++info.num_edges;
1396  info.edge_elem_types.insert(edge->type());
1397 
1398  if (verbosity > 1)
1399  info.bbox.union_with(edge->loose_bounding_box());
1400  }
1401 
1402  // Add edgeset info
1403  if (!global || this->processor_id() == 0)
1404  {
1405  oss << "\n " << (global ? "" : "Local ") << "Mesh Edgesets:\n";
1406  if (edgeset_ids.empty())
1407  oss << " None\n";
1408  }
1409 
1410  const auto & edgeset_name_map = this->get_boundary_info().get_edgeset_name_map();
1411  for (const auto id : edgeset_ids)
1412  {
1413  EdgesetInfo & info = edgeset_info_map[id];
1414 
1415  // Reduce the local information for this edgeset if required
1416  if (global)
1417  {
1418  this->comm().sum(info.num_edges);
1419  this->comm().set_union(info.edge_elem_types, 0);
1420  if (verbosity > 1)
1421  {
1422  this->comm().min(info.bbox.min());
1423  this->comm().min(info.bbox.max());
1424  }
1425  }
1426 
1427  const bool has_name = edgeset_name_map.count(id) && edgeset_name_map.at(id).size();
1428  const std::string name = has_name ? edgeset_name_map.at(id) : "";
1429  if (global)
1430  libmesh_assert(this->comm().verify(name));
1431 
1432  if (global ? this->processor_id() == 0 : info.num_edges > 0)
1433  {
1434  oss << " Edgeset " << id;
1435  if (has_name)
1436  oss << " (" << name << ")";
1437  oss << ", " << info.num_edges << " " << (global ? "" : "local ") << "edges ("
1438  << elem_type_helper(info.edge_elem_types) << ")\n";
1439 
1440  if (verbosity > 1)
1441  {
1442  oss << " " << (global ? "Bounding" : "Local bounding") << " box minimum: "
1443  << info.bbox.min() << "\n"
1444  << " " << (global ? "Bounding" : "Local bounding") << " box maximum: "
1445  << info.bbox.max() << "\n"
1446  << " " << (global ? "Bounding" : "Local bounding") << " box delta: "
1447  << (info.bbox.max() - info.bbox.min()) << "\n";
1448  }
1449  }
1450  }
1451 
1452  // Reduce the block IDs and block names
1453  std::set<subdomain_id_type> subdomains;
1454  for (const Elem * elem : this->active_element_ptr_range())
1455  if (include_object(*elem))
1456  subdomains.insert(elem->subdomain_id());
1457  if (global)
1458  this->comm().set_union(subdomains);
1459 
1460  // Accumulate local information for each subdomain
1461  struct SubdomainInfo
1462  {
1463  std::size_t num_elems = 0;
1464  Real volume = 0;
1465  std::set<int> elem_types;
1466  std::set<dof_id_type> active_node_ids;
1467 #ifdef LIBMESH_ENABLE_AMR
1468  std::size_t num_active_elems = 0;
1469 #endif
1470  BoundingBox bbox;
1471  };
1472  std::map<subdomain_id_type, SubdomainInfo> subdomain_info_map;
1473  for (const Elem * elem : this->element_ptr_range())
1474  if (include_object(*elem))
1475  {
1476  SubdomainInfo & info = subdomain_info_map[elem->subdomain_id()];
1477 
1478  ++info.num_elems;
1479  info.elem_types.insert(elem->type());
1480 
1481 #ifdef LIBMESH_ENABLE_AMR
1482  if (elem->active())
1483  ++info.num_active_elems;
1484 #endif
1485 
1486  for (const Node & node : elem->node_ref_range())
1487  if (include_object(node) && node.active())
1488  info.active_node_ids.insert(node.id());
1489 
1490  if (verbosity > 1 && elem->active())
1491  {
1492  info.volume += elem->volume();
1493  info.bbox.union_with(elem->loose_bounding_box());
1494  }
1495  }
1496 
1497  // Add subdomain info
1498  oss << "\n " << (global ? "" : "Local ") << "Mesh Subdomains:\n";
1499  const auto & subdomain_name_map = this->get_subdomain_name_map();
1500  for (const auto id : subdomains)
1501  {
1502  SubdomainInfo & info = subdomain_info_map[id];
1503 
1504  auto num_active_nodes = info.active_node_ids.size();
1505 
1506  // Reduce the information for this subdomain if needed
1507  if (global)
1508  {
1509  this->comm().sum(info.num_elems);
1510 #ifdef LIBMESH_ENABLE_AMR
1511  this->comm().sum(info.num_active_elems);
1512 #endif
1513  this->comm().sum(num_active_nodes);
1514  this->comm().set_union(info.elem_types, 0);
1515  if (verbosity > 1)
1516  {
1517  this->comm().min(info.bbox.min());
1518  this->comm().max(info.bbox.max());
1519  this->comm().sum(info.volume);
1520  }
1521  }
1522  if (verbosity > 1)
1523  volume += info.volume;
1524 
1525  const bool has_name = subdomain_name_map.count(id);
1526  const std::string name = has_name ? subdomain_name_map.at(id) : "";
1527  if (global)
1528  libmesh_assert(this->comm().verify(name));
1529 
1530  if (!global || this->processor_id() == 0)
1531  {
1532  oss << " Subdomain " << id;
1533  if (has_name)
1534  oss << " (" << name << ")";
1535  oss << ": " << info.num_elems << " " << (global ? "" : "local ") << "elems "
1536  << "(" << elem_type_helper(info.elem_types);
1537 #ifdef LIBMESH_ENABLE_AMR
1538  oss << ", " << info.num_active_elems << " active";
1539 #endif
1540  oss << "), " << num_active_nodes << " " << (global ? "" : "local ") << "active nodes\n";
1541  if (verbosity > 1)
1542  {
1543  oss << " " << (global ? "Volume" : "Local volume") << ": " << info.volume << "\n";
1544  oss << " " << (global ? "Bounding" : "Local bounding") << " box minimum: "
1545  << info.bbox.min() << "\n"
1546  << " " << (global ? "Bounding" : "Local bounding") << " box maximum: "
1547  << info.bbox.max() << "\n"
1548  << " " << (global ? "Bounding" : "Local bounding") << " box delta: "
1549  << (info.bbox.max() - info.bbox.min()) << "\n";
1550  }
1551  }
1552  }
1553 
1554  oss << " " << (global ? "Global" : "Local") << " mesh volume = " << volume << "\n";
1555 
1556  }
1557 
1558  return oss.str();
1559 }
1560 
1561 
1562 void MeshBase::print_info(std::ostream & os, const unsigned int verbosity /* = 0 */, const bool global /* = true */) const
1563 {
1564  os << this->get_info(verbosity, global)
1565  << std::endl;
1566 }
1567 
1568 
1569 std::ostream & operator << (std::ostream & os, const MeshBase & m)
1570 {
1571  m.print_info(os);
1572  return os;
1573 }
1574 
1575 
1576 void MeshBase::partition (const unsigned int n_parts)
1577 {
1578  // If we get here and we have unpartitioned elements, we need that
1579  // fixed.
1580  if (this->n_unpartitioned_elem() > 0)
1581  {
1582  libmesh_assert (partitioner().get());
1583  libmesh_assert (this->is_serial());
1584  partitioner()->partition (*this, n_parts);
1585  }
1586  // A nullptr partitioner or a skip_partitioning(true) call or a
1587  // skip_noncritical_partitioning(true) call means don't repartition;
1588  // skip_noncritical_partitioning() checks all these.
1589  else if (!skip_noncritical_partitioning())
1590  {
1591  partitioner()->partition (*this, n_parts);
1592  }
1593  else
1594  {
1595  // Adaptive coarsening may have "orphaned" nodes on processors
1596  // whose elements no longer share them. We need to check for
1597  // and possibly fix that.
1599 
1600  // Make sure locally cached partition count is correct
1601  this->recalculate_n_partitions();
1602 
1603  // Make sure any other locally cached data is correct
1604  this->update_post_partitioning();
1605  }
1606 }
1607 
1608 void MeshBase::all_second_order (const bool full_ordered)
1609 {
1610  this->all_second_order_range(this->element_ptr_range(), full_ordered);
1611 }
1612 
1614 {
1615  this->all_complete_order_range(this->element_ptr_range());
1616 }
1617 
1619 {
1620  // This requires an inspection on every processor
1621  parallel_object_only();
1622 
1623  unsigned int max_proc_id=0;
1624 
1625  for (const auto & elem : this->active_local_element_ptr_range())
1626  max_proc_id = std::max(max_proc_id, static_cast<unsigned int>(elem->processor_id()));
1627 
1628  // The number of partitions is one more than the max processor ID.
1629  _n_parts = max_proc_id+1;
1630 
1631  this->comm().max(_n_parts);
1632 
1633  return _n_parts;
1634 }
1635 
1636 
1637 
1638 std::unique_ptr<PointLocatorBase> MeshBase::sub_point_locator () const
1639 {
1640  // If there's no master point locator, then we need one.
1641  if (_point_locator.get() == nullptr)
1642  {
1643  // PointLocator construction may not be safe within threads
1645 
1646  // And it may require parallel communication
1647  parallel_object_only();
1648 
1649 #ifdef LIBMESH_ENABLE_NANOFLANN_POINTLOCATOR
1651 #else
1653 #endif
1654 
1656  _point_locator->set_close_to_point_tol(_point_locator_close_to_point_tol);
1657  }
1658 
1659  // Otherwise there was a master point locator, and we can grab a
1660  // sub-locator easily.
1661  return
1662 #ifdef LIBMESH_ENABLE_NANOFLANN_POINTLOCATOR
1664 #else
1666 #endif
1667 }
1668 
1669 
1670 
1672 {
1673  _point_locator.reset(nullptr);
1674 }
1675 
1676 
1677 
1679 {
1680  _count_lower_dim_elems_in_point_locator = count_lower_dim_elems;
1681 }
1682 
1683 
1684 
1686 {
1688 }
1689 
1690 
1691 
1693 {
1694  return _block_id_to_name[id];
1695 }
1696 
1697 const std::string & MeshBase::subdomain_name(subdomain_id_type id) const
1698 {
1699  // An empty string to return when no matching subdomain name is found
1700  static const std::string empty;
1701 
1702  if (const auto iter = _block_id_to_name.find(id);
1703  iter == _block_id_to_name.end())
1704  return empty;
1705  else
1706  return iter->second;
1707 }
1708 
1709 
1710 
1711 
1713 {
1714  // Linear search over the map values.
1715  for (const auto & [sbd_id, sbd_name] : _block_id_to_name)
1716  if (sbd_name == name)
1717  return sbd_id;
1718 
1719  // If we made it here without returning, we don't have a subdomain
1720  // with the requested name, so return Elem::invalid_subdomain_id.
1722 }
1723 
1724 #ifdef LIBMESH_ENABLE_DEPRECATED
1726 {
1727  libmesh_deprecated();
1728 
1729  this->cache_elem_data();
1730 }
1731 #endif // LIBMESH_ENABLE_DEPRECATED
1732 
1734 {
1735  // This requires an inspection on every processor
1736  parallel_object_only();
1737 
1738  // Need to clear containers first in case all elements of a
1739  // particular dimension/order/subdomain have been deleted.
1740  _elem_dims.clear();
1741  _elem_default_orders.clear();
1742  _mesh_subdomains.clear();
1744 
1745  for (const auto & elem : this->active_element_ptr_range())
1746  {
1747  _elem_dims.insert(cast_int<unsigned char>(elem->dim()));
1748  _elem_default_orders.insert(elem->default_order());
1749  _mesh_subdomains.insert(elem->subdomain_id());
1751  static_cast<Order>
1752  (std::min(static_cast<int>(_supported_nodal_order),
1753  static_cast<int>(elem->supported_nodal_order())));
1754  }
1755 
1756  if (!this->is_serial())
1757  {
1758  // Some different dimension/order/subdomain elements may only live
1759  // on other processors
1760  this->comm().set_union(_elem_dims);
1762  this->comm().min(_supported_nodal_order);
1763  this->comm().set_union(_mesh_subdomains);
1764  }
1765 
1766  // If the largest element dimension found is larger than the current
1767  // _spatial_dimension, increase _spatial_dimension.
1768  unsigned int max_dim = this->mesh_dimension();
1769  if (max_dim > _spatial_dimension)
1770  _spatial_dimension = cast_int<unsigned char>(max_dim);
1771 
1772  // _spatial_dimension may need to increase from 1->2 or 2->3 if the
1773  // mesh is full of 1D elements but they are not x-aligned, or the
1774  // mesh is full of 2D elements but they are not in the x-y plane.
1775  // If the mesh is x-aligned or x-y planar, we will end up checking
1776  // every node's coordinates and not breaking out of the loop
1777  // early...
1778  if (_spatial_dimension < LIBMESH_DIM)
1779  {
1780  for (const auto & node : this->node_ptr_range())
1781  {
1782  // Note: the exact floating point comparison is intentional,
1783  // we don't want to get tripped up by tolerances.
1784  if ((*node)(0) != 0. && _spatial_dimension < 1)
1785  _spatial_dimension = 1;
1786 
1787  if ((*node)(1) != 0. && _spatial_dimension < 2)
1788  {
1789  _spatial_dimension = 2;
1790 #if LIBMESH_DIM == 2
1791  // If libmesh is compiled in 2D mode, this is the
1792  // largest spatial dimension possible so we can break
1793  // out.
1794  break;
1795 #endif
1796  }
1797 
1798 #if LIBMESH_DIM > 2
1799  if ((*node)(2) != 0.)
1800  {
1801  // Spatial dimension can't get any higher than this, so
1802  // we can break out.
1803  _spatial_dimension = 3;
1804  break;
1805  }
1806 #endif
1807  }
1808  }
1809 }
1810 
1812 {
1813  // This requires an inspection on every processor
1814  parallel_object_only();
1815 
1816  // Check if the mesh contains mixed dimensions. If so, then we may
1817  // have interior parents to set. Otherwise return.
1818  if (this->elem_dimensions().size() == 1)
1819  return;
1820 
1821  // Do we have interior parent pointers going to a different mesh?
1822  // If so then we'll still check to make sure that's the only place
1823  // they go, so we can libmesh_not_implemented() if not.
1824  const bool separate_interior_mesh = (&(this->interior_mesh()) != this);
1825 
1826  // This map will be used to set interior parents
1827  std::unordered_map<dof_id_type, std::vector<dof_id_type>> node_to_elem;
1828 
1829  for (const auto & elem : this->element_ptr_range())
1830  {
1831  // Populating the node_to_elem map, same as MeshTools::build_nodes_to_elem_map
1832  for (auto n : make_range(elem->n_vertices()))
1833  {
1834  libmesh_assert_less (elem->id(), this->max_elem_id());
1835 
1836  node_to_elem[elem->node_id(n)].push_back(elem->id());
1837  }
1838  }
1839 
1840  // Automatically set interior parents
1841  for (const auto & element : this->element_ptr_range())
1842  {
1843  // Ignore an 3D element or an element that already has an interior parent
1844  if (element->dim()>=LIBMESH_DIM || element->interior_parent())
1845  continue;
1846 
1847  // Start by generating a SET of elements that are dim+1 to the current
1848  // element at each vertex of the current element, thus ignoring interior nodes.
1849  // If one of the SET of elements is empty, then we will not have an interior parent
1850  // since an interior parent must be connected to all vertices of the current element
1851  std::vector<std::set<dof_id_type>> neighbors( element->n_vertices() );
1852 
1853  bool found_interior_parents = false;
1854 
1855  for (auto n : make_range(element->n_vertices()))
1856  {
1857  std::vector<dof_id_type> & element_ids = node_to_elem[element->node_id(n)];
1858  for (const auto & eid : element_ids)
1859  if (this->elem_ref(eid).dim() == element->dim()+1)
1860  neighbors[n].insert(eid);
1861 
1862  if (neighbors[n].size()>0)
1863  {
1864  found_interior_parents = true;
1865  }
1866  else
1867  {
1868  // We have found an empty set, no reason to continue
1869  // Ensure we set this flag to false before the break since it could have
1870  // been set to true for previous vertex
1871  found_interior_parents = false;
1872  break;
1873  }
1874  }
1875 
1876  // If we have successfully generated a set of elements for each vertex, we will compare
1877  // the set for vertex 0 will the sets for the vertices until we find a id that exists in
1878  // all sets. If found, this is our an interior parent id. The interior parent id found
1879  // will be the lowest element id if there is potential for multiple interior parents.
1880  if (found_interior_parents)
1881  {
1882  std::set<dof_id_type> & neighbors_0 = neighbors[0];
1883  for (const auto & interior_parent_id : neighbors_0)
1884  {
1885  found_interior_parents = false;
1886  for (auto n : make_range(1u, element->n_vertices()))
1887  {
1888  if (neighbors[n].count(interior_parent_id))
1889  {
1890  found_interior_parents = true;
1891  }
1892  else
1893  {
1894  found_interior_parents = false;
1895  break;
1896  }
1897  }
1898 
1899  if (found_interior_parents)
1900  {
1901  element->set_interior_parent(this->elem_ptr(interior_parent_id));
1902  break;
1903  }
1904  }
1905 
1906  // Do we have a mixed dimensional mesh that contains some of
1907  // its own interior parents, but we already expect to have
1908  // interior parents on a different mesh? That's going to
1909  // take some work to support if anyone needs it.
1910  if (separate_interior_mesh)
1911  libmesh_not_implemented_msg
1912  ("interior_parent() values in multiple meshes are unsupported.");
1913  }
1914  }
1915 }
1916 
1917 
1918 
1920 {
1922  if (_point_locator)
1923  {
1924  if (val > 0.)
1925  _point_locator->set_close_to_point_tol(val);
1926  else
1927  _point_locator->unset_close_to_point_tol();
1928  }
1929 }
1930 
1931 
1932 
1934 {
1936 }
1937 
1938 
1939 
1941 {
1942  const std::size_t new_size = _elem_integer_names.size();
1943  for (auto elem : this->element_ptr_range())
1944  elem->add_extra_integers(new_size, _elem_integer_default_values);
1945 }
1946 
1947 
1948 
1950 {
1951  const std::size_t new_size = _node_integer_names.size();
1952  for (auto node : this->node_ptr_range())
1953  node->add_extra_integers(new_size, _node_integer_default_values);
1954 }
1955 
1956 
1957 std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
1959 {
1960  std::pair<std::vector<unsigned int>, std::vector<unsigned int>> returnval;
1961  returnval.first = this->add_elem_integers(other._elem_integer_names, true, &other._elem_integer_default_values);
1962  returnval.second = this->add_node_integers(other._node_integer_names, true, &other._node_integer_default_values);
1963  return returnval;
1964 }
1965 
1966 
1967 
1968 void
1970 {
1971  // Now that all the DofObject moving is done, we can move the GhostingFunctor objects
1972  // which include the _default_ghosting,_ghosting_functors and _shared_functors. We also need
1973  // to set the mesh object associated with these functors to the assignee mesh.
1974 
1975  // _default_ghosting
1976  _default_ghosting = std::move(other_mesh._default_ghosting);
1977  _default_ghosting->set_mesh(this);
1978 
1979  // _ghosting_functors
1980  _ghosting_functors = std::move(other_mesh._ghosting_functors);
1981 
1982  for (const auto gf : _ghosting_functors )
1983  {
1984  gf->set_mesh(this);
1985  }
1986 
1987  // _shared_functors
1988  _shared_functors = std::move(other_mesh._shared_functors);
1989 
1990  for (const auto & sf : _shared_functors )
1991  {
1992  (sf.second)->set_mesh(this);
1993  }
1994 
1995  // _constraint_rows
1996  _constraint_rows = std::move(other_mesh._constraint_rows);
1997 
1998  if (other_mesh.partitioner())
1999  _partitioner = std::move(other_mesh.partitioner());
2000 }
2001 
2002 
2003 void
2005 {
2006  this->_spatial_dimension = other_mesh._spatial_dimension;
2007  this->_elem_dims = other_mesh._elem_dims;
2008  this->_elem_default_orders = other_mesh._elem_default_orders;
2010  this->_mesh_subdomains = other_mesh._mesh_subdomains;
2011 }
2012 
2013 
2014 bool MeshBase::nodes_and_elements_equal(const MeshBase & other_mesh) const
2015 {
2016  for (const auto & other_node : other_mesh.node_ptr_range())
2017  {
2018  const Node * node = this->query_node_ptr(other_node->id());
2019  if (!node)
2020  return false;
2021  if (*other_node != *node)
2022  return false;
2023  }
2024  for (const auto & node : this->node_ptr_range())
2025  if (!other_mesh.query_node_ptr(node->id()))
2026  return false;
2027 
2028  for (const auto & other_elem : other_mesh.element_ptr_range())
2029  {
2030  const Elem * elem = this->query_elem_ptr(other_elem->id());
2031  if (!elem)
2032  return false;
2033  if (!other_elem->topologically_equal(*elem))
2034  return false;
2035  }
2036  for (const auto & elem : this->element_ptr_range())
2037  if (!other_mesh.query_elem_ptr(elem->id()))
2038  return false;
2039 
2040  return true;
2041 }
2042 
2043 
2045 {
2046  dof_id_type n_local_rows=0, n_unpartitioned_rows=0;
2047  for (const auto & [node, node_constraints] : _constraint_rows)
2048  {
2049  // Unpartitioned nodes
2050  if (node->processor_id() == DofObject::invalid_processor_id)
2051  n_unpartitioned_rows++;
2052  else if (node->processor_id() == this->processor_id())
2053  n_local_rows++;
2054  }
2055 
2056  this->comm().sum(n_local_rows);
2057 
2058  return n_unpartitioned_rows + n_local_rows;
2059 }
2060 
2061 
2062 void
2064 {
2065  LOG_SCOPE("copy_constraint_rows(mesh)", "MeshBase");
2066 
2067  _constraint_rows.clear();
2068 
2069  const auto & other_constraint_rows = other_mesh.get_constraint_rows();
2070  for (const auto & [other_node, other_node_constraints] : other_constraint_rows)
2071  {
2072  const Node * const our_node = this->node_ptr(other_node->id());
2073  constraint_rows_mapped_type our_node_constraints;
2074  for (const auto & [other_inner_key_pair, constraint_value] : other_node_constraints)
2075  {
2076  const auto & [other_elem, local_node_id] = other_inner_key_pair;
2077  const Elem * const our_elem = this->elem_ptr(other_elem->id());
2078  our_node_constraints.emplace_back(std::make_pair(our_elem, local_node_id), constraint_value);
2079  }
2080  _constraint_rows[our_node] = std::move(our_node_constraints);
2081  }
2082 }
2083 
2084 
2085 template <typename T>
2086 void
2088  bool precondition_constraint_operator)
2089 {
2090  LOG_SCOPE("copy_constraint_rows(mat)", "MeshBase");
2091 
2092  this->_constraint_rows.clear();
2093 
2094  // We're not going to support doing this distributed yet; it'd be
2095  // pointless unless we temporarily had a linear partitioning to
2096  // better match the constraint operator.
2097  MeshSerializer serialize(*this);
2098 
2099  // Our current mesh should already reflect the desired assembly space
2100  libmesh_error_msg_if(this->n_nodes() != constraint_operator.m(),
2101  "Constraint operator matrix with " <<
2102  constraint_operator.m() <<
2103  "rows does not match this mesh with " <<
2104  this->n_nodes() << " nodes");
2105 
2106  // First, find what new unconstrained DoFs we need to add. We can't
2107  // iterate over columns in a SparseMatrix, so we'll iterate over
2108  // rows and keep track of columns.
2109 
2110  // If we have nodes that will work unconstrained, keep track of
2111  // their node ids and corresponding column indices.
2112  // existing_unconstrained_nodes[column_id] = node_id
2113  std::map<dof_id_type, dof_id_type> existing_unconstrained_columns;
2114  std::set<dof_id_type> existing_unconstrained_nodes;
2115 
2116  // In case we need new nodes, keep track of their columns.
2117  // columns[j][k] will be the kth row index and value of column j
2118  typedef
2119  std::unordered_map<dof_id_type,
2120  std::vector<std::pair<dof_id_type, Real>>>
2121  columns_type;
2122  columns_type columns(constraint_operator.n());
2123 
2124  // If we need to precondition the constraint operator (e.g. it's an
2125  // unpreconditioned extraction operator for a Flex IGA matrix),
2126  // we'll want to keep track of the sum of each column, because we'll
2127  // be dividing each column by that sum (Jacobi preconditioning on
2128  // the right, which then leads to symmetric preconditioning on a
2129  // physics Jacobian).
2130  std::unordered_map<dof_id_type, Real> column_sums;
2131 
2132  // Work in parallel, though we'll have to sync shortly
2133  for (auto i : make_range(constraint_operator.row_start(),
2134  constraint_operator.row_stop()))
2135  {
2136  std::vector<numeric_index_type> indices;
2137  std::vector<T> values;
2138 
2139  constraint_operator.get_row(i, indices, values);
2140  libmesh_assert_equal_to(indices.size(), values.size());
2141 
2142  if (indices.size() == 1 &&
2143  values[0] == T(1))
2144  {
2145  // If we have multiple simple Ui=Uj constraints, let the
2146  // first one be our "unconstrained" node and let the others
2147  // be constrained to it.
2148  if (existing_unconstrained_columns.find(indices[0]) !=
2149  existing_unconstrained_columns.end())
2150  {
2151  const auto j = indices[0];
2152  columns[j].emplace_back(i, 1);
2153  }
2154  else
2155  {
2156  existing_unconstrained_nodes.insert(i);
2157  existing_unconstrained_columns.emplace(indices[0],i);
2158  }
2159  }
2160  else
2161  for (auto jj : index_range(indices))
2162  {
2163  const auto j = indices[jj];
2164  const Real coef = libmesh_real(values[jj]);
2165  libmesh_assert_equal_to(coef, values[jj]);
2166  columns[j].emplace_back(i, coef);
2167  }
2168  }
2169 
2170  // Merge data from different processors' slabs of the matrix
2171  this->comm().set_union(existing_unconstrained_nodes);
2172  this->comm().set_union(existing_unconstrained_columns);
2173 
2174  std::vector<columns_type> all_columns;
2175  this->comm().allgather(columns, all_columns);
2176 
2177  columns.clear();
2178  for (auto p : index_range(all_columns))
2179  for (auto & [j, subcol] : all_columns[p])
2180  for (auto [i, v] : subcol)
2181  columns[j].emplace_back(i,v);
2182 
2183  // Keep track of elements on which unconstrained nodes exist, and
2184  // their local node indices.
2185  // node_to_elem_ptrs[node] = [elem_id, local_node_num]
2186  std::unordered_map<const Node *, std::pair<dof_id_type, unsigned int>> node_to_elem_ptrs;
2187 
2188  // Find elements attached to any existing nodes that will stay
2189  // unconstrained. We'll also build a subdomain set here so we don't
2190  // have to assert that the mesh is already prepared before we pick a
2191  // new subdomain for any NodeElems we need to add.
2192  std::set<subdomain_id_type> subdomain_ids;
2193  for (const Elem * elem : this->element_ptr_range())
2194  {
2195  subdomain_ids.insert(elem->subdomain_id());
2196  for (auto n : make_range(elem->n_nodes()))
2197  {
2198  const Node * node = elem->node_ptr(n);
2199  if (existing_unconstrained_nodes.count(node->id()))
2200  node_to_elem_ptrs.emplace(node, std::make_pair(elem->id(), n));
2201  }
2202  }
2203 
2204  const subdomain_id_type new_sbd_id = *subdomain_ids.rbegin() + 1;
2205 
2206  for (auto j : make_range(constraint_operator.n()))
2207  {
2208  // If we already have a good node for this then we're done
2209  if (existing_unconstrained_columns.count(j))
2210  continue;
2211 
2212  // Get a half-decent spot to place a new NodeElem for
2213  // unconstrained DoF(s) here. Getting a *fully*-decent spot
2214  // would require finding a Moore-Penrose pseudoinverse, and I'm
2215  // not going to do that, but scaling a transpose will at least
2216  // get us a little uniqueness to make visualization reasonable.
2217  Point newpt;
2218  Real total_scaling = 0;
2219  unsigned int total_entries = 0;
2220 
2221  // We'll get a decent initial pid choice here too, if only to
2222  // aid in later repartitioning.
2223  std::map<processor_id_type, int> pids;
2224 
2225  auto & column = columns[j];
2226  for (auto [i, r] : column)
2227  {
2228  Node & constrained_node = this->node_ref(i);
2229  const Point constrained_pt = constrained_node;
2230  newpt += r*constrained_pt;
2231  total_scaling += r;
2232  ++total_entries;
2233  ++pids[constrained_node.processor_id()];
2234  }
2235 
2236  if (precondition_constraint_operator)
2237  column_sums[j] = total_scaling;
2238 
2239  libmesh_error_msg_if
2240  (!total_entries,
2241  "Empty column " << j <<
2242  " found in constraint operator matrix");
2243 
2244  // If we have *cancellation* here then we can end up dividing by
2245  // zero; try just evenly scaling across all constrained node
2246  // points instead.
2247  if (total_scaling > TOLERANCE)
2248  newpt /= total_scaling;
2249  else
2250  newpt /= total_entries;
2251 
2252  Node *n = this->add_point(newpt);
2253  std::unique_ptr<Elem> elem = Elem::build(NODEELEM);
2254  elem->set_node(0, n);
2255  elem->subdomain_id() = new_sbd_id;
2256 
2257  Elem * added_elem = this->add_elem(std::move(elem));
2258  this->_elem_dims.insert(0);
2259  this->_elem_default_orders.insert(added_elem->default_order());
2260  this->_supported_nodal_order =
2261  static_cast<Order>
2262  (std::min(static_cast<int>(this->_supported_nodal_order),
2263  static_cast<int>(added_elem->supported_nodal_order())));
2264  this->_mesh_subdomains.insert(new_sbd_id);
2265  node_to_elem_ptrs.emplace(n, std::make_pair(added_elem->id(), 0));
2266  existing_unconstrained_columns.emplace(j,n->id());
2267 
2268  // Repartition the new objects *after* adding them, so a
2269  // DistributedMesh doesn't get confused and think you're not
2270  // adding them on all processors at once.
2271  int n_pids = 0;
2272  for (auto [pid, count] : pids)
2273  if (count >= n_pids)
2274  {
2275  n_pids = count;
2276  added_elem->processor_id() = pid;
2277  n->processor_id() = pid;
2278  }
2279  }
2280 
2281  // Calculate constraint rows in an indexed form that's easy for us
2282  // to allgather
2283  std::unordered_map<dof_id_type,
2284  std::vector<std::pair<std::pair<dof_id_type, unsigned int>,Real>>>
2285  indexed_constraint_rows;
2286 
2287  for (auto i : make_range(constraint_operator.row_start(),
2288  constraint_operator.row_stop()))
2289  {
2290  if (existing_unconstrained_nodes.count(i))
2291  continue;
2292 
2293  std::vector<numeric_index_type> indices;
2294  std::vector<T> values;
2295 
2296  constraint_operator.get_row(i, indices, values);
2297 
2298  std::vector<std::pair<std::pair<dof_id_type, unsigned int>, Real>> constraint_row;
2299 
2300  for (auto jj : index_range(indices))
2301  {
2302  const dof_id_type node_id =
2303  existing_unconstrained_columns[indices[jj]];
2304 
2305  Node & constraining_node = this->node_ref(node_id);
2306 
2307  libmesh_assert(node_to_elem_ptrs.count(&constraining_node));
2308 
2309  auto p = node_to_elem_ptrs[&constraining_node];
2310 
2311  Real coef = libmesh_real(values[jj]);
2312  libmesh_assert_equal_to(coef, values[jj]);
2313 
2314  // If we're preconditioning and we created a nodeelem then
2315  // we can scale the meaning of that nodeelem's value to give
2316  // us a better-conditioned matrix after the constraints are
2317  // applied.
2318  if (precondition_constraint_operator)
2319  if (auto sum_it = column_sums.find(indices[jj]);
2320  sum_it != column_sums.end())
2321  {
2322  const Real scaling = sum_it->second;
2323 
2324  if (scaling > TOLERANCE)
2325  coef /= scaling;
2326  }
2327 
2328  constraint_row.emplace_back(std::make_pair(p, coef));
2329  }
2330 
2331  indexed_constraint_rows.emplace(i, std::move(constraint_row));
2332  }
2333 
2334  this->comm().set_union(indexed_constraint_rows);
2335 
2336  // Add constraint rows as mesh constraint rows
2337  for (auto & [node_id, indexed_row] : indexed_constraint_rows)
2338  {
2339  Node * constrained_node = this->node_ptr(node_id);
2340 
2341  constraint_rows_mapped_type constraint_row;
2342 
2343  for (auto [p, coef] : indexed_row)
2344  {
2345  const Elem * elem = this->elem_ptr(p.first);
2346  constraint_row.emplace_back
2347  (std::make_pair(std::make_pair(elem, p.second), coef));
2348  }
2349 
2350  this->_constraint_rows.emplace(constrained_node,
2351  std::move(constraint_row));
2352  }
2353 }
2354 
2355 
2356 void MeshBase::print_constraint_rows(std::ostream & os,
2357  bool print_nonlocal) const
2358 {
2359  parallel_object_only();
2360 
2361  std::string local_constraints =
2362  this->get_local_constraints(print_nonlocal);
2363 
2364  if (this->processor_id())
2365  {
2366  this->comm().send(0, local_constraints);
2367  }
2368  else
2369  {
2370  os << "Processor 0:\n";
2371  os << local_constraints;
2372 
2373  for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
2374  {
2375  this->comm().receive(p, local_constraints);
2376  os << "Processor " << p << ":\n";
2377  os << local_constraints;
2378  }
2379  }
2380 }
2381 
2382 
2383 
2384 std::string MeshBase::get_local_constraints(bool print_nonlocal) const
2385 {
2386  std::ostringstream os;
2387 
2388  if (print_nonlocal)
2389  os << "All ";
2390  else
2391  os << "Local ";
2392 
2393  os << "Mesh Constraint Rows:"
2394  << std::endl;
2395 
2396  for (const auto & [node, row] : _constraint_rows)
2397  {
2398  const bool local = (node->processor_id() == this->processor_id());
2399 
2400  // Skip non-local dofs if requested
2401  if (!print_nonlocal && !local)
2402  continue;
2403 
2404  os << "Constraints for " << (local ? "Local" : "Ghost") << " Node " << node->id()
2405  << ": \t";
2406 
2407  for (const auto & [elem_and_node, coef] : row)
2408  os << " ((" << elem_and_node.first->id() << ',' << elem_and_node.second << "), " << coef << ")\t";
2409 
2410  os << std::endl;
2411  }
2412 
2413  return os.str();
2414 }
2415 
2416 
2417 
2418 
2419 // Explicit instantiations for our template function
2420 template LIBMESH_EXPORT void
2421 MeshBase::copy_constraint_rows(const SparseMatrix<Real> & constraint_operator,
2422  bool precondition_constraint_operator);
2423 
2424 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
2425 template LIBMESH_EXPORT void
2426 MeshBase::copy_constraint_rows(const SparseMatrix<Complex> & constraint_operator,
2427  bool precondition_constraint_operator);
2428 #endif
2429 
2430 
2431 } // 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
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:1995
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:1950
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:1703
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:1797
unsigned int n_threads()
Definition: libmesh_base.h:96
Real get_point_locator_close_to_point_tol() const
Definition: mesh_base.C:1933
dof_id_type n_elem_on_proc(const processor_id_type proc) const
Definition: mesh_base.C:1036
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:2033
std::unique_ptr< PointLocatorBase > sub_point_locator() const
Definition: mesh_base.C:1638
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:1811
dof_id_type n_active_elem_on_proc(const processor_id_type proc) const
Definition: mesh_base.C:1049
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:1932
std::vector< std::pair< std::pair< const Elem *, unsigned int >, Real > > constraint_rows_mapped_type
Definition: mesh_base.h:1697
subdomain_id_type n_local_subdomains() const
Definition: mesh_base.C:1011
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:2004
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:2019
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:1685
void remove_ghosting_functor(GhostingFunctor &ghosting_functor)
Removes a functor which was previously added to the set of ghosting functors.
Definition: mesh_base.C:948
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:2104
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:2063
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:1938
unsigned char _spatial_dimension
The "spatial dimension" of the Mesh.
Definition: mesh_base.h:2027
std::unique_ptr< BoundaryInfo > boundary_info
This class holds the boundary information.
Definition: mesh_base.h:1830
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:1963
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:2020
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:2045
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:2021
dof_id_type n_constraint_rows() const
Definition: mesh_base.C:2044
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:2014
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:1712
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.
std::set< GhostingFunctor * > _ghosting_functors
The list of all GhostingFunctor objects to be used when distributing a DistributedMesh.
Definition: mesh_base.h:2086
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:1925
unsigned int _n_parts
The number of partitions the mesh has.
Definition: mesh_base.h:1878
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:1891
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:1694
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:1058
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:1884
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:1562
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:1733
std::unique_ptr< Partitioner > _partitioner
A partitioner to use at each prepare_for_use().
Definition: mesh_base.h:1919
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:1671
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:959
std::set< Order > _elem_default_orders
We cache the (default) order of the geometric elements present in the mesh.
Definition: mesh_base.h:1984
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:1955
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:1958
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:1990
std::set< unsigned char > _elem_dims
We cache the dimension of the elements present in the mesh.
Definition: mesh_base.h:1977
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:1692
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:1023
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:2384
void size_node_extra_integers()
Size extra-integer arrays of all nodes in the mesh.
Definition: mesh_base.C:1949
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:2092
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:1345
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:2051
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:1618
std::set< GhostingFunctor * >::const_iterator ghosting_functors_begin() const
Beginning of range of ghosting functors.
Definition: mesh_base.h:1291
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:2039
subdomain_id_type n_subdomains() const
Definition: mesh_base.C:997
std::string get_info(const unsigned int verbosity=0, const bool global=true) const
Definition: mesh_base.C:1082
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:1969
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:1613
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:1698
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:2356
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:1896
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:1678
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:1608
dof_id_type n_active_sub_elem() const
Same as n_sub_elem(), but only counts active elements.
Definition: mesh_base.C:1070
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:2077
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:1943
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:1725
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:1905
virtual void redistribute()
Redistribute elements between processors.
Definition: mesh_base.C:985
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:2110
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:1911
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:1970
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:1940
virtual numeric_index_type n() const =0
std::set< GhostingFunctor * >::const_iterator ghosting_functors_end() const
End of range of ghosting functors.
Definition: mesh_base.h:1297
void add_ghosting_functor(GhostingFunctor &ghosting_functor)
Adds a functor which can specify ghosting requirements for use on distributed meshes.
Definition: mesh_base.h:1267
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:1919
void set_union(T &data, const unsigned int root_id) const