libMesh
unstructured_mesh.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 it 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 // Local includes
21 #include "libmesh/boundary_info.h"
22 #include "libmesh/ghosting_functor.h"
23 #include "libmesh/ghost_point_neighbors.h"
24 #include "libmesh/unstructured_mesh.h"
25 #include "libmesh/libmesh_logging.h"
26 #include "libmesh/elem.h"
27 #include "libmesh/mesh_tools.h" // For n_levels
28 #include "libmesh/parallel.h"
29 #include "libmesh/remote_elem.h"
30 #include "libmesh/namebased_io.h"
31 #include "libmesh/partitioner.h"
32 #include "libmesh/enum_order.h"
33 #include "libmesh/mesh_communication.h"
34 #include "libmesh/enum_to_string.h"
35 #include "libmesh/mesh_serializer.h"
36 #include "libmesh/utility.h"
37 
38 #ifdef LIBMESH_HAVE_NANOFLANN
39 #include "libmesh/nanoflann.hpp"
40 #endif
41 
42 // C++ includes
43 #include <algorithm> // std::all_of
44 #include <fstream>
45 #include <iomanip>
46 #include <map>
47 #include <sstream>
48 #include <unordered_map>
49 
50 namespace {
51 
52 using namespace libMesh;
53 
54 // Helper functions for all_second_order, all_complete_order
55 
56 std::map<std::vector<dof_id_type>, Node *>::iterator
57 map_hi_order_node(unsigned int hon,
58  const Elem & hi_elem,
59  std::map<std::vector<dof_id_type>, Node *> & adj_vertices_to_ho_nodes)
60 {
61  /*
62  * form a vector that will hold the node id's of
63  * the vertices that are adjacent to the nth
64  * higher-order node.
65  */
66  const unsigned int n_adjacent_vertices =
68 
69  std::vector<dof_id_type> adjacent_vertices_ids(n_adjacent_vertices);
70 
71  for (unsigned int v=0; v<n_adjacent_vertices; v++)
72  adjacent_vertices_ids[v] =
73  hi_elem.node_id( hi_elem.second_order_adjacent_vertex(hon,v) );
74 
75  /*
76  * \p adjacent_vertices_ids is now in order of the current
77  * side. sort it, so that comparisons with the
78  * \p adjacent_vertices_ids created through other elements'
79  * sides can match
80  */
81  std::sort(adjacent_vertices_ids.begin(),
82  adjacent_vertices_ids.end());
83 
84  // Does this set of vertices already have a mid-node added? If not
85  // we'll want to add it.
86  return adj_vertices_to_ho_nodes.try_emplace(adjacent_vertices_ids, nullptr).first;
87 }
88 
89 void transfer_elem(Elem & lo_elem,
90  std::unique_ptr<Elem> hi_elem,
91 #ifdef LIBMESH_ENABLE_UNIQUE_ID
92  unique_id_type max_unique_id,
93  unique_id_type max_new_nodes_per_elem,
94 #endif
96  std::map<std::vector<dof_id_type>, Node *> & adj_vertices_to_ho_nodes,
97  std::unordered_map<Elem *, std::vector<Elem *>> & exterior_children_of)
98 {
99  libmesh_assert_equal_to (lo_elem.n_vertices(), hi_elem->n_vertices());
100 
101  const processor_id_type my_pid = mesh.processor_id();
102  const processor_id_type lo_pid = lo_elem.processor_id();
103 
104  /*
105  * Now handle the additional higher-order nodes. This
106  * is simply handled through a map that remembers
107  * the already-added nodes. This map maps the global
108  * ids of the vertices (that uniquely define this
109  * higher-order node) to the new node.
110  * Notation: hon = high-order node
111  */
112  const unsigned int hon_begin = lo_elem.n_nodes();
113  const unsigned int hon_end = hi_elem->n_nodes();
114 
115  for (unsigned int hon=hon_begin; hon<hon_end; hon++)
116  {
117  auto pos = map_hi_order_node(hon, *hi_elem, adj_vertices_to_ho_nodes);
118 
119  // no, not added yet
120  if (!pos->second)
121  {
122  const auto & adjacent_vertices_ids = pos->first;
123 
124  /*
125  * for this set of vertices, there is no
126  * second_order node yet. Add it.
127  *
128  * compute the location of the new node as
129  * the average over the adjacent vertices.
130  */
131  Point new_location = 0;
132  for (dof_id_type vertex_id : adjacent_vertices_ids)
133  new_location += mesh.point(vertex_id);
134 
135  new_location /= static_cast<Real>(adjacent_vertices_ids.size());
136 
137  /* Add the new point to the mesh.
138  *
139  * If we are on a serialized mesh, then we're doing this
140  * all in sync, and the node processor_id will be
141  * consistent between processors.
142  *
143  * If we are on a distributed mesh, we can fix
144  * inconsistent processor ids later, but only if every
145  * processor gives new nodes a *locally* consistent
146  * processor id, so we'll give the new node the
147  * processor id of an adjacent element for now and then
148  * we'll update that later if appropriate.
149  */
150  Node * hi_node = mesh.add_point
151  (new_location, DofObject::invalid_id, lo_pid);
152 
153  /* Come up with a unique unique_id for a potentially new
154  * node. On a distributed mesh we don't yet know what
155  * processor_id will definitely own it, so we can't let
156  * the pid determine the unique_id. But we're not
157  * adding unpartitioned nodes in sync, so we can't let
158  * the mesh autodetermine a unique_id for a new
159  * unpartitioned node either. So we have to pick unique
160  * unique_id values manually.
161  *
162  * We don't have to pick the *same* unique_id value as
163  * will be picked on other processors, though; we'll
164  * sync up each node later. We just need to make sure
165  * we don't duplicate any unique_id that might be chosen
166  * by the same process elsewhere.
167  */
168 #ifdef LIBMESH_ENABLE_UNIQUE_ID
169  unique_id_type new_unique_id = max_unique_id +
170  max_new_nodes_per_elem * lo_elem.id() +
171  hon - hon_begin;
172 
173  hi_node->set_unique_id(new_unique_id);
174 #endif
175 
176  /*
177  * insert the new node with its defining vertex
178  * set into the map, and relocate pos to this
179  * new entry, so that the hi_elem can use
180  * \p pos for inserting the node
181  */
182  pos->second = hi_node;
183 
184  hi_elem->set_node(hon, hi_node);
185  }
186  // yes, already added.
187  else
188  {
189  Node * hi_node = pos->second;
190  libmesh_assert(hi_node);
191  libmesh_assert_equal_to(mesh.node_ptr(hi_node->id()), hi_node);
192 
193  hi_elem->set_node(hon, hi_node);
194 
195  // We need to ensure that the processor who should own a
196  // node *knows* they own the node. And because
197  // Node::choose_processor_id() may depend on Node id,
198  // which may not yet be authoritative, we still have to
199  // use a dumb-but-id-independent partitioning heuristic.
200  processor_id_type chosen_pid =
201  std::min (hi_node->processor_id(), lo_pid);
202 
203  // Plus, if we just discovered that we own this node,
204  // then on a distributed mesh we need to make sure to
205  // give it a valid id, not just a placeholder id!
206  if (!mesh.is_replicated() &&
207  hi_node->processor_id() != my_pid &&
208  chosen_pid == my_pid)
209  mesh.own_node(*hi_node);
210 
211  hi_node->processor_id() = chosen_pid;
212  }
213  }
214 
215  /*
216  * find_neighbors relies on remote_elem neighbor links being
217  * properly maintained. Our own code here relies on ordinary
218  * neighbor links being properly maintained, so let's just keep
219  * everything up to date.
220  */
221  for (auto s : lo_elem.side_index_range())
222  {
223  Elem * neigh = lo_elem.neighbor_ptr(s);
224  if (!neigh)
225  continue;
226 
227  if (neigh != remote_elem)
228  {
229  // We don't support AMR even outside our own range yet.
230  libmesh_assert_equal_to (neigh->level(), 0);
231 
232  const unsigned int ns = neigh->which_neighbor_am_i(&lo_elem);
233  libmesh_assert_not_equal_to(ns, libMesh::invalid_uint);
234 
235  neigh->set_neighbor(ns, hi_elem.get());
236  }
237 
238  hi_elem->set_neighbor(s, neigh);
239  }
240 
247  Elem * interior_p = lo_elem.interior_parent();
248  if (interior_p)
249  hi_elem->set_interior_parent(interior_p);
250 
251  if (auto parent_exterior_it = exterior_children_of.find(interior_p);
252  parent_exterior_it != exterior_children_of.end())
253  {
254  auto & exteriors = parent_exterior_it->second;
255  for (std::size_t i : index_range(exteriors))
256  if (exteriors[i] == &lo_elem)
257  {
258  exteriors[i] = hi_elem.get();
259  break;
260  }
261  }
262 
267  if (auto exterior_it = exterior_children_of.find(&lo_elem);
268  exterior_it != exterior_children_of.end())
269  {
270  for (Elem * exterior_elem : exterior_it->second)
271  {
272  libmesh_assert(exterior_elem->interior_parent() == &lo_elem);
273  exterior_elem->set_interior_parent(hi_elem.get());
274  }
275  }
276 
289  (mesh.get_boundary_info(), &lo_elem, hi_elem.get());
290 
291  /*
292  * The new second-order element is ready.
293  * Inserting it into the mesh will replace and delete
294  * the first-order element.
295  */
296  hi_elem->set_id(lo_elem.id());
297 #ifdef LIBMESH_ENABLE_UNIQUE_ID
298  hi_elem->set_unique_id(lo_elem.unique_id());
299 #endif
300 
301  const unsigned int nei = lo_elem.n_extra_integers();
302  hi_elem->add_extra_integers(nei);
303  for (unsigned int i=0; i != nei; ++i)
304  hi_elem->set_extra_integer(i, lo_elem.get_extra_integer(i));
305 
306  hi_elem->inherit_data_from(lo_elem);
307 
308  mesh.insert_elem(std::move(hi_elem));
309 }
310 
311 
312 template <typename ElemTypeConverter>
313 void
314 all_increased_order_range (UnstructuredMesh & mesh,
316  const unsigned int max_new_nodes_per_elem,
317  const ElemTypeConverter & elem_type_converter)
318 {
319  // This function must be run on all processors at once
320  timpi_parallel_only(mesh.comm());
321 
322  /*
323  * The maximum number of new higher-order nodes we might be adding,
324  * for use when picking unique unique_id values later. This variable
325  * is not used unless unique ids are enabled, so libmesh_ignore() it
326  * to avoid warnings in that case.
327  */
328  libmesh_ignore(max_new_nodes_per_elem);
329 
330  /*
331  * The mesh should at least be consistent enough for us to add new
332  * nodes consistently.
333  */
336 
337  /*
338  * If the mesh is empty then we have nothing to do
339  */
340  if (!mesh.n_elem())
341  return;
342 
343  // If every element in the range _on every proc_ is already of the
344  // requested higher order then we have nothing to do. However, if
345  // any proc has some lower-order elements in the range, then _all_
346  // processors need to continue this function because it is
347  // parallel_only().
348  //
349  // Note: std::all_of() returns true for an empty range, which can
350  // happen for example in the DistributedMesh case when there are
351  // more processors than elements. In the case of an empty range we
352  // therefore set already_second_order to true on that proc.
353  auto is_higher_order = [&elem_type_converter](const Elem * elem) {
354  ElemType old_type = elem->type();
355  ElemType new_type = elem_type_converter(old_type);
356  return old_type == new_type;
357  };
358 
359  bool already_higher_order =
360  std::all_of(range.begin(), range.end(), is_higher_order);
361 
362  // Check with other processors and possibly return early
363  mesh.comm().min(already_higher_order);
364  if (already_higher_order)
365  return;
366 
367  /*
368  * this map helps in identifying higher order
369  * nodes. Namely, a higher-order node:
370  * - edge node
371  * - face node
372  * - bubble node
373  * is uniquely defined through a set of adjacent
374  * vertices. This set of adjacent vertices is
375  * used to identify already added higher-order
376  * nodes. We are safe to use node id's since we
377  * make sure that these are correctly numbered.
378  *
379  * We lazily use an ordered map here to avoid having to implement a
380  * good hash for vector<dof_id_type>
381  */
382  std::map<std::vector<dof_id_type>, Node *> adj_vertices_to_ho_nodes;
383 
384  /*
385  * This map helps us reset any interior_parent() values from the
386  * lower order element to its higher order replacement. Unlike with
387  * neighbor pointers, we don't have backlinks here, so we have to
388  * iterate over the mesh to track forward links.
389  */
390  std::unordered_map<Elem *, std::vector<Elem *>> exterior_children_of;
391 
392  /*
393  * max_new_nodes_per_elem is the maximum number of new higher order
394  * nodes we might be adding, for use when picking unique unique_id
395  * values later. This variable is not used unless unique ids are
396  * enabled.
397  */
398 #ifdef LIBMESH_ENABLE_UNIQUE_ID
399  unique_id_type max_unique_id = mesh.parallel_max_unique_id();
400 #endif
401 
408  dof_id_type n_unpartitioned_elem = 0,
409  n_partitioned_elem = 0;
410 
425  auto track_if_necessary = [&adj_vertices_to_ho_nodes,
426  &exterior_children_of,
427  &elem_type_converter](Elem * elem) {
428  if (elem && elem != remote_elem)
429  {
430  if (elem->default_order() != FIRST)
431  for (unsigned int hon : make_range(elem->n_vertices(), elem->n_nodes()))
432  {
433  auto pos = map_hi_order_node(hon, *elem, adj_vertices_to_ho_nodes);
434  pos->second = elem->node_ptr(hon);
435  }
436 
437  const ElemType old_type = elem->type();
438  const ElemType new_type = elem_type_converter(old_type);
439  if (old_type != new_type)
440  exterior_children_of.emplace(elem, std::vector<Elem *>());
441  }
442  };
443 
444  // If we're in the common case then just track everything; otherwise
445  // find point neighbors to track
446  if (range.begin() == mesh.elements_begin() &&
447  range.end() == mesh.elements_end())
448  {
449  for (auto & elem : range)
450  track_if_necessary(elem);
451  }
452  else
453  {
454  GhostingFunctor::map_type point_neighbor_elements;
455 
456  GhostPointNeighbors point_neighbor_finder(mesh);
457  point_neighbor_finder(range.begin(), range.end(),
458  mesh.n_processors(),
459  point_neighbor_elements);
460 
461  for (auto & [elem, coupling_map] : point_neighbor_elements)
462  {
463  libmesh_ignore(coupling_map);
464  track_if_necessary(const_cast<Elem *>(elem));
465  }
466  }
467 
472  for (auto & elem : mesh.element_ptr_range())
473  if (auto exterior_map_it = exterior_children_of.find(elem->interior_parent());
474  exterior_map_it != exterior_children_of.end())
475  exterior_map_it->second.push_back(elem);
476 
483  for (auto & lo_elem : range)
484  {
485  // Now we can skip the elements in the range that are already
486  // higher-order.
487  const ElemType old_type = lo_elem->type();
488  const ElemType new_type = elem_type_converter(old_type);
489 
490  if (old_type == new_type)
491  continue;
492 
493  // this does _not_ work for refined elements
494  libmesh_assert_equal_to (lo_elem->level(), 0);
495 
497  ++n_unpartitioned_elem;
498  else
499  ++n_partitioned_elem;
500 
501  /*
502  * Build the higher-order equivalent; add to
503  * the new_elements list.
504  */
505  auto ho_elem = Elem::build (new_type);
506 
507  libmesh_assert_equal_to (lo_elem->n_vertices(), ho_elem->n_vertices());
508 
509  /*
510  * By definition the initial nodes of the lower and higher order
511  * element are identically numbered. Transfer these.
512  */
513  for (unsigned int v=0, lnn=lo_elem->n_nodes(); v < lnn; v++)
514  ho_elem->set_node(v, lo_elem->node_ptr(v));
515 
516  transfer_elem(*lo_elem, std::move(ho_elem),
517 #ifdef LIBMESH_ENABLE_UNIQUE_ID
518  max_unique_id, max_new_nodes_per_elem,
519 #endif
520  mesh, adj_vertices_to_ho_nodes,
521  exterior_children_of);
522  } // end for (auto & lo_elem : range)
523 
524  // we can clear the map at this point.
525  adj_vertices_to_ho_nodes.clear();
526 
527 #ifdef LIBMESH_ENABLE_UNIQUE_ID
528  const unique_id_type new_max_unique_id = max_unique_id +
529  max_new_nodes_per_elem * mesh.n_elem();
530  mesh.set_next_unique_id(new_max_unique_id);
531 #endif
532 
533  // On a DistributedMesh our ghost node processor ids may be bad,
534  // the ids of nodes touching remote elements may be inconsistent,
535  // unique_ids of newly added non-local nodes remain unset, and our
536  // partitioning of new nodes may not be well balanced.
537  //
538  // make_nodes_parallel_consistent() will fix all this.
539  if (!mesh.is_replicated())
540  {
541  dof_id_type max_unpartitioned_elem = n_unpartitioned_elem;
542  mesh.comm().max(max_unpartitioned_elem);
543  if (max_unpartitioned_elem)
544  {
545  // We'd better be effectively serialized here. In theory we
546  // could support more complicated cases but in practice we
547  // only support "completely partitioned" and/or "serialized"
548  if (!mesh.comm().verify(n_unpartitioned_elem) ||
549  !mesh.comm().verify(n_partitioned_elem) ||
550  !mesh.is_serial())
551  libmesh_not_implemented();
552  }
553  else
554  {
556  }
557  }
558 
559  // renumber nodes, repartition nodes, etc. We may no longer need a
560  // find_neighbors() here since we're keeping neighbor links intact
561  // ourselves, *except* that if we're not already prepared we may
562  // have user code that was expecting this call to prepare neighbors.
563  const bool old_find_neighbors = mesh.allow_find_neighbors();
564  if (mesh.is_prepared())
565  mesh.allow_find_neighbors(false);
567  mesh.allow_find_neighbors(old_find_neighbors);
568 }
569 
570 
571 } // anonymous namespace
572 
573 
574 namespace libMesh
575 {
576 
577 // This class adapts a vector of Nodes (represented by a pair of a Point and a dof_id_type)
578 // for use in a nanoflann KD-Tree
579 
581 {
582 private:
583  const std::vector<std::pair<Point, dof_id_type>> _nodes;
584 
585 public:
586  VectorOfNodesAdaptor(const std::vector<std::pair<Point, dof_id_type>> & nodes) :
587  _nodes(nodes)
588  {}
589 
593  inline size_t kdtree_get_point_count() const { return _nodes.size(); }
594 
600  inline Real kdtree_get_pt(const size_t idx, int dim) const
601  {
602  libmesh_assert_less (idx, _nodes.size());
603  libmesh_assert_less (dim, 3);
604 
605  const Point & p(_nodes[idx].first);
606 
607  if (dim==0) return p(0);
608  if (dim==1) return p(1);
609  return p(2);
610  }
611 
612  /*
613  * Optional bounding-box computation
614  */
615  template <class BBOX>
616  bool kdtree_get_bbox(BBOX & /* bb */) const { return false; }
617 };
618 
619 
620 // ------------------------------------------------------------
621 // UnstructuredMesh class member functions
623  unsigned char d) :
624  MeshBase (comm_in,d)
625 {
627 }
628 
629 
630 
632  MeshBase (other_mesh)
633 {
635 }
636 
637 
638 
640  const bool skip_find_neighbors,
641  dof_id_type element_id_offset,
642  dof_id_type node_id_offset,
644 #ifdef LIBMESH_ENABLE_UNIQUE_ID
645  unique_id_offset
646 #endif
647  ,
648  std::unordered_map<subdomain_id_type, subdomain_id_type> *
649  id_remapping)
650 {
651  LOG_SCOPE("copy_nodes_and_elements()", "UnstructuredMesh");
652 
653  std::pair<std::vector<unsigned int>, std::vector<unsigned int>>
654  extra_int_maps = this->merge_extra_integer_names(other_mesh);
655 
656  const unsigned int n_old_node_ints = extra_int_maps.second.size(),
657  n_new_node_ints = _node_integer_names.size(),
658  n_old_elem_ints = extra_int_maps.first.size(),
659  n_new_elem_ints = _elem_integer_names.size();
660 
661  // If we are partitioned into fewer parts than the incoming mesh has
662  // processors to handle, then we need to "wrap" the other Mesh's
663  // processor ids to fit within our range. This can happen, for
664  // example, while stitching meshes with small numbers of elements in
665  // parallel...
666  bool wrap_proc_ids = (this->n_processors() <
667  other_mesh.n_partitions());
668 
669  // We're assuming the other mesh has proper element number ordering,
670  // so that we add parents before their children, and that the other
671  // mesh is consistently partitioned.
672 #ifdef DEBUG
674  MeshTools::libmesh_assert_valid_procids<Node>(other_mesh);
675 #endif
676 
677  //Copy in Nodes
678  {
679  //Preallocate Memory if necessary
680  this->reserve_nodes(other_mesh.n_nodes());
681 
682  for (const auto & oldn : other_mesh.node_ptr_range())
683  {
684  processor_id_type added_pid = cast_int<processor_id_type>
685  (wrap_proc_ids ? oldn->processor_id() % this->n_processors() : oldn->processor_id());
686 
687  // Add new nodes in old node Point locations
688  Node * newn =
689  this->add_point(*oldn,
690  oldn->id() + node_id_offset,
691  added_pid);
692 
693  newn->add_extra_integers(n_new_node_ints);
694  for (unsigned int i = 0; i != n_old_node_ints; ++i)
695  newn->set_extra_integer(extra_int_maps.second[i],
696  oldn->get_extra_integer(i));
697 
698 #ifdef LIBMESH_ENABLE_UNIQUE_ID
699  newn->set_unique_id(oldn->unique_id() + unique_id_offset);
700 #endif
701  }
702  }
703 
704  //Copy in Elements
705  {
706  //Preallocate Memory if necessary
707  this->reserve_elem(other_mesh.n_elem());
708 
709  // Declare a map linking old and new elements, needed to copy the neighbor lists
710  typedef std::unordered_map<const Elem *, Elem *> map_type;
711  map_type old_elems_to_new_elems, ip_map;
712 
713  // Loop over the elements
714  for (const auto & old : other_mesh.element_ptr_range())
715  {
716  // Build a new element
717  Elem * newparent = old->parent() ?
718  this->elem_ptr(old->parent()->id() + element_id_offset) :
719  nullptr;
720  auto el = old->disconnected_clone();
721  el->set_parent(newparent);
722 
723  subdomain_id_type sbd_id = old->subdomain_id();
724  if (id_remapping)
725  {
726  auto remapping_it = id_remapping->find(sbd_id);
727  if (remapping_it != id_remapping->end())
728  sbd_id = remapping_it->second;
729  }
730  el->subdomain_id() = sbd_id;
731 
732  // Hold off on trying to set the interior parent because we may actually
733  // add lower dimensional elements before their interior parents
734  if (el->dim() < other_mesh.mesh_dimension() && old->interior_parent())
735  ip_map[old] = el.get();
736 
737 #ifdef LIBMESH_ENABLE_AMR
738  if (old->has_children())
739  for (unsigned int c = 0, nc = old->n_children(); c != nc; ++c)
740  if (old->child_ptr(c) == remote_elem)
741  el->add_child(const_cast<RemoteElem *>(remote_elem), c);
742 
743  //Create the parent's child pointers if necessary
744  if (newparent)
745  {
746  unsigned int oldc = old->parent()->which_child_am_i(old);
747  newparent->add_child(el.get(), oldc);
748  }
749 
750  // Copy the refinement flags
751  el->set_refinement_flag(old->refinement_flag());
752 
753  // Use hack_p_level since we may not have sibling elements
754  // added yet
755  el->hack_p_level(old->p_level());
756 
757  el->set_p_refinement_flag(old->p_refinement_flag());
758 #endif // #ifdef LIBMESH_ENABLE_AMR
759 
760  //Assign all the nodes
761  for (auto i : el->node_index_range())
762  el->set_node(i,
763  this->node_ptr(old->node_id(i) + node_id_offset));
764 
765  // And start it off with the same processor id (mod _n_parts).
766  el->processor_id() = cast_int<processor_id_type>
767  (wrap_proc_ids ? old->processor_id() % this->n_processors() : old->processor_id());
768 
769  // Give it the same element and unique ids
770  el->set_id(old->id() + element_id_offset);
771 
772  el->add_extra_integers(n_new_elem_ints);
773  for (unsigned int i = 0; i != n_old_elem_ints; ++i)
774  el->set_extra_integer(extra_int_maps.first[i],
775  old->get_extra_integer(i));
776 
777 #ifdef LIBMESH_ENABLE_UNIQUE_ID
778  el->set_unique_id(old->unique_id() + unique_id_offset);
779 #endif
780 
781  //Hold onto it
782  if (!skip_find_neighbors)
783  {
784  for (auto s : old->side_index_range())
785  if (old->neighbor_ptr(s) == remote_elem)
786  el->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
787  this->add_elem(std::move(el));
788  }
789  else
790  {
791  Elem * new_el = this->add_elem(std::move(el));
792  old_elems_to_new_elems[old] = new_el;
793  }
794  }
795 
796  for (auto & elem_pair : ip_map)
797  elem_pair.second->set_interior_parent(
798  this->elem_ptr(elem_pair.first->interior_parent()->id() + element_id_offset));
799 
800  // Loop (again) over the elements to fill in the neighbors
801  if (skip_find_neighbors)
802  {
803  old_elems_to_new_elems[remote_elem] = const_cast<RemoteElem*>(remote_elem);
804 
805  for (const auto & old_elem : other_mesh.element_ptr_range())
806  {
807  Elem * new_elem = old_elems_to_new_elems[old_elem];
808  for (auto s : old_elem->side_index_range())
809  {
810  const Elem * old_neighbor = old_elem->neighbor_ptr(s);
811  Elem * new_neighbor = old_elems_to_new_elems[old_neighbor];
812  new_elem->set_neighbor(s, new_neighbor);
813  }
814  }
815  }
816  }
817 
818 #ifdef LIBMESH_ENABLE_UNIQUE_ID
819  // We set the unique ids of nodes after adding them to the mesh such that our value of
820  // _next_unique_id may be wrong. So we amend that here
821  this->set_next_unique_id(other_mesh.parallel_max_unique_id() + unique_id_offset + 1);
822 #endif
823 
824  // Finally, partially prepare the new Mesh for use.
825  // This is for backwards compatibility, so we don't want to prepare
826  // everything.
827  //
828  // Keep the same numbering and partitioning and distribution status
829  // for now, but save our original policies to restore later.
830  const bool allowed_renumbering = this->allow_renumbering();
831  const bool allowed_find_neighbors = this->allow_find_neighbors();
832  const bool allowed_elem_removal = this->allow_remote_element_removal();
833  this->allow_renumbering(false);
834  this->allow_remote_element_removal(false);
835  this->allow_find_neighbors(!skip_find_neighbors);
836 
837  // We should generally be able to skip *all* partitioning here
838  // because we're only adding one already-consistent mesh to another.
839  const bool skipped_partitioning = this->skip_partitioning();
840  this->skip_partitioning(true);
841 
842  const bool was_prepared = this->is_prepared();
843  this->prepare_for_use();
844 
845  //But in the long term, don't change our policies.
846  this->allow_find_neighbors(allowed_find_neighbors);
847  this->allow_renumbering(allowed_renumbering);
848  this->allow_remote_element_removal(allowed_elem_removal);
849  this->skip_partitioning(skipped_partitioning);
850 
851  // That prepare_for_use() call marked us as prepared, but we
852  // specifically avoided some important preparation, so we might not
853  // actually be prepared now.
854  if (skip_find_neighbors ||
855  !was_prepared || !other_mesh.is_prepared())
856  this->set_isnt_prepared();
857 }
858 
859 
860 
862 {
863  // this->clear (); // Nothing to clear at this level
864 
865  libmesh_exceptionless_assert (!libMesh::closed());
866 }
867 
868 
869 
870 
871 
872 void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
873  const bool reset_current_list)
874 {
875  // We might actually want to run this on an empty mesh
876  // (e.g. the boundary mesh for a nonexistent bcid!)
877  // libmesh_assert_not_equal_to (this->n_nodes(), 0);
878  // libmesh_assert_not_equal_to (this->n_elem(), 0);
879 
880  // This function must be run on all processors at once
881  parallel_object_only();
882 
883  LOG_SCOPE("find_neighbors()", "Mesh");
884 
885  //TODO:[BSK] This should be removed later?!
886  if (reset_current_list)
887  for (const auto & e : this->element_ptr_range())
888  for (auto s : e->side_index_range())
889  if (e->neighbor_ptr(s) != remote_elem || reset_remote_elements)
890  e->set_neighbor(s, nullptr);
891 
892  // Find neighboring elements by first finding elements
893  // with identical side keys and then check to see if they
894  // are neighbors
895  {
896  // data structures -- Use the hash_multimap if available
897  typedef dof_id_type key_type;
898  typedef std::pair<Elem *, unsigned char> val_type;
899  typedef std::unordered_multimap<key_type, val_type> map_type;
900 
901  // A map from side keys to corresponding elements & side numbers
902  map_type side_to_elem_map;
903 
904  // Pull objects out of the loop to reduce heap operations
905  std::unique_ptr<Elem> my_side, their_side;
906 
907  for (const auto & element : this->element_ptr_range())
908  {
909  for (auto ms : element->side_index_range())
910  {
911  next_side:
912  // If we haven't yet found a neighbor on this side, try.
913  // Even if we think our neighbor is remote, that
914  // information may be out of date.
915  if (element->neighbor_ptr(ms) == nullptr ||
916  element->neighbor_ptr(ms) == remote_elem)
917  {
918  // Get the key for the side of this element. Use the
919  // low_order_key so we can find neighbors in
920  // mixed-order meshes if necessary.
921  const dof_id_type key = element->low_order_key(ms);
922 
923  // Look for elements that have an identical side key
924  auto bounds = side_to_elem_map.equal_range(key);
925 
926  // May be multiple keys, check all the possible
927  // elements which _might_ be neighbors.
928  if (bounds.first != bounds.second)
929  {
930  // Get the side for this element
931  element->side_ptr(my_side, ms);
932 
933  // Look at all the entries with an equivalent key
934  while (bounds.first != bounds.second)
935  {
936  // Get the potential element
937  Elem * neighbor = bounds.first->second.first;
938 
939  // Get the side for the neighboring element
940  const unsigned int ns = bounds.first->second.second;
941  neighbor->side_ptr(their_side, ns);
942  //libmesh_assert(my_side.get());
943  //libmesh_assert(their_side.get());
944 
945  // If found a match with my side
946  //
947  // In 1D, since parents and children have an
948  // equal side (i.e. a node) we need to check
949  // for matching level() to avoid setting our
950  // neighbor pointer to any of our neighbor's
951  // descendants.
952  if ((*my_side == *their_side) &&
953  (element->level() == neighbor->level()))
954  {
955  // So share a side. Is this a mixed pair
956  // of subactive and active/ancestor
957  // elements?
958  // If not, then we're neighbors.
959  // If so, then the subactive's neighbor is
960 
961  if (element->subactive() ==
962  neighbor->subactive())
963  {
964  // an element is only subactive if it has
965  // been coarsened but not deleted
966  element->set_neighbor (ms,neighbor);
967  neighbor->set_neighbor(ns,element);
968  }
969  else if (element->subactive())
970  {
971  element->set_neighbor(ms,neighbor);
972  }
973  else if (neighbor->subactive())
974  {
975  neighbor->set_neighbor(ns,element);
976  }
977  side_to_elem_map.erase (bounds.first);
978 
979  // get out of this nested crap
980  goto next_side;
981  }
982 
983  ++bounds.first;
984  }
985  }
986 
987  // didn't find a match...
988  // Build the map entry for this element
989  side_to_elem_map.emplace
990  (key, std::make_pair(element, cast_int<unsigned char>(ms)));
991  }
992  }
993  }
994  }
995 
996 #ifdef LIBMESH_ENABLE_AMR
997 
1025  const unsigned int n_levels = MeshTools::n_levels(*this);
1026  for (unsigned int level = 1; level < n_levels; ++level)
1027  {
1028  for (auto & current_elem : as_range(level_elements_begin(level),
1029  level_elements_end(level)))
1030  {
1031  libmesh_assert(current_elem);
1032  Elem * parent = current_elem->parent();
1033  libmesh_assert(parent);
1034  const unsigned int my_child_num = parent->which_child_am_i(current_elem);
1035 
1036  for (auto s : current_elem->side_index_range())
1037  {
1038  if (current_elem->neighbor_ptr(s) == nullptr ||
1039  (current_elem->neighbor_ptr(s) == remote_elem &&
1040  parent->is_child_on_side(my_child_num, s)))
1041  {
1042  Elem * neigh = parent->neighbor_ptr(s);
1043 
1044  // If neigh was refined and had non-subactive children
1045  // made remote earlier, then our current elem should
1046  // actually have one of those remote children as a
1047  // neighbor
1048  if (neigh &&
1049  (neigh->ancestor() ||
1050  // If neigh has subactive children which should have
1051  // matched as neighbors of the current element but
1052  // did not, then those likewise must be remote
1053  // children.
1054  (current_elem->subactive() && neigh->has_children() &&
1055  (neigh->level()+1) == current_elem->level())))
1056  {
1057 #ifdef DEBUG
1058  // Let's make sure that "had children made remote"
1059  // situation is actually the case
1060  libmesh_assert(neigh->has_children());
1061  bool neigh_has_remote_children = false;
1062  for (auto & child : neigh->child_ref_range())
1063  if (&child == remote_elem)
1064  neigh_has_remote_children = true;
1065  libmesh_assert(neigh_has_remote_children);
1066 
1067  // And let's double-check that we don't have
1068  // a remote_elem neighboring an active local element
1069  if (current_elem->active())
1070  libmesh_assert_not_equal_to (current_elem->processor_id(),
1071  this->processor_id());
1072 #endif // DEBUG
1073  neigh = const_cast<RemoteElem *>(remote_elem);
1074  }
1075  // If neigh and current_elem are more than one level
1076  // apart, figuring out whether we have a remote
1077  // neighbor here becomes much harder.
1078  else if (neigh && (current_elem->subactive() &&
1079  neigh->has_children()))
1080  {
1081  // Find the deepest descendant of neigh which
1082  // we could consider for a neighbor. If we run
1083  // out of neigh children, then that's our
1084  // neighbor. If we find a potential neighbor
1085  // with remote_children and we don't find any
1086  // potential neighbors among its non-remote
1087  // children, then our neighbor must be remote.
1088  while (neigh != remote_elem &&
1089  neigh->has_children())
1090  {
1091  bool found_neigh = false;
1092  for (unsigned int c = 0, nc = neigh->n_children();
1093  !found_neigh && c != nc; ++c)
1094  {
1095  Elem * child = neigh->child_ptr(c);
1096  if (child == remote_elem)
1097  continue;
1098  for (auto ncn : child->neighbor_ptr_range())
1099  {
1100  if (ncn != remote_elem &&
1101  ncn->is_ancestor_of(current_elem))
1102  {
1103  neigh = ncn;
1104  found_neigh = true;
1105  break;
1106  }
1107  }
1108  }
1109  if (!found_neigh)
1110  neigh = const_cast<RemoteElem *>(remote_elem);
1111  }
1112  }
1113  current_elem->set_neighbor(s, neigh);
1114 #ifdef DEBUG
1115  if (neigh != nullptr && neigh != remote_elem)
1116  // We ignore subactive elements here because
1117  // we don't care about neighbors of subactive element.
1118  if ((!neigh->active()) && (!current_elem->subactive()))
1119  {
1120  libMesh::err << "On processor " << this->processor_id()
1121  << std::endl;
1122  libMesh::err << "Bad element ID = " << current_elem->id()
1123  << ", Side " << s << ", Bad neighbor ID = " << neigh->id() << std::endl;
1124  libMesh::err << "Bad element proc_ID = " << current_elem->processor_id()
1125  << ", Bad neighbor proc_ID = " << neigh->processor_id() << std::endl;
1126  libMesh::err << "Bad element size = " << current_elem->hmin()
1127  << ", Bad neighbor size = " << neigh->hmin() << std::endl;
1128  libMesh::err << "Bad element center = " << current_elem->vertex_average()
1129  << ", Bad neighbor center = " << neigh->vertex_average() << std::endl;
1130  libMesh::err << "ERROR: "
1131  << (current_elem->active()?"Active":"Ancestor")
1132  << " Element at level "
1133  << current_elem->level() << std::endl;
1134  libMesh::err << "with "
1135  << (parent->active()?"active":
1136  (parent->subactive()?"subactive":"ancestor"))
1137  << " parent share "
1138  << (neigh->subactive()?"subactive":"ancestor")
1139  << " neighbor at level " << neigh->level()
1140  << std::endl;
1141  NameBasedIO(*this).write ("bad_mesh.gmv");
1142  libmesh_error_msg("Problematic mesh written to bad_mesh.gmv.");
1143  }
1144 #endif // DEBUG
1145  }
1146  }
1147 
1148  // We can skip to the next element if we're full-dimension
1149  // and therefore don't have any interior parents
1150  if (current_elem->dim() >= LIBMESH_DIM)
1151  continue;
1152 
1153  // We have no interior parents unless we can find one later
1154  current_elem->set_interior_parent(nullptr);
1155 
1156  Elem * pip = parent->interior_parent();
1157 
1158  if (!pip)
1159  continue;
1160 
1161  // If there's no interior_parent children, whether due to a
1162  // remote element or a non-conformity, then there's no
1163  // children to search.
1164  if (pip == remote_elem || pip->active())
1165  {
1166  current_elem->set_interior_parent(pip);
1167  continue;
1168  }
1169 
1170  // For node comparisons we'll need a sensible tolerance
1171  Real node_tolerance = current_elem->hmin() * TOLERANCE;
1172 
1173  // Otherwise our interior_parent should be a child of our
1174  // parent's interior_parent.
1175  for (auto & child : pip->child_ref_range())
1176  {
1177  // If we have a remote_elem, that might be our
1178  // interior_parent. We'll set it provisionally now and
1179  // keep trying to find something better.
1180  if (&child == remote_elem)
1181  {
1182  current_elem->set_interior_parent
1183  (const_cast<RemoteElem *>(remote_elem));
1184  continue;
1185  }
1186 
1187  bool child_contains_our_nodes = true;
1188  for (auto & n : current_elem->node_ref_range())
1189  {
1190  bool child_contains_this_node = false;
1191  for (auto & cn : child.node_ref_range())
1192  if (cn.absolute_fuzzy_equals
1193  (n, node_tolerance))
1194  {
1195  child_contains_this_node = true;
1196  break;
1197  }
1198  if (!child_contains_this_node)
1199  {
1200  child_contains_our_nodes = false;
1201  break;
1202  }
1203  }
1204  if (child_contains_our_nodes)
1205  {
1206  current_elem->set_interior_parent(&child);
1207  break;
1208  }
1209  }
1210 
1211  // We should have found *some* interior_parent at this
1212  // point, whether semilocal or remote.
1213  libmesh_assert(current_elem->interior_parent());
1214  }
1215  }
1216 
1217 #endif // AMR
1218 
1219 
1220 #ifdef DEBUG
1222  !reset_remote_elements);
1224 #endif
1225 }
1226 
1227 
1228 
1229 void UnstructuredMesh::read (const std::string & name,
1230  void *,
1231  bool skip_renumber_nodes_and_elements,
1232  bool skip_find_neighbors)
1233 {
1234  // Set the skip_renumber_nodes_and_elements flag on all processors
1235  // if necessary.
1236  // This ensures that renumber_nodes_and_elements is *not* called
1237  // during prepare_for_use() for certain types of mesh files.
1238  // This is required in cases where there is an associated solution
1239  // file which expects a certain ordering of the nodes.
1240  if (name.rfind(".gmv") == name.size() - 4)
1241  this->allow_renumbering(false);
1242 
1243  NameBasedIO(*this).read(name);
1244 
1245  if (skip_renumber_nodes_and_elements)
1246  {
1247  // Use MeshBase::allow_renumbering() yourself instead.
1248  libmesh_deprecated();
1249  this->allow_renumbering(false);
1250  }
1251 
1252  // Done reading the mesh. Now prepare it for use.
1253  const bool old_allow_find_neighbors = this->allow_find_neighbors();
1254  this->allow_find_neighbors(!skip_find_neighbors);
1255  this->prepare_for_use();
1256  this->allow_find_neighbors(old_allow_find_neighbors);
1257 }
1258 
1259 
1260 
1261 void UnstructuredMesh::write (const std::string & name) const
1262 {
1263  LOG_SCOPE("write()", "Mesh");
1264 
1265  NameBasedIO(*this).write(name);
1266 }
1267 
1268 
1269 
1270 void UnstructuredMesh::write (const std::string & name,
1271  const std::vector<Number> & v,
1272  const std::vector<std::string> & vn) const
1273 {
1274  LOG_SCOPE("write()", "Mesh");
1275 
1276  NameBasedIO(*this).write_nodal_data(name, v, vn);
1277 }
1278 
1279 
1280 
1281 
1282 
1284  const processor_id_type pid) const
1285 {
1286 
1287  // Issue a warning if the number the number of processors
1288  // currently available is less that that requested for
1289  // partitioning. This is not necessarily an error since
1290  // you may run on one processor and still partition the
1291  // mesh into several partitions.
1292 #ifdef DEBUG
1293  if (this->n_processors() < pid)
1294  {
1295  libMesh::out << "WARNING: You are creating a "
1296  << "mesh for a processor id (="
1297  << pid
1298  << ") greater than "
1299  << "the number of processors available for "
1300  << "the calculation. (="
1301  << this->n_processors()
1302  << ")."
1303  << std::endl;
1304  }
1305 #endif
1306 
1307  this->create_submesh (pid_mesh,
1308  this->active_pid_elements_begin(pid),
1309  this->active_pid_elements_end(pid));
1310 }
1311 
1312 
1313 
1314 
1315 
1316 
1317 
1319  const const_element_iterator & it,
1320  const const_element_iterator & it_end) const
1321 {
1322  // Just in case the subdomain_mesh already has some information
1323  // in it, get rid of it.
1324  new_mesh.clear();
1325 
1326  // If we're not serial, our submesh isn't either.
1327  // There are no remote elements to delete on an empty mesh, but
1328  // calling the method to do so marks the mesh as parallel.
1329  if (!this->is_serial())
1330  new_mesh.delete_remote_elements();
1331 
1332  // Fail if (*this == new_mesh), we cannot create a submesh inside ourself!
1333  // This may happen if the user accidentally passes the original mesh into
1334  // this function! We will check this by making sure we did not just
1335  // clear ourself.
1336  libmesh_assert_not_equal_to (this->n_nodes(), 0);
1337  libmesh_assert_not_equal_to (this->n_elem(), 0);
1338 
1339  // Container to catch boundary IDs handed back by BoundaryInfo
1340  std::vector<boundary_id_type> bc_ids;
1341 
1342  // Put any extra integers on the new mesh too
1343  new_mesh.merge_extra_integer_names(*this);
1344  const unsigned int n_node_ints = _node_integer_names.size();
1345 
1346  for (const auto & old_elem : as_range(it, it_end))
1347  {
1348  // Add an equivalent element type to the new_mesh.
1349  // disconnected_clone() copies ids, extra element integers, etc.
1350  auto uelem = old_elem->disconnected_clone();
1351  Elem * new_elem = new_mesh.add_elem(std::move(uelem));
1352  libmesh_assert(new_elem);
1353 
1354  // Loop over the nodes on this element.
1355  for (auto n : old_elem->node_index_range())
1356  {
1357  const dof_id_type this_node_id = old_elem->node_id(n);
1358 
1359  // Add this node to the new mesh if it's not there already
1360  if (!new_mesh.query_node_ptr(this_node_id))
1361  {
1362  Node * newn =
1363  new_mesh.add_point (old_elem->point(n),
1364  this_node_id,
1365  old_elem->node_ptr(n)->processor_id());
1366 
1367  newn->add_extra_integers(n_node_ints);
1368  for (unsigned int i = 0; i != n_node_ints; ++i)
1369  newn->set_extra_integer(i, old_elem->node_ptr(n)->get_extra_integer(i));
1370 
1371 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1372  newn->set_unique_id(old_elem->node_ptr(n)->unique_id());
1373 #endif
1374  }
1375 
1376  // Define this element's connectivity on the new mesh
1377  new_elem->set_node(n, new_mesh.node_ptr(this_node_id));
1378  }
1379 
1380  // Maybe add boundary conditions for this element
1381  for (auto s : old_elem->side_index_range())
1382  {
1383  this->get_boundary_info().boundary_ids(old_elem, s, bc_ids);
1384  new_mesh.get_boundary_info().add_side (new_elem, s, bc_ids);
1385  }
1386  } // end loop over elements
1387 
1388  // Prepare the new_mesh for use
1389  new_mesh.prepare_for_use();
1390 }
1391 
1392 
1393 
1394 #ifdef LIBMESH_ENABLE_AMR
1396 {
1397  LOG_SCOPE ("contract()", "Mesh");
1398 
1399  // Flag indicating if this call actually changes the mesh
1400  bool mesh_changed = false;
1401 
1402 #ifdef DEBUG
1403  for (const auto & elem : this->element_ptr_range())
1404  libmesh_assert(elem->active() || elem->subactive() || elem->ancestor());
1405 #endif
1406 
1407  // Loop over the elements.
1408  for (auto & elem : this->element_ptr_range())
1409  {
1410  // Delete all the subactive ones
1411  if (elem->subactive())
1412  {
1413  // No level-0 element should be subactive.
1414  // Note that we CAN'T test elem->level(), as that
1415  // touches elem->parent()->dim(), and elem->parent()
1416  // might have already been deleted!
1417  libmesh_assert(elem->parent());
1418 
1419  // Delete the element
1420  // This just sets a pointer to nullptr, and doesn't
1421  // invalidate any iterators
1422  this->delete_elem(elem);
1423 
1424  // the mesh has certainly changed
1425  mesh_changed = true;
1426  }
1427  else
1428  {
1429  // Compress all the active ones
1430  if (elem->active())
1431  elem->contract();
1432  else
1433  libmesh_assert (elem->ancestor());
1434  }
1435  }
1436 
1437  // Strip any newly-created nullptr voids out of the element array
1439 
1440  // FIXME: Need to understand why deleting subactive children
1441  // invalidates the point locator. For now we will clear it explicitly
1442  this->clear_point_locator();
1443 
1444  // Allow our GhostingFunctor objects to reinit if necessary.
1445  for (auto & gf : as_range(this->ghosting_functors_begin(),
1446  this->ghosting_functors_end()))
1447  {
1448  libmesh_assert(gf);
1449  gf->mesh_reinit();
1450  }
1451 
1452  return mesh_changed;
1453 }
1454 #endif // #ifdef LIBMESH_ENABLE_AMR
1455 
1456 
1457 
1459 {
1460  LOG_SCOPE("all_first_order()", "Mesh");
1461 
1465  std::vector<bool> node_touched_by_me(this->max_node_id(), false);
1466 
1467  // Loop over the high-ordered elements.
1468  // First make sure they _are_ indeed high-order, and then replace
1469  // them with an equivalent first-order element.
1470  for (auto & so_elem : element_ptr_range())
1471  {
1472  libmesh_assert(so_elem);
1473 
1474  /*
1475  * build the first-order equivalent, add to
1476  * the new_elements list.
1477  */
1478  auto lo_elem = Elem::build
1480  (so_elem->type()), so_elem->parent());
1481 
1482  const unsigned short n_sides = so_elem->n_sides();
1483 
1484  for (unsigned short s=0; s != n_sides; ++s)
1485  if (so_elem->neighbor_ptr(s) == remote_elem)
1486  lo_elem->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
1487 
1488 #ifdef LIBMESH_ENABLE_AMR
1489  /*
1490  * Reset the parent links of any child elements
1491  */
1492  if (so_elem->has_children())
1493  for (unsigned int c = 0, nc = so_elem->n_children(); c != nc; ++c)
1494  {
1495  Elem * child = so_elem->child_ptr(c);
1496  if (child != remote_elem)
1497  child->set_parent(lo_elem.get());
1498  lo_elem->add_child(child, c);
1499  }
1500 
1501  /*
1502  * Reset the child link of any parent element
1503  */
1504  if (so_elem->parent())
1505  {
1506  unsigned int c =
1507  so_elem->parent()->which_child_am_i(so_elem);
1508  lo_elem->parent()->replace_child(lo_elem.get(), c);
1509  }
1510 
1511  /*
1512  * Copy as much data to the new element as makes sense
1513  */
1514  lo_elem->set_p_level(so_elem->p_level());
1515  lo_elem->set_refinement_flag(so_elem->refinement_flag());
1516  lo_elem->set_p_refinement_flag(so_elem->p_refinement_flag());
1517 #endif
1518 
1519  libmesh_assert_equal_to (lo_elem->n_vertices(), so_elem->n_vertices());
1520 
1521  /*
1522  * By definition the vertices of the linear and
1523  * second order element are identically numbered.
1524  * transfer these.
1525  */
1526  for (unsigned int v=0, snv=so_elem->n_vertices(); v < snv; v++)
1527  {
1528  lo_elem->set_node(v, so_elem->node_ptr(v));
1529  node_touched_by_me[lo_elem->node_id(v)] = true;
1530  }
1531 
1532  /*
1533  * find_neighbors relies on remote_elem neighbor links being
1534  * properly maintained.
1535  */
1536  for (unsigned short s=0; s != n_sides; s++)
1537  {
1538  if (so_elem->neighbor_ptr(s) == remote_elem)
1539  lo_elem->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
1540  }
1541 
1549  (this->get_boundary_info(), so_elem, lo_elem.get());
1550 
1551  /*
1552  * The new first-order element is ready.
1553  * Inserting it into the mesh will replace and delete
1554  * the second-order element.
1555  */
1556  lo_elem->set_id(so_elem->id());
1557 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1558  lo_elem->set_unique_id(so_elem->unique_id());
1559 #endif
1560 
1561  const unsigned int nei = so_elem->n_extra_integers();
1562  lo_elem->add_extra_integers(nei);
1563  for (unsigned int i=0; i != nei; ++i)
1564  lo_elem->set_extra_integer(i, so_elem->get_extra_integer(i));
1565 
1566  lo_elem->inherit_data_from(*so_elem);
1567 
1568  this->insert_elem(std::move(lo_elem));
1569  }
1570 
1571  // Deleting nodes does not invalidate iterators, so this is safe.
1572  for (const auto & node : this->node_ptr_range())
1573  if (!node_touched_by_me[node->id()])
1574  this->delete_node(node);
1575 
1576  // If crazy people applied boundary info to non-vertices and then
1577  // deleted those non-vertices, we should make sure their boundary id
1578  // caches are correct.
1580 
1581  // On hanging nodes that used to also be second order nodes, we
1582  // might now have an invalid nodal processor_id()
1584 
1585  // delete or renumber nodes if desired
1586  this->prepare_for_use();
1587 }
1588 
1589 
1590 
1591 void
1593  const bool full_ordered)
1594 {
1595  LOG_SCOPE("all_second_order_range()", "Mesh");
1596 
1597  /*
1598  * The maximum number of new second order nodes we might be adding,
1599  * for use when picking unique unique_id values later. This variable
1600  * is not used unless unique ids are enabled.
1601  */
1602  unsigned int max_new_nodes_per_elem;
1603 
1604  /*
1605  * For speed-up of the \p add_point() method, we
1606  * can reserve memory. Guess the number of additional
1607  * nodes based on the element spatial dimensions and the
1608  * total number of nodes in the mesh as an upper bound.
1609  */
1610  switch (this->mesh_dimension())
1611  {
1612  case 1:
1613  /*
1614  * in 1D, there can only be order-increase from Edge2
1615  * to Edge3. Something like 1/2 of n_nodes() have
1616  * to be added
1617  */
1618  max_new_nodes_per_elem = 3 - 2;
1619  this->reserve_nodes(static_cast<unsigned int>
1620  (1.5*static_cast<double>(this->n_nodes())));
1621  break;
1622 
1623  case 2:
1624  /*
1625  * in 2D, either refine from Tri3 to Tri6 (double the nodes)
1626  * or from Quad4 to Quad8 (again, double) or Quad9 (2.25 that much)
1627  */
1628  max_new_nodes_per_elem = 9 - 4;
1629  this->reserve_nodes(static_cast<unsigned int>
1630  (2*static_cast<double>(this->n_nodes())));
1631  break;
1632 
1633 
1634  case 3:
1635  /*
1636  * in 3D, either refine from Tet4 to Tet10 (factor = 2.5) up to
1637  * Hex8 to Hex27 (something > 3). Since in 3D there _are_ already
1638  * quite some nodes, and since we do not want to overburden the memory by
1639  * a too conservative guess, use the lower bound
1640  */
1641  max_new_nodes_per_elem = 27 - 8;
1642  this->reserve_nodes(static_cast<unsigned int>
1643  (2.5*static_cast<double>(this->n_nodes())));
1644  break;
1645 
1646  default:
1647  // Hm?
1648  libmesh_error_msg("Unknown mesh dimension " << this->mesh_dimension());
1649  }
1650 
1651  // All the real work is done in the helper function
1652  all_increased_order_range(*this, range, max_new_nodes_per_elem,
1653  [full_ordered](ElemType t) {
1654  return Elem::second_order_equivalent_type(t, full_ordered);
1655  });
1656 }
1657 
1658 
1659 
1661 {
1662  LOG_SCOPE("all_complete_order()", "Mesh");
1663 
1664  /*
1665  * The maximum number of new higher-order nodes we might be adding,
1666  * for use when picking unique unique_id values later. This variable
1667  * is not used unless unique ids are enabled.
1668  */
1669  unsigned int max_new_nodes_per_elem;
1670 
1671  /*
1672  * for speed-up of the \p add_point() method, we
1673  * can reserve memory. Guess the number of additional
1674  * nodes based on the element spatial dimensions and the
1675  * total number of nodes in the mesh as an upper bound.
1676  */
1677  switch (this->mesh_dimension())
1678  {
1679  case 1:
1680  /*
1681  * in 1D, there can only be order-increase from Edge2
1682  * to Edge3. Something like 1/2 of n_nodes() have
1683  * to be added
1684  */
1685  max_new_nodes_per_elem = 3 - 2;
1686  this->reserve_nodes(static_cast<unsigned int>
1687  (1.5*static_cast<double>(this->n_nodes())));
1688  break;
1689 
1690  case 2:
1691  /*
1692  * in 2D, we typically refine from Tri6 to Tri7 (1.1667 times
1693  * the nodes) but might refine from Quad4 to Quad9
1694  * (2.25 times the nodes)
1695  */
1696  max_new_nodes_per_elem = 9 - 4;
1697  this->reserve_nodes(static_cast<unsigned int>
1698  (2*static_cast<double>(this->n_nodes())));
1699  break;
1700 
1701 
1702  case 3:
1703  /*
1704  * in 3D, we typically refine from Tet10 to Tet14 (factor = 1.4)
1705  * but may go Hex8 to Hex27 (something > 3). Since in 3D there
1706  * _are_ already quite some nodes, and since we do not want to
1707  * overburden the memory by a too conservative guess, use a
1708  * moderate bound
1709  */
1710  max_new_nodes_per_elem = 27 - 8;
1711  this->reserve_nodes(static_cast<unsigned int>
1712  (2.5*static_cast<double>(this->n_nodes())));
1713  break;
1714 
1715  default:
1716  // Hm?
1717  libmesh_error_msg("Unknown mesh dimension " << this->mesh_dimension());
1718  }
1719 
1720  // All the real work is done in the helper function
1721  all_increased_order_range(*this, range, max_new_nodes_per_elem,
1722  [](ElemType t) {
1724  });
1725 }
1726 
1727 
1728 std::size_t
1730  boundary_id_type this_mesh_boundary_id,
1731  boundary_id_type other_mesh_boundary_id,
1732  Real tol,
1733  bool clear_stitched_boundary_ids,
1734  bool verbose,
1735  bool use_binary_search,
1736  bool enforce_all_nodes_match_on_boundaries,
1737  bool merge_boundary_nodes_all_or_nothing,
1738  bool remap_subdomain_ids)
1739 {
1740  LOG_SCOPE("stitch_meshes()", "UnstructuredMesh");
1741  return stitching_helper(&other_mesh,
1742  this_mesh_boundary_id,
1743  other_mesh_boundary_id,
1744  tol,
1745  clear_stitched_boundary_ids,
1746  verbose,
1747  use_binary_search,
1748  enforce_all_nodes_match_on_boundaries,
1749  true,
1750  merge_boundary_nodes_all_or_nothing,
1751  remap_subdomain_ids);
1752 }
1753 
1754 
1755 std::size_t
1757  boundary_id_type boundary_id_2,
1758  Real tol,
1759  bool clear_stitched_boundary_ids,
1760  bool verbose,
1761  bool use_binary_search,
1762  bool enforce_all_nodes_match_on_boundaries,
1763  bool merge_boundary_nodes_all_or_nothing)
1764 
1765 {
1766  return stitching_helper(nullptr,
1767  boundary_id_1,
1768  boundary_id_2,
1769  tol,
1770  clear_stitched_boundary_ids,
1771  verbose,
1772  use_binary_search,
1773  enforce_all_nodes_match_on_boundaries,
1774  true,
1775  merge_boundary_nodes_all_or_nothing,
1776  false);
1777 }
1778 
1779 
1780 std::size_t
1782  boundary_id_type this_mesh_boundary_id,
1783  boundary_id_type other_mesh_boundary_id,
1784  Real tol,
1785  bool clear_stitched_boundary_ids,
1786  bool verbose,
1787  bool use_binary_search,
1788  bool enforce_all_nodes_match_on_boundaries,
1789  bool skip_find_neighbors,
1790  bool merge_boundary_nodes_all_or_nothing,
1791  bool remap_subdomain_ids)
1792 {
1793 #ifdef DEBUG
1794  // We rely on neighbor links here
1796 #endif
1797 
1798  // We can't even afford any unset neighbor links here.
1799  if (!this->is_prepared())
1800  this->find_neighbors();
1801 
1802  // FIXME: make distributed mesh support efficient.
1803  // Yes, we currently suck.
1804  MeshSerializer serialize(*this);
1805 
1806  // *Badly*.
1807  std::unique_ptr<MeshSerializer> serialize_other;
1808  if (other_mesh)
1809  serialize_other = std::make_unique<MeshSerializer>
1810  (*const_cast<MeshBase *>(other_mesh));
1811 
1812  std::map<dof_id_type, dof_id_type> node_to_node_map, other_to_this_node_map; // The second is the inverse map of the first
1813  std::map<dof_id_type, std::vector<dof_id_type>> node_to_elems_map;
1814 
1815  typedef dof_id_type key_type;
1816  typedef std::pair<const Elem *, unsigned char> val_type;
1817  typedef std::pair<key_type, val_type> key_val_pair;
1818  typedef std::unordered_multimap<key_type, val_type> map_type;
1819  // Mapping between all side keys in this mesh and elements+side numbers relevant to the boundary in this mesh as well.
1820  map_type side_to_elem_map;
1821 
1822  // If there is only one mesh (i.e. other_mesh == nullptr), then loop over this mesh twice
1823  if (!other_mesh)
1824  {
1825  other_mesh = this;
1826  }
1827 
1828  if ((this_mesh_boundary_id != BoundaryInfo::invalid_id) &&
1829  (other_mesh_boundary_id != BoundaryInfo::invalid_id))
1830  {
1831  LOG_SCOPE("stitch_meshes node merging", "UnstructuredMesh");
1832 
1833  // While finding nodes on the boundary, also find the minimum edge length
1834  // of all faces on both boundaries. This will later be used in relative
1835  // distance checks when stitching nodes.
1836  Real h_min = std::numeric_limits<Real>::max();
1837  bool h_min_updated = false;
1838 
1839  // Loop below fills in these sets for the two meshes.
1840  std::set<dof_id_type> this_boundary_node_ids, other_boundary_node_ids;
1841 
1842  // Pull objects out of the loop to reduce heap operations
1843  std::unique_ptr<const Elem> side;
1844 
1845  {
1846  // Make temporary fixed-size arrays for loop
1847  boundary_id_type id_array[2] = {this_mesh_boundary_id, other_mesh_boundary_id};
1848  std::set<dof_id_type> * set_array[2] = {&this_boundary_node_ids, &other_boundary_node_ids};
1849  const MeshBase * mesh_array[2] = {this, other_mesh};
1850 
1851  for (unsigned i=0; i<2; ++i)
1852  {
1853  // First we deal with node boundary IDs. We only enter
1854  // this loop if we have at least one nodeset. Note that we
1855  // do not attempt to make an h_min determination here.
1856  // The h_min determination is done while looping over the
1857  // Elems and checking their sides and edges for boundary
1858  // information, below.
1859  if (mesh_array[i]->get_boundary_info().n_nodeset_conds() > 0)
1860  {
1861  // build_node_list() returns a vector of (node-id, bc-id) tuples
1862  for (const auto & t : mesh_array[i]->get_boundary_info().build_node_list())
1863  {
1864  boundary_id_type node_bc_id = std::get<1>(t);
1865  if (node_bc_id == id_array[i])
1866  {
1867  dof_id_type this_node_id = std::get<0>(t);
1868  set_array[i]->insert( this_node_id );
1869  }
1870  }
1871  }
1872 
1873  // Container to catch boundary IDs passed back from BoundaryInfo.
1874  std::vector<boundary_id_type> bc_ids;
1875 
1876  // Pointers to boundary NodeElems encountered while looping over the entire Mesh
1877  // and checking side and edge boundary ids. The Nodes associated with NodeElems
1878  // may be in a boundary nodeset, but not connected to any other Elems. In this
1879  // case, we also consider the "minimum node separation distance" amongst all
1880  // NodeElems when determining the relevant h_min value for this mesh.
1881  std::vector<const Elem *> boundary_node_elems;
1882 
1883  for (auto & el : mesh_array[i]->element_ptr_range())
1884  {
1885  // Now check whether elem has a face on the specified boundary
1886  for (auto side_id : el->side_index_range())
1887  if (el->neighbor_ptr(side_id) == nullptr)
1888  {
1889  // Get *all* boundary IDs on this side, not just the first one!
1890  mesh_array[i]->get_boundary_info().boundary_ids (el, side_id, bc_ids);
1891 
1892  if (std::find(bc_ids.begin(), bc_ids.end(), id_array[i]) != bc_ids.end())
1893  {
1894  el->build_side_ptr(side, side_id);
1895  for (auto & n : side->node_ref_range())
1896  set_array[i]->insert(n.id());
1897 
1898  h_min = std::min(h_min, side->hmin());
1899  h_min_updated = true;
1900 
1901  // This side is on the boundary, add its information to side_to_elem
1902  if (skip_find_neighbors && (i==0))
1903  {
1904  key_type key = el->low_order_key(side_id);
1905  val_type val;
1906  val.first = el;
1907  val.second = cast_int<unsigned char>(side_id);
1908 
1909  key_val_pair kvp;
1910  kvp.first = key;
1911  kvp.second = val;
1912  side_to_elem_map.insert (kvp);
1913  }
1914  }
1915 
1916  // Also, check the edges on this side. We don't have to worry about
1917  // updating neighbor info in this case since elements don't store
1918  // neighbor info on edges.
1919  for (auto edge_id : el->edge_index_range())
1920  {
1921  if (el->is_edge_on_side(edge_id, side_id))
1922  {
1923  // Get *all* boundary IDs on this edge, not just the first one!
1924  mesh_array[i]->get_boundary_info().edge_boundary_ids (el, edge_id, bc_ids);
1925 
1926  if (std::find(bc_ids.begin(), bc_ids.end(), id_array[i]) != bc_ids.end())
1927  {
1928  std::unique_ptr<const Elem> edge (el->build_edge_ptr(edge_id));
1929  for (auto & n : edge->node_ref_range())
1930  set_array[i]->insert( n.id() );
1931 
1932  h_min = std::min(h_min, edge->hmin());
1933  h_min_updated = true;
1934  }
1935  }
1936  } // end for (edge_id)
1937  } // end if (side == nullptr)
1938 
1939  // Alternatively, is this a boundary NodeElem? If so,
1940  // add it to a list of NodeElems that will later be
1941  // used to set h_min based on the minimum node
1942  // separation distance between all pairs of boundary
1943  // NodeElems.
1944  if (el->type() == NODEELEM)
1945  {
1946  mesh_array[i]->get_boundary_info().boundary_ids(el->node_ptr(0), bc_ids);
1947  if (std::find(bc_ids.begin(), bc_ids.end(), id_array[i]) != bc_ids.end())
1948  {
1949  boundary_node_elems.push_back(el);
1950 
1951  // Debugging:
1952  // libMesh::out << "Elem " << el->id() << " is a NodeElem on boundary " << id_array[i] << std::endl;
1953  }
1954  }
1955  } // end for (el)
1956 
1957  // Compute the minimum node separation distance amongst
1958  // all boundary NodeElem pairs.
1959  {
1960  const auto N = boundary_node_elems.size();
1961  for (auto node_elem_i : make_range(N))
1962  for (auto node_elem_j : make_range(node_elem_i+1, N))
1963  {
1964  Real node_sep =
1965  (boundary_node_elems[node_elem_i]->point(0) - boundary_node_elems[node_elem_j]->point(0)).norm();
1966 
1967  // We only want to consider non-coincident
1968  // boundary NodeElem pairs when determining the
1969  // minimum node separation distance.
1970  if (node_sep > 0.)
1971  {
1972  h_min = std::min(h_min, node_sep);
1973  h_min_updated = true;
1974  }
1975  } // end for (node_elem_j)
1976  } // end minimum NodeElem separation scope
1977  } // end for (i)
1978  } // end scope
1979 
1980  if (verbose)
1981  {
1982  libMesh::out << "In UnstructuredMesh::stitch_meshes:\n"
1983  << "This mesh has " << this_boundary_node_ids.size()
1984  << " nodes on boundary `"
1985  << this->get_boundary_info().get_sideset_name(this_mesh_boundary_id)
1986  << "' (" << this_mesh_boundary_id << ").\n"
1987  << "Other mesh has " << other_boundary_node_ids.size()
1988  << " nodes on boundary `"
1989  << this->get_boundary_info().get_sideset_name(other_mesh_boundary_id)
1990  << "' (" << other_mesh_boundary_id << ").\n";
1991 
1992  if (h_min_updated)
1993  {
1994  libMesh::out << "Minimum edge length on both surfaces is " << h_min << ".\n";
1995  }
1996  else
1997  {
1998  libMesh::out << "No minimum edge length determined on specified surfaces." << std::endl;
1999  }
2000  }
2001 
2002  // At this point, if h_min==0 it means that there were at least two coincident
2003  // nodes on the surfaces being stitched, and we don't currently support that case.
2004  // (It might be possible to support, but getting it exactly right would be tricky
2005  // and probably not worth the extra complications to the "normal" case.)
2006  libmesh_error_msg_if(h_min < std::numeric_limits<Real>::epsilon(),
2007  "Coincident nodes detected on source and/or target "
2008  "surface, stitching meshes is not possible.");
2009 
2010  // We require nanoflann for the "binary search" (really kd-tree)
2011  // option to work. If it's not available, turn that option off,
2012  // warn the user, and fall back on the N^2 search algorithm.
2013  if (use_binary_search)
2014  {
2015 #ifndef LIBMESH_HAVE_NANOFLANN
2016  use_binary_search = false;
2017  libmesh_warning("The use_binary_search option in the "
2018  "UnstructuredMesh stitching algorithms requires nanoflann "
2019  "support. Falling back on N^2 search algorithm.");
2020 #endif
2021  }
2022 
2023  if (!this_boundary_node_ids.empty())
2024  {
2025  if (use_binary_search)
2026  {
2027 #ifdef LIBMESH_HAVE_NANOFLANN
2028  typedef nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<Real, VectorOfNodesAdaptor>,
2029  VectorOfNodesAdaptor, 3, std::size_t> kd_tree_t;
2030 
2031  // Create the dataset needed to build the kd tree with nanoflann
2032  std::vector<std::pair<Point, dof_id_type>> this_mesh_nodes(this_boundary_node_ids.size());
2033 
2034  for (auto [it, ctr] = std::make_tuple(this_boundary_node_ids.begin(), 0u);
2035  it != this_boundary_node_ids.end(); ++it, ++ctr)
2036  {
2037  this_mesh_nodes[ctr].first = this->point(*it);
2038  this_mesh_nodes[ctr].second = *it;
2039  }
2040 
2041  VectorOfNodesAdaptor vec_nodes_adaptor(this_mesh_nodes);
2042 
2043  kd_tree_t this_kd_tree(3, vec_nodes_adaptor, 10);
2044  this_kd_tree.buildIndex();
2045 
2046  // Storage for nearest neighbor in the loop below
2047  std::size_t ret_index;
2048  Real ret_dist_sqr;
2049 
2050  // Loop over other mesh. For each node, find its nearest neighbor in this mesh, and fill in the maps.
2051  for (const auto & node_id : other_boundary_node_ids)
2052  {
2053  const auto & p = other_mesh->point(node_id);
2054  const Real query_pt[] = {p(0), p(1), p(2)};
2055  this_kd_tree.knnSearch(&query_pt[0], 1, &ret_index, &ret_dist_sqr);
2056 
2057  // TODO: here we should use the user's specified tolerance
2058  // and the previously determined value of h_min in the
2059  // distance comparison, not just TOLERANCE^2.
2060  if (ret_dist_sqr < TOLERANCE*TOLERANCE)
2061  {
2062  node_to_node_map[this_mesh_nodes[ret_index].second] = node_id;
2063  other_to_this_node_map[node_id] = this_mesh_nodes[ret_index].second;
2064  }
2065  }
2066 
2067  // If the two maps don't have the same size, it means one
2068  // node in this mesh is the nearest neighbor of several
2069  // nodes in other mesh. Since the stitching is ambiguous in
2070  // this case, we throw an error.
2071  libmesh_error_msg_if(node_to_node_map.size() != other_to_this_node_map.size(),
2072  "Error: Found multiple matching nodes in stitch_meshes");
2073 #endif
2074  }
2075  else // !use_binary_search
2076  {
2077  // In the unlikely event that two meshes composed entirely of
2078  // NodeElems are being stitched together, we will not have
2079  // selected a valid h_min value yet, and the distance
2080  // comparison below will be true for essentially any two
2081  // nodes. In this case we simply fall back on an absolute
2082  // distance check.
2083  if (!h_min_updated)
2084  {
2085  libmesh_warning("No valid h_min value was found, falling back on "
2086  "absolute distance check in the N^2 search algorithm.");
2087  h_min = 1.;
2088  }
2089 
2090  // Otherwise, use a simple N^2 search to find the closest matching points. This can be helpful
2091  // in the case that we have tolerance issues which cause mismatch between the two surfaces
2092  // that are being stitched.
2093  for (const auto & this_node_id : this_boundary_node_ids)
2094  {
2095  Node & this_node = this->node_ref(this_node_id);
2096 
2097  bool found_matching_nodes = false;
2098 
2099  for (const auto & other_node_id : other_boundary_node_ids)
2100  {
2101  const Node & other_node = other_mesh->node_ref(other_node_id);
2102 
2103  Real node_distance = (this_node - other_node).norm();
2104 
2105  if (node_distance < tol*h_min)
2106  {
2107  // Make sure we didn't already find a matching node!
2108  libmesh_error_msg_if(found_matching_nodes,
2109  "Error: Found multiple matching nodes in stitch_meshes");
2110 
2111  node_to_node_map[this_node_id] = other_node_id;
2112  other_to_this_node_map[other_node_id] = this_node_id;
2113 
2114  found_matching_nodes = true;
2115  }
2116  }
2117  }
2118  }
2119  }
2120 
2121  // Build up the node_to_elems_map, using only one loop over other_mesh
2122  for (auto & el : other_mesh->element_ptr_range())
2123  {
2124  // For each node on the element, find the corresponding node
2125  // on "this" Mesh, 'this_node_id', if it exists, and push
2126  // the current element ID back onto node_to_elems_map[this_node_id].
2127  // For that we will use the reverse mapping we created at
2128  // the same time as the forward mapping.
2129  for (auto & n : el->node_ref_range())
2130  if (const auto it = other_to_this_node_map.find(/*other_node_id=*/n.id());
2131  it != other_to_this_node_map.end())
2132  node_to_elems_map[/*this_node_id=*/it->second].push_back( el->id() );
2133  }
2134 
2135  if (verbose)
2136  {
2137  libMesh::out << "In UnstructuredMesh::stitch_meshes:\n"
2138  << "Found " << node_to_node_map.size()
2139  << " matching nodes.\n"
2140  << std::endl;
2141  }
2142 
2143  if (enforce_all_nodes_match_on_boundaries)
2144  {
2145  std::size_t n_matching_nodes = node_to_node_map.size();
2146  std::size_t this_mesh_n_nodes = this_boundary_node_ids.size();
2147  std::size_t other_mesh_n_nodes = other_boundary_node_ids.size();
2148  libmesh_error_msg_if((n_matching_nodes != this_mesh_n_nodes) || (n_matching_nodes != other_mesh_n_nodes),
2149  "Error: We expected the number of nodes to match.");
2150  }
2151 
2152  if (merge_boundary_nodes_all_or_nothing)
2153  {
2154  std::size_t n_matching_nodes = node_to_node_map.size();
2155  std::size_t this_mesh_n_nodes = this_boundary_node_ids.size();
2156  std::size_t other_mesh_n_nodes = other_boundary_node_ids.size();
2157  if ((n_matching_nodes != this_mesh_n_nodes) || (n_matching_nodes != other_mesh_n_nodes))
2158  {
2159  if (verbose)
2160  {
2161  libMesh::out << "Skipping node merging in "
2162  "UnstructuredMesh::stitch_meshes because not "
2163  "all boundary nodes were matched."
2164  << std::endl;
2165  }
2166  node_to_node_map.clear();
2167  other_to_this_node_map.clear();
2168  node_to_elems_map.clear();
2169  }
2170  }
2171  }
2172  else
2173  {
2174  if (verbose)
2175  {
2176  libMesh::out << "Skip node merging in UnstructuredMesh::stitch_meshes:" << std::endl;
2177  }
2178  }
2179 
2180  dof_id_type node_delta = this->max_node_id();
2181  dof_id_type elem_delta = this->max_elem_id();
2182 
2183  unique_id_type unique_delta =
2184 #ifdef LIBMESH_ENABLE_UNIQUE_ID
2185  this->parallel_max_unique_id();
2186 #else
2187  0;
2188 #endif
2189 
2190  // If other_mesh != nullptr, then we have to do a bunch of work
2191  // in order to copy it to this mesh
2192  if (this!=other_mesh)
2193  {
2194  LOG_SCOPE("stitch_meshes copying", "UnstructuredMesh");
2195 
2196  // Increment the node_to_node_map and node_to_elems_map
2197  // to account for id offsets
2198  for (auto & pr : node_to_node_map)
2199  pr.second += node_delta;
2200 
2201  for (auto & pr : node_to_elems_map)
2202  for (auto & entry : pr.second)
2203  entry += elem_delta;
2204 
2205  // We run into problems when the libMesh subdomain standard (the
2206  // id defines the subdomain; the name was an afterthought) and
2207  // the MOOSE standard (the name defines the subdomain; the id
2208  // might be autogenerated) clash.
2209  //
2210  // Subdomain ids with the same name in both meshes are surely
2211  // meant to represent the same subdomain. We can just merge
2212  // them.
2213  //
2214  // Subdomain ids which don't have a name in either mesh are
2215  // almost surely meant to represent the same subdomain. We'll
2216  // just merge them.
2217  //
2218  // Subdomain ids with different names in different meshes, or
2219  // names with different ids in different meshes, are trickier.
2220  // For backwards compatibility we default to the old "just copy
2221  // all the subdomain ids over" behavior, but if requested we'll
2222  // remap any ids that appear to be clear conflicts, and we'll
2223  // scream and die if we see any ids that are ambiguous due to
2224  // being named in one mesh but not the other.
2225  std::unordered_map<subdomain_id_type, subdomain_id_type> id_remapping;
2226  if (remap_subdomain_ids)
2227  {
2228  const auto & this_map = this->get_subdomain_name_map();
2229  const auto & other_map = other_mesh->get_subdomain_name_map();
2230  std::unordered_map<std::string, subdomain_id_type> other_map_reversed;
2231  for (auto & [sid, sname] : other_map)
2232  other_map_reversed.emplace(sname, sid);
2233 
2234  std::unordered_map<std::string, subdomain_id_type> this_map_reversed;
2235  for (auto & [sid, sname] : this_map)
2236  this_map_reversed.emplace(sname, sid);
2237 
2238  // We don't require either mesh to be prepared, but that
2239  // means we need to check for subdomains manually.
2240  auto get_subdomains = [](const MeshBase & mesh) {
2241  std::set<subdomain_id_type> all_subdomains;
2242  for (auto & el : mesh.element_ptr_range())
2243  all_subdomains.insert(el->subdomain_id());
2244  return all_subdomains;
2245  };
2246 
2247  const auto this_subdomains = get_subdomains(*this);
2248  const auto other_subdomains = get_subdomains(*other_mesh);
2249 
2250  for (auto & [sid, sname] : this_map)
2251  {
2252  // The same name with the same id means we're fine. The
2253  // same name with another id means we remap their id to
2254  // ours
2255  if (const auto other_reverse_it = other_map_reversed.find(sname);
2256  other_reverse_it != other_map_reversed.end() && other_reverse_it->second != sid)
2257  id_remapping[other_reverse_it->second] = sid;
2258 
2259  // The same id with a different name, we'll get to
2260  // later. The same id without any name means we don't
2261  // know what the user wants.
2262  if (other_subdomains.count(sid) && !other_map.count(sid))
2263  libmesh_error_msg("Can't safely stitch with a mesh sharing subdomain id "
2264  << sid << " but not subdomain name " << sname);
2265  }
2266 
2267  subdomain_id_type next_free_id = 0;
2268  // We might try to stitch empty meshes ...
2269  if (!this_subdomains.empty())
2270  next_free_id = *this_subdomains.rbegin() + 1;
2271  if (!other_subdomains.empty())
2272  next_free_id =
2273  std::max(next_free_id,
2274  cast_int<subdomain_id_type>
2275  (*other_subdomains.rbegin() + 1));
2276 
2277  for (auto & [sid, sname] : other_map)
2278  {
2279  // At this point we've figured out any remapping
2280  // necessary for an sname that we share. And we don't
2281  // need to remap any sid we don't share.
2282  if (!this_map_reversed.count(sname))
2283  {
2284  // But if we don't have this sname and we do have this
2285  // sid then we can't just merge into that.
2286  if (this_subdomains.count(sid))
2287  {
2288  // If we have this sid with no name, we don't
2289  // know what the user wants.
2290  if (!this_map.count(sid))
2291  libmesh_error_msg("Can't safely stitch with a mesh sharing subdomain id "
2292  << sid << " but under subdomain name " << sname);
2293 
2294  // We have this sid under a different name, so
2295  // we just need to give the other elements a new
2296  // id.
2297 
2298  // Users might have done crazy things with id
2299  // choice so let's make sure they didn't get too
2300  // crazy.
2301  libmesh_error_msg_if ((!this_subdomains.empty() &&
2302  next_free_id < *this_subdomains.rbegin()) ||
2303  (!other_subdomains.empty() &&
2304  next_free_id < *other_subdomains.rbegin()),
2305  "Subdomain id overflow");
2306 
2307  id_remapping[sid] = next_free_id++;
2308  this->subdomain_name(next_free_id) = sname;
2309  }
2310  // If we don't have this subdomain id, well, we're
2311  // about to, so we should have its name too.
2312  else
2313  this->subdomain_name(sid) = sname;
2314  }
2315  }
2316  }
2317 
2318  // Copy mesh data. If we skip the call to find_neighbors(), the lists
2319  // of neighbors will be copied verbatim from the other mesh
2320  this->copy_nodes_and_elements(*other_mesh, skip_find_neighbors,
2321  elem_delta, node_delta,
2322  unique_delta, &id_remapping);
2323 
2324  // Copy BoundaryInfo from other_mesh too. We do this via the
2325  // list APIs rather than element-by-element for speed.
2326  BoundaryInfo & boundary = this->get_boundary_info();
2327  const BoundaryInfo & other_boundary = other_mesh->get_boundary_info();
2328 
2329  for (const auto & t : other_boundary.build_node_list())
2330  boundary.add_node(std::get<0>(t) + node_delta,
2331  std::get<1>(t));
2332 
2333  for (const auto & t : other_boundary.build_side_list())
2334  boundary.add_side(std::get<0>(t) + elem_delta,
2335  std::get<1>(t),
2336  std::get<2>(t));
2337 
2338  for (const auto & t : other_boundary.build_edge_list())
2339  boundary.add_edge(std::get<0>(t) + elem_delta,
2340  std::get<1>(t),
2341  std::get<2>(t));
2342 
2343  for (const auto & t : other_boundary.build_shellface_list())
2344  boundary.add_shellface(std::get<0>(t) + elem_delta,
2345  std::get<1>(t),
2346  std::get<2>(t));
2347 
2348  const auto & other_ns_id_to_name = other_boundary.get_nodeset_name_map();
2349  auto & ns_id_to_name = boundary.set_nodeset_name_map();
2350  ns_id_to_name.insert(other_ns_id_to_name.begin(), other_ns_id_to_name.end());
2351 
2352  const auto & other_ss_id_to_name = other_boundary.get_sideset_name_map();
2353  auto & ss_id_to_name = boundary.set_sideset_name_map();
2354  ss_id_to_name.insert(other_ss_id_to_name.begin(), other_ss_id_to_name.end());
2355 
2356  const auto & other_es_id_to_name = other_boundary.get_edgeset_name_map();
2357  auto & es_id_to_name = boundary.set_edgeset_name_map();
2358  es_id_to_name.insert(other_es_id_to_name.begin(), other_es_id_to_name.end());
2359 
2360  // Merge other_mesh's elemset information with ours. Throw an
2361  // error if this and other_mesh have overlapping elemset codes
2362  // that refer to different elemset ids.
2363  std::vector<dof_id_type> this_elemset_codes = this->get_elemset_codes();
2364  MeshBase::elemset_type this_id_set_to_fill, other_id_set_to_fill;
2365  for (const auto & elemset_code : other_mesh->get_elemset_codes())
2366  {
2367  // Get the elemset ids for this elemset_code on other_mesh
2368  other_mesh->get_elemsets(elemset_code, other_id_set_to_fill);
2369 
2370  // Check that this elemset code does not already exist
2371  // in this mesh, or if it does, that it has the same elemset
2372  // ids associated with it.
2373  //
2374  // Note: get_elemset_codes() is guaranteed to return a
2375  // sorted vector, so we can binary search in it.
2376  auto it = Utility::binary_find(this_elemset_codes.begin(),
2377  this_elemset_codes.end(),
2378  elemset_code);
2379 
2380  if (it != this_elemset_codes.end())
2381  {
2382  // This mesh has the same elemset code. Does it refer to
2383  // the same elemset ids?
2384  this->get_elemsets(elemset_code, this_id_set_to_fill);
2385 
2386  // Throw an error if they don't match, otherwise we
2387  // don't need to do anything
2388  libmesh_error_msg_if(other_id_set_to_fill != this_id_set_to_fill,
2389  "Attempted to stitch together meshes with conflicting elemset codes.");
2390  }
2391  else
2392  {
2393  // Add other_mesh's elemset code to this mesh
2394  this->add_elemset_code(elemset_code, other_id_set_to_fill);
2395  }
2396  }
2397  } // end if (other_mesh)
2398 
2399  // Finally, we need to "merge" the overlapping nodes
2400  // We do this by iterating over node_to_elems_map and updating
2401  // the elements so that they "point" to the nodes that came
2402  // from this mesh, rather than from other_mesh.
2403  // Then we iterate over node_to_node_map and delete the
2404  // duplicate nodes that came from other_mesh.
2405 
2406  {
2407  LOG_SCOPE("stitch_meshes node updates", "UnstructuredMesh");
2408 
2409  // Container to catch boundary IDs passed back from BoundaryInfo.
2410  std::vector<boundary_id_type> bc_ids;
2411 
2412  for (const auto & [target_node_id, elem_vec] : node_to_elems_map)
2413  {
2414  dof_id_type other_node_id = node_to_node_map[target_node_id];
2415  Node & target_node = this->node_ref(target_node_id);
2416 
2417  std::size_t n_elems = elem_vec.size();
2418  for (std::size_t i=0; i<n_elems; i++)
2419  {
2420  dof_id_type elem_id = elem_vec[i];
2421  Elem * el = this->elem_ptr(elem_id);
2422 
2423  // find the local node index that we want to update
2424  unsigned int local_node_index = el->local_node(other_node_id);
2425  libmesh_assert_not_equal_to(local_node_index, libMesh::invalid_uint);
2426 
2427  // We also need to copy over the nodeset info here,
2428  // because the node will get deleted below
2429  this->get_boundary_info().boundary_ids(el->node_ptr(local_node_index), bc_ids);
2430  el->set_node(local_node_index, &target_node);
2431  this->get_boundary_info().add_node(&target_node, bc_ids);
2432  }
2433  }
2434  }
2435 
2436  {
2437  LOG_SCOPE("stitch_meshes node deletion", "UnstructuredMesh");
2438  for (const auto & [other_node_id, this_node_id] : node_to_node_map)
2439  {
2440  // In the case that this==other_mesh, the two nodes might be the same (e.g. if
2441  // we're stitching a "sliver"), hence we need to skip node deletion in that case.
2442  if ((this == other_mesh) && (this_node_id == other_node_id))
2443  continue;
2444 
2445  this->delete_node( this->node_ptr(this_node_id) );
2446  }
2447  }
2448 
2449  // If find_neighbors() wasn't called in prepare_for_use(), we need to
2450  // manually loop once more over all elements adjacent to the stitched boundary
2451  // and fix their lists of neighbors.
2452  // This is done according to the following steps:
2453  // 1. Loop over all copied elements adjacent to the boundary using node_to_elems_map (trying to avoid duplicates)
2454  // 2. Look at all their sides with a nullptr neighbor and update them using side_to_elem_map if necessary
2455  // 3. Update the corresponding side in side_to_elem_map as well
2456  if (skip_find_neighbors)
2457  {
2458  LOG_SCOPE("stitch_meshes neighbor fixes", "UnstructuredMesh");
2459 
2460  // Pull objects out of the loop to reduce heap operations
2461  std::unique_ptr<const Elem> my_side, their_side;
2462 
2463  std::set<dof_id_type> fixed_elems;
2464  for (const auto & pr : node_to_elems_map)
2465  {
2466  std::size_t n_elems = pr.second.size();
2467  for (std::size_t i=0; i<n_elems; i++)
2468  {
2469  dof_id_type elem_id = pr.second[i];
2470  if (!fixed_elems.count(elem_id))
2471  {
2472  Elem * el = this->elem_ptr(elem_id);
2473  fixed_elems.insert(elem_id);
2474  for (auto s : el->side_index_range())
2475  {
2476  if (el->neighbor_ptr(s) == nullptr)
2477  {
2478  key_type key = el->low_order_key(s);
2479  auto bounds = side_to_elem_map.equal_range(key);
2480 
2481  if (bounds.first != bounds.second)
2482  {
2483  // Get the side for this element
2484  el->side_ptr(my_side, s);
2485 
2486  // Look at all the entries with an equivalent key
2487  while (bounds.first != bounds.second)
2488  {
2489  // Get the potential element
2490  Elem * neighbor = const_cast<Elem *>(bounds.first->second.first);
2491 
2492  // Get the side for the neighboring element
2493  const unsigned int ns = bounds.first->second.second;
2494  neighbor->side_ptr(their_side, ns);
2495  //libmesh_assert(my_side.get());
2496  //libmesh_assert(their_side.get());
2497 
2498  // If found a match with my side
2499  //
2500  // We need special tests here for 1D:
2501  // since parents and children have an equal
2502  // side (i.e. a node), we need to check
2503  // ns != ms, and we also check level() to
2504  // avoid setting our neighbor pointer to
2505  // any of our neighbor's descendants
2506  if ((*my_side == *their_side) &&
2507  (el->level() == neighbor->level()) &&
2508  ((el->dim() != 1) || (ns != s)))
2509  {
2510  // So share a side. Is this a mixed pair
2511  // of subactive and active/ancestor
2512  // elements?
2513  // If not, then we're neighbors.
2514  // If so, then the subactive's neighbor is
2515 
2516  if (el->subactive() ==
2517  neighbor->subactive())
2518  {
2519  // an element is only subactive if it has
2520  // been coarsened but not deleted
2521  el->set_neighbor (s,neighbor);
2522  neighbor->set_neighbor(ns,el);
2523  }
2524  else if (el->subactive())
2525  {
2526  el->set_neighbor(s,neighbor);
2527  }
2528  else if (neighbor->subactive())
2529  {
2530  neighbor->set_neighbor(ns,el);
2531  }
2532  // It's OK to invalidate the
2533  // bounds.first iterator here,
2534  // as we are immediately going
2535  // to break out of this while
2536  // loop. bounds.first will
2537  // therefore not be used for
2538  // anything else.
2539  side_to_elem_map.erase (bounds.first);
2540  break;
2541  }
2542 
2543  ++bounds.first;
2544  }
2545  }
2546  }
2547  }
2548  }
2549  }
2550  }
2551  }
2552 
2553  const bool old_allow_find_neighbors = this->allow_find_neighbors();
2554  const bool old_allow_remote_element_removal = this->allow_remote_element_removal();
2555  this->allow_find_neighbors(!skip_find_neighbors);
2556  this->allow_remote_element_removal(false);
2557  this->prepare_for_use();
2558  this->allow_find_neighbors(old_allow_find_neighbors);
2559  this->allow_remote_element_removal(old_allow_remote_element_removal);
2560 
2561  // After the stitching, we may want to clear boundary IDs from element
2562  // faces that are now internal to the mesh
2563  if (clear_stitched_boundary_ids)
2564  {
2565  LOG_SCOPE("stitch_meshes clear bcids", "UnstructuredMesh");
2566 
2568  this_mesh_boundary_id, other_mesh_boundary_id, /*clear_nodeset_data=*/true);
2569  }
2570 
2571  // Return the number of nodes which were merged.
2572  return node_to_node_map.size();
2573 }
2574 
2575 
2576 } // namespace libMesh
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
OStreamProxy err
IndexType begin() const
Definition: simple_range.h:41
void set_p_level(const unsigned int p)
Sets the value of the p-refinement level for the element.
bool closed()
Checks that the library has been closed.
Definition: libmesh.C:283
ElemType
Defines an enum for geometric element types.
The SimpleRange templated class is intended to make it easy to construct ranges from pairs of iterato...
Definition: simple_range.h:36
This class supports simple reads and writes in any libMesh-supported format, by dispatching to one of...
Definition: namebased_io.h:44
unique_id_type & set_unique_id()
Definition: dof_object.h:858
virtual void reserve_nodes(const dof_id_type nn)=0
Reserves space for a known number of nodes.
std::size_t stitch_surfaces(boundary_id_type boundary_id_1, boundary_id_type boundary_id_2, Real tol=TOLERANCE, bool clear_stitched_boundary_ids=false, bool verbose=true, bool use_binary_search=true, bool enforce_all_nodes_match_on_boundaries=false, bool merge_boundary_nodes_all_or_nothing=false)
Similar to stitch_meshes, except that we stitch two adjacent surfaces within this mesh...
const Elem * parent() const
Definition: elem.h:3030
bool is_prepared() const
Definition: mesh_base.h:198
virtual void all_complete_order_range(const SimpleRange< element_iterator > &range) override
Converts a (conforming, non-refined) mesh with linear elements into a mesh with "complete" order elem...
virtual void copy_nodes_and_elements(const MeshBase &other_mesh, const bool skip_find_neighbors=false, dof_id_type element_id_offset=0, dof_id_type node_id_offset=0, unique_id_type unique_id_offset=0, std::unordered_map< subdomain_id_type, subdomain_id_type > *id_remapping=nullptr)
Deep copy of nodes and elements from another mesh object (used by subclass copy constructors and by m...
In parallel meshes where a ghost element has neighbors which do not exist on the local processor...
Definition: remote_elem.h:59
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2558
A Node is like a Point, but with more information.
Definition: node.h:52
bool is_ancestor_of(const Elem *descendant) const
Definition: elem.h:3009
static ElemType complete_order_equivalent_type(const ElemType et)
Definition: elem.C:3319
virtual unique_id_type parallel_max_unique_id() const =0
void libmesh_assert_valid_amr_interior_parents(const MeshBase &mesh)
A function for verifying that any interior_parent pointers on elements are consistent with AMR (paren...
Definition: mesh_tools.C:1361
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
void set_parent(Elem *p)
Sets the pointer to the element&#39;s parent.
Definition: elem.h:3046
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
const Elem * interior_parent() const
Definition: elem.C:1186
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
virtual void read(const std::string &mesh_file) override
This method implements reading a mesh from a specified file.
Definition: namebased_io.C:69
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2710
The definition of the const_element_iterator struct.
Definition: mesh_base.h:2216
static constexpr Real TOLERANCE
unsigned int dim
static void set_node_processor_ids(MeshBase &mesh)
This function is called after partitioning to set the processor IDs for the nodes.
Definition: partitioner.C:852
virtual bool is_child_on_side(const unsigned int c, const unsigned int s) const =0
size_t kdtree_get_point_count() const
Must return the number of data points.
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly ecreated (or read) mesh for use.
Definition: mesh_base.C:759
void set_isnt_prepared()
Tells this we have done some operation where we should no longer consider ourself prepared...
Definition: mesh_base.h:204
virtual dof_id_type low_order_key(const unsigned int s) const =0
const boundary_id_type side_id
VectorOfNodesAdaptor(const std::vector< std::pair< Point, dof_id_type >> &nodes)
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
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
void add_child(Elem *elem)
Adds a child pointer to the array of children of this element.
Definition: elem.C:2044
void set_refinement_flag(const RefinementState rflag)
Sets the value of the refinement flag for the element.
Definition: elem.h:3218
unique_id_type unique_id() const
Definition: dof_object.h:844
const Parallel::Communicator & comm() const
virtual unsigned int n_children() const =0
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
Fills a user-provided std::vector with the boundary ids associated with Node node.
std::map< const Elem *, const CouplingMatrix *, CompareDofObjectsByPIDAndThenID > map_type
What elements do we care about and what variables do we care about on each element?
virtual void own_node(Node &)
Takes ownership of node n on this partition of a distributed mesh, by setting n.processor_id() to thi...
Definition: mesh_base.h:740
virtual void all_first_order() override
Converts a mesh with higher-order elements into a mesh with linear elements.
virtual void set_next_unique_id(unique_id_type id)=0
Sets the next available unique id to be used.
virtual void find_neighbors(const bool reset_remote_elements=false, const bool reset_current_list=true) override
Other functions from MeshBase requiring re-definition.
The libMesh namespace provides an interface to certain functionality in the library.
void create_pid_mesh(UnstructuredMesh &pid_mesh, const processor_id_type pid) const
Generates a new mesh containing all the elements which are assigned to processor pid.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
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
virtual void all_second_order_range(const SimpleRange< element_iterator > &range, const bool full_ordered=true) override
Converts a (conforming, non-refined) mesh with linear elements into a mesh with second-order elements...
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.
void set_interior_parent(Elem *p)
Sets the pointer to the element&#39;s interior_parent.
Definition: elem.C:1248
uint8_t processor_id_type
Definition: id_types.h:104
This is the MeshBase class.
Definition: mesh_base.h:75
SimpleRange< ChildRefIter > child_ref_range()
Returns a range with all children of a parent element, usable in range-based for loops.
Definition: elem.h:2346
dof_id_type & set_id()
Definition: dof_object.h:836
void build_side_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &side_list, std::vector< boundary_id_type > &bc_id_list) const
Creates a list of element numbers, sides, and ids for those sides.
bool ancestor() const
Definition: elem.C:2010
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
virtual bool contract() override
Delete subactive (i.e.
std::size_t stitching_helper(const MeshBase *other_mesh, boundary_id_type boundary_id_1, boundary_id_type boundary_id_2, Real tol, bool clear_stitched_boundary_ids, bool verbose, bool use_binary_search, bool enforce_all_nodes_match_on_boundaries, bool skip_find_neighbors, bool merge_boundary_nodes_all_or_nothing, bool remap_subdomain_ids)
Helper function for stitch_meshes and stitch_surfaces that does the mesh stitching.
uint8_t processor_id_type
std::map< boundary_id_type, std::string > & set_sideset_name_map()
void replace_child(Elem *elem, unsigned int c)
Replaces the child pointer at the specified index in the child array.
Definition: elem.C:2090
processor_id_type n_processors() const
virtual bool is_serial() const
Definition: mesh_base.h:211
void libmesh_ignore(const Args &...)
void build_node_list(std::vector< dof_id_type > &node_id_list, std::vector< boundary_id_type > &bc_id_list) const
Creates a list of nodes and ids for those nodes.
void add_node(const Node *node, const boundary_id_type id)
Add Node node with boundary id id to the boundary information data structures.
const std::map< boundary_id_type, std::string > & get_nodeset_name_map() const
int8_t boundary_id_type
Definition: id_types.h:51
static const boundary_id_type invalid_id
Number used for internal use.
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
This is the MeshCommunication class.
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
virtual Real hmin() const
Definition: elem.C:702
std::size_t stitch_meshes(const MeshBase &other_mesh, boundary_id_type this_mesh_boundary, boundary_id_type other_mesh_boundary, Real tol=TOLERANCE, bool clear_stitched_boundary_ids=false, bool verbose=true, bool use_binary_search=true, bool enforce_all_nodes_match_on_boundaries=false, bool merge_boundary_nodes_all_or_nothing=false, bool remap_subdomain_ids=false)
Stitch other_mesh to this mesh so that this mesh is the union of the two meshes.
virtual unsigned int n_nodes() const =0
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
The UnstructuredMesh class is derived from the MeshBase class.
unsigned int which_neighbor_am_i(const Elem *e) const
This function tells you which neighbor e is.
Definition: elem.h:2919
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
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
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:802
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
bool allow_find_neighbors() const
Definition: mesh_base.h:1204
virtual void write(const std::string &mesh_file) override
This method implements writing a mesh to a specified file.
Definition: namebased_io.C:303
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
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
void libmesh_assert_valid_amr_elem_ids(const MeshBase &mesh)
A function for verifying that ids of elements are correctly sorted for AMR (parents have lower ids th...
Definition: mesh_tools.C:1341
std::string & subdomain_name(subdomain_id_type id)
Definition: mesh_base.C:1692
void allow_find_neighbors(bool allow)
If false is passed then this mesh will no longer work to find element neighbors when being prepared f...
Definition: mesh_base.h:1203
virtual std::unique_ptr< Elem > disconnected_clone() const
Definition: elem.C:412
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int n) const
Definition: elem.C:3040
auto norm(const T &a) -> decltype(std::abs(a))
Definition: tensor_tools.h:74
bool allow_remote_element_removal() const
Definition: mesh_base.h:1213
void set_neighbor(const unsigned int i, Elem *n)
Assigns n as the neighbor.
Definition: elem.h:2618
void build_shellface_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &shellface_list, std::vector< boundary_id_type > &bc_id_list) const
Creates a list of element numbers, shellfaces, and boundary ids for those shellfaces.
virtual void write_nodal_data(const std::string &, const std::vector< Number > &, const std::vector< std::string > &) override
This method implements writing a mesh with nodal data to a specified file where the nodal data and va...
Definition: namebased_io.C:464
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
static ElemType second_order_equivalent_type(const ElemType et, const bool full_ordered=true)
Definition: elem.C:3135
unsigned int which_child_am_i(const Elem *e) const
Definition: elem.h:3192
bool skip_partitioning() const
Definition: mesh_base.h:1256
unsigned int n_partitions() const
Definition: mesh_base.h:1345
virtual const Elem * elem_ptr(const dof_id_type i) const =0
UnstructuredMesh(const Parallel::Communicator &comm_in, unsigned char dim=1)
Constructor.
Real kdtree_get_pt(const size_t idx, int dim) const
virtual Elem * insert_elem(Elem *e)=0
Insert elem e to the element array, preserving its id and replacing/deleting any existing element wit...
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2598
std::set< GhostingFunctor * >::const_iterator ghosting_functors_begin() const
Beginning of range of ghosting functors.
Definition: mesh_base.h:1291
std::set< elemset_id_type > elemset_type
Typedef for the "set" container used to store elemset ids.
Definition: mesh_base.h:318
unsigned int level() const
Definition: elem.h:3074
virtual unsigned int n_vertices() const =0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const std::vector< std::pair< Point, dof_id_type > > _nodes
virtual std::unique_ptr< Elem > side_ptr(unsigned int i)=0
virtual ~UnstructuredMesh()
Destructor.
timpi_pure bool verify(const T &r) const
void copy_boundary_ids(const BoundaryInfo &old_boundary_info, const Elem *const old_elem, const Elem *const new_elem)
void max(const T &r, T &o, Request &req) const
virtual unsigned short dim() const =0
IndexType end() const
Definition: simple_range.h:43
Temporarily serialize a DistributedMesh for non-distributed-mesh capable code paths.
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2507
std::map< boundary_id_type, std::string > & set_nodeset_name_map()
void build_edge_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &edge_list, std::vector< boundary_id_type > &bc_id_list) const
Creates a list of element numbers, edges, and boundary ids for those edges.
void clear_stitched_boundary_side_ids(boundary_id_type sideset_id, boundary_id_type other_sideset_id, bool clear_nodeset_data=false)
Clear sideset information along a stitched mesh interface.
void make_nodes_parallel_consistent(MeshBase &)
Copy processor_ids and ids on ghost nodes from their local processors.
OStreamProxy out
virtual bool is_replicated() const
Definition: mesh_base.h:233
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
Add side side of element number elem with boundary id id to the boundary information data structure...
const std::string & get_sideset_name(boundary_id_type id) const
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const
Definition: elem.C:3048
void add_shellface(const dof_id_type elem, const unsigned short int shellface, const boundary_id_type id)
Add shell face shellface of element number elem with boundary id id to the boundary information data ...
bool subactive() const
Definition: elem.h:2959
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 add_extra_integers(const unsigned int n_integers)
Assigns a set of extra integers to this DofObject.
Definition: dof_object.C:490
ForwardIterator binary_find(ForwardIterator first, ForwardIterator last, const T &value)
The STL provides std::binary_search() which returns true or false depending on whether the searched-f...
Definition: utility.h:265
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:276
unsigned int n_extra_integers() const
Returns how many extra integers are associated to the DofObject.
Definition: dof_object.h:1170
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:596
virtual const Point & point(const dof_id_type i) const =0
unsigned int local_node(const dof_id_type i) const
Definition: elem.h:2487
void edge_boundary_ids(const Elem *const elem, const unsigned short int edge, std::vector< boundary_id_type > &vec_to_fill) const
SimpleRange< NeighborPtrIter > neighbor_ptr_range()
Returns a range with all neighbors of an element, usable in range-based for loops.
Definition: elem.h:3503
void set_p_refinement_flag(const RefinementState pflag)
Sets the value of the p-refinement flag for the element.
Definition: elem.h:3234
bool allow_renumbering() const
Definition: mesh_base.h:1197
void inherit_data_from(const Elem &src)
A helper function for copying generic element data (mapping, subdomain, processor) from an element to...
Definition: elem.h:3339
void libmesh_assert_valid_neighbors(const MeshBase &mesh, bool assert_valid_remote_elems=true)
A function for verifying that neighbor connectivity is correct (each element is a neighbor of or desc...
Definition: mesh_tools.C:2205
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
virtual void read(const std::string &name, void *mesh_data=nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false) override
Reads the file specified by name.
virtual dof_id_type max_node_id() const =0
virtual dof_id_type n_elem() const =0
virtual const Node * node_ptr(const dof_id_type i) const =0
processor_id_type processor_id() const
bool active() const
Definition: elem.h:2941
bool kdtree_get_bbox(BBOX &) const
static ElemType first_order_equivalent_type(const ElemType et)
Definition: elem.C:3066
processor_id_type processor_id() const
Definition: dof_object.h:905
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:2475
virtual void reserve_elem(const dof_id_type ne)=0
Reserves space for a known number of elements.
uint8_t unique_id_type
Definition: id_types.h:86
bool has_children() const
Definition: elem.h:2979
virtual void write(const std::string &name) const override
Write the file specified by name.
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
void set_extra_integer(const unsigned int index, const dof_id_type value)
Sets the value on this object of the extra integer associated with index, which should have been obta...
Definition: dof_object.h:1086
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...
virtual dof_id_type n_nodes() const =0
This class implements the original default geometry ghosting requirements in libMesh: point neighbors...
void create_submesh(UnstructuredMesh &new_mesh, const const_element_iterator &it, const const_element_iterator &it_end) const
Constructs a mesh called "new_mesh" from the current mesh by iterating over the elements between it a...
Point vertex_average() const
Definition: elem.C:688
dof_id_type get_extra_integer(const unsigned int index) const
Gets the value on this object of the extra integer associated with index, which should have been obta...
Definition: dof_object.h:1102
const Elem * child_ptr(unsigned int i) const
Definition: elem.h:3163
std::set< GhostingFunctor * >::const_iterator ghosting_functors_end() const
End of range of ghosting functors.
Definition: mesh_base.h:1297
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int i, const unsigned int j)
A useful inline function which replaces the macros used previously.
uint8_t dof_id_type
Definition: id_types.h:67
std::map< boundary_id_type, std::string > & set_edgeset_name_map()
void add_edge(const dof_id_type elem, const unsigned short int edge, const boundary_id_type id)
Add edge edge of element number elem with boundary id id to the boundary information data structure...
const RemoteElem * remote_elem
Definition: remote_elem.C:57