unstructured_mesh.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 // C++ includes
21 #include <fstream>
22 #include <sstream>
23 #include <iomanip>
24 
25 // C includes
26 #include <sys/types.h> // for pid_t
27 #include <unistd.h> // for getpid(), unlink()
28 
29 // Local includes
30 #include "libmesh/boundary_info.h"
34 #include "libmesh/elem.h"
35 #include "libmesh/mesh_tools.h" // For n_levels
36 #include "libmesh/parallel.h"
37 #include "libmesh/remote_elem.h"
38 
39 // For most I/O
40 #include "libmesh/namebased_io.h"
41 
42 #include LIBMESH_INCLUDE_UNORDERED_MAP
43 
44 
45 namespace libMesh
46 {
47 
48 
49 // ------------------------------------------------------------
50 // UnstructuredMesh class member functions
52  unsigned char d) :
53  MeshBase (comm_in,d)
54 {
56 }
57 
58 
59 
60 #ifndef LIBMESH_DISABLE_COMMWORLD
62  MeshBase (d)
63 {
64  libmesh_deprecated();
66 }
67 #endif
68 
69 
70 
72  const bool skip_find_neighbors)
73 {
74  LOG_SCOPE("copy_nodes_and_elements()", "UnstructuredMesh");
75 
76  // We're assuming our subclass data needs no copy
77  libmesh_assert_equal_to (_n_parts, other_mesh._n_parts);
78  libmesh_assert_equal_to (_is_prepared, other_mesh._is_prepared);
79 
80  // We're assuming the other mesh has proper element number ordering,
81  // so that we add parents before their children.
82 #ifdef DEBUG
84 #endif
85 
86  //Copy in Nodes
87  {
88  //Preallocate Memory if necessary
89  this->reserve_nodes(other_mesh.n_nodes());
90 
91  const_node_iterator it = other_mesh.nodes_begin();
92  const_node_iterator end = other_mesh.nodes_end();
93 
94  for (; it != end; ++it)
95  {
96  const Node * oldn = *it;
97 
98  // Add new nodes in old node Point locations
99 #ifdef LIBMESH_ENABLE_UNIQUE_ID
100  Node * newn =
101 #endif
102  this->add_point(*oldn, oldn->id(), oldn->processor_id());
103 
104 #ifdef LIBMESH_ENABLE_UNIQUE_ID
105  newn->set_unique_id() = oldn->unique_id();
106 #endif
107  }
108  }
109 
110  //Copy in Elements
111  {
112  //Preallocate Memory if necessary
113  this->reserve_elem(other_mesh.n_elem());
114 
115  // Declare a map linking old and new elements, needed to copy the neighbor lists
116  typedef LIBMESH_BEST_UNORDERED_MAP<const Elem *, Elem *> map_type;
117  map_type old_elems_to_new_elems;
118 
119  // Loop over the elements
122 
123  for (; it != end; ++it)
124  {
125  //Look at the old element
126  const Elem * old = *it;
127  //Build a new element
128  Elem * newparent = old->parent() ?
129  this->elem_ptr(old->parent()->id()) : libmesh_nullptr;
130  UniquePtr<Elem> ap = Elem::build(old->type(), newparent);
131  Elem * el = ap.release();
132 
133  el->subdomain_id() = old->subdomain_id();
134 
135  for (unsigned int s=0; s != old->n_sides(); ++s)
136  if (old->neighbor_ptr(s) == remote_elem)
137  el->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
138 
139 #ifdef LIBMESH_ENABLE_AMR
140  if (old->has_children())
141  for (unsigned int c=0; c != old->n_children(); ++c)
142  if (old->child_ptr(c) == remote_elem)
143  el->add_child(const_cast<RemoteElem *>(remote_elem), c);
144 
145  //Create the parent's child pointers if necessary
146  if (newparent)
147  {
148  unsigned int oldc = old->parent()->which_child_am_i(old);
149  newparent->add_child(el, oldc);
150  }
151 
152  // Copy the refinement flags
153  el->set_refinement_flag(old->refinement_flag());
154 
155  // Use hack_p_level since we may not have sibling elements
156  // added yet
157  el->hack_p_level(old->p_level());
158 
159  el->set_p_refinement_flag(old->p_refinement_flag());
160 #endif // #ifdef LIBMESH_ENABLE_AMR
161 
162  //Assign all the nodes
163  for (unsigned int i=0;i<el->n_nodes();i++)
164  el->set_node(i) = this->node_ptr(old->node_id(i));
165 
166  // And start it off in the same subdomain
167  el->processor_id() = old->processor_id();
168 
169  // Give it the same element and unique ids
170  el->set_id(old->id());
171 
172 #ifdef LIBMESH_ENABLE_UNIQUE_ID
173  el->set_unique_id() = old->unique_id();
174 #endif
175 
176  //Hold onto it
177  if (!skip_find_neighbors)
178  {
179  this->add_elem(el);
180  }
181  else
182  {
183  Elem * new_el = this->add_elem(el);
184  old_elems_to_new_elems[old] = new_el;
185  }
186 
187  // Add the link between the original element and this copy to the map
188  if (skip_find_neighbors)
189  old_elems_to_new_elems[old] = el;
190  }
191 
192  // Loop (again) over the elements to fill in the neighbors
193  if (skip_find_neighbors)
194  {
195  it = other_mesh.elements_begin();
196  for (; it != end; ++it)
197  {
198  Elem * old_elem = *it;
199  Elem * new_elem = old_elems_to_new_elems[old_elem];
200  for (unsigned int s=0; s != old_elem->n_neighbors(); ++s)
201  {
202  const Elem * old_neighbor = old_elem->neighbor_ptr(s);
203  Elem * new_neighbor = old_elems_to_new_elems[old_neighbor];
204  new_elem->set_neighbor(s, new_neighbor);
205  }
206  }
207  }
208  }
209 
210  //Finally prepare the new Mesh for use. Keep the same numbering and
211  //partitioning for now.
212  this->allow_renumbering(false);
213  this->allow_remote_element_removal(false);
214  this->skip_partitioning(true);
215 
216  this->prepare_for_use(false, skip_find_neighbors);
217 
218  //But in the long term, use the same renumbering and partitioning
219  //policies as our source mesh.
220  this->allow_renumbering(other_mesh.allow_renumbering());
222  this->skip_partitioning(other_mesh.skip_partitioning());
223 }
224 
225 
226 
228 {
229  // this->clear (); // Nothing to clear at this level
230 
231  libmesh_exceptionless_assert (!libMesh::closed());
232 }
233 
234 
235 
236 
237 
238 void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
239  const bool reset_current_list)
240 {
241  // We might actually want to run this on an empty mesh
242  // (e.g. the boundary mesh for a nonexistent bcid!)
243  // libmesh_assert_not_equal_to (this->n_nodes(), 0);
244  // libmesh_assert_not_equal_to (this->n_elem(), 0);
245 
246  // This function must be run on all processors at once
247  parallel_object_only();
248 
249  LOG_SCOPE("find_neighbors()", "Mesh");
250 
251  const element_iterator el_end = this->elements_end();
252 
253  //TODO:[BSK] This should be removed later?!
254  if (reset_current_list)
255  for (element_iterator el = this->elements_begin(); el != el_end; ++el)
256  {
257  Elem * e = *el;
258  for (unsigned int s=0; s<e->n_neighbors(); s++)
259  if (e->neighbor_ptr(s) != remote_elem ||
260  reset_remote_elements)
262  }
263 
264  // Find neighboring elements by first finding elements
265  // with identical side keys and then check to see if they
266  // are neighbors
267  {
268  // data structures -- Use the hash_multimap if available
269  typedef unsigned int key_type;
270  typedef std::pair<Elem *, unsigned char> val_type;
271  typedef std::pair<key_type, val_type> key_val_pair;
272 
273  typedef LIBMESH_BEST_UNORDERED_MULTIMAP<key_type, val_type> map_type;
274 
275  // A map from side keys to corresponding elements & side numbers
276  map_type side_to_elem_map;
277 
278 
279 
280  for (element_iterator el = this->elements_begin(); el != el_end; ++el)
281  {
282  Elem * element = *el;
283 
284  for (unsigned char ms=0; ms<element->n_neighbors(); ms++)
285  {
286  next_side:
287  // If we haven't yet found a neighbor on this side, try.
288  // Even if we think our neighbor is remote, that
289  // information may be out of date.
290  if (element->neighbor_ptr(ms) == libmesh_nullptr ||
291  element->neighbor_ptr(ms) == remote_elem)
292  {
293  // Get the key for the side of this element
294  const unsigned int key = element->key(ms);
295 
296  // Look for elements that have an identical side key
297  std::pair <map_type::iterator, map_type::iterator>
298  bounds = side_to_elem_map.equal_range(key);
299 
300  // May be multiple keys, check all the possible
301  // elements which _might_ be neighbors.
302  if (bounds.first != bounds.second)
303  {
304  // Get the side for this element
305  const UniquePtr<Elem> my_side(element->side_ptr(ms));
306 
307  // Look at all the entries with an equivalent key
308  while (bounds.first != bounds.second)
309  {
310  // Get the potential element
311  Elem * neighbor = bounds.first->second.first;
312 
313  // Get the side for the neighboring element
314  const unsigned int ns = bounds.first->second.second;
315  const UniquePtr<Elem> their_side(neighbor->side_ptr(ns));
316  //libmesh_assert(my_side.get());
317  //libmesh_assert(their_side.get());
318 
319  // If found a match with my side
320  //
321  // We need special tests here for 1D:
322  // since parents and children have an equal
323  // side (i.e. a node), we need to check
324  // ns != ms, and we also check level() to
325  // avoid setting our neighbor pointer to
326  // any of our neighbor's descendants
327  if ((*my_side == *their_side) &&
328  (element->level() == neighbor->level()) &&
329  ((element->dim() != 1) || (ns != ms)))
330  {
331  // So share a side. Is this a mixed pair
332  // of subactive and active/ancestor
333  // elements?
334  // If not, then we're neighbors.
335  // If so, then the subactive's neighbor is
336 
337  if (element->subactive() ==
338  neighbor->subactive())
339  {
340  // an element is only subactive if it has
341  // been coarsened but not deleted
342  element->set_neighbor (ms,neighbor);
343  neighbor->set_neighbor(ns,element);
344  }
345  else if (element->subactive())
346  {
347  element->set_neighbor(ms,neighbor);
348  }
349  else if (neighbor->subactive())
350  {
351  neighbor->set_neighbor(ns,element);
352  }
353  side_to_elem_map.erase (bounds.first);
354 
355  // get out of this nested crap
356  goto next_side;
357  }
358 
359  ++bounds.first;
360  }
361  }
362 
363  // didn't find a match...
364  // Build the map entry for this element
365  key_val_pair kvp;
366 
367  kvp.first = key;
368  kvp.second.first = element;
369  kvp.second.second = ms;
370 
371  // use the lower bound as a hint for
372  // where to put it.
373 #if defined(LIBMESH_HAVE_UNORDERED_MAP) || defined(LIBMESH_HAVE_TR1_UNORDERED_MAP) || defined(LIBMESH_HAVE_HASH_MAP) || defined(LIBMESH_HAVE_EXT_HASH_MAP)
374  side_to_elem_map.insert (kvp);
375 #else
376  side_to_elem_map.insert (bounds.first,kvp);
377 #endif
378  }
379  }
380  }
381  }
382 
383 #ifdef LIBMESH_ENABLE_AMR
384 
412  const unsigned int n_levels = MeshTools::n_levels(*this);
413  for (unsigned int level = 1; level < n_levels; ++level)
414  {
415  element_iterator end = this->level_elements_end(level);
416  for (element_iterator el = this->level_elements_begin(level);
417  el != end; ++el)
418  {
419  Elem * current_elem = *el;
420  libmesh_assert(current_elem);
421  Elem * parent = current_elem->parent();
422  libmesh_assert(parent);
423  const unsigned int my_child_num = parent->which_child_am_i(current_elem);
424 
425  for (unsigned int s=0; s < current_elem->n_neighbors(); s++)
426  {
427  if (current_elem->neighbor_ptr(s) == libmesh_nullptr ||
428  (current_elem->neighbor_ptr(s) == remote_elem &&
429  parent->is_child_on_side(my_child_num, s)))
430  {
431  Elem * neigh = parent->neighbor_ptr(s);
432 
433  // If neigh was refined and had non-subactive children
434  // made remote earlier, then our current elem should
435  // actually have one of those remote children as a
436  // neighbor
437  if (neigh &&
438  (neigh->ancestor() ||
439  // If neigh has subactive children which should have
440  // matched as neighbors of the current element but
441  // did not, then those likewise must be remote
442  // children.
443  (current_elem->subactive() && neigh->has_children() &&
444  (neigh->level()+1) == current_elem->level())))
445  {
446 #ifdef DEBUG
447  // Let's make sure that "had children made remote"
448  // situation is actually the case
449  libmesh_assert(neigh->has_children());
450  bool neigh_has_remote_children = false;
451  for (unsigned int c = 0; c != neigh->n_children(); ++c)
452  {
453  if (neigh->child_ptr(c) == remote_elem)
454  neigh_has_remote_children = true;
455  }
456  libmesh_assert(neigh_has_remote_children);
457 
458  // And let's double-check that we don't have
459  // a remote_elem neighboring an active local element
460  if (current_elem->active())
461  libmesh_assert_not_equal_to (current_elem->processor_id(),
462  this->processor_id());
463 #endif // DEBUG
464  neigh = const_cast<RemoteElem *>(remote_elem);
465  }
466  // If neigh and current_elem are more than one level
467  // apart, figuring out whether we have a remote
468  // neighbor here becomes much harder.
469  else if (neigh && (current_elem->subactive() &&
470  neigh->has_children()))
471  {
472  // Find the deepest descendant of neigh which
473  // we could consider for a neighbor. If we run
474  // out of neigh children, then that's our
475  // neighbor. If we find a potential neighbor
476  // with remote_children and we don't find any
477  // potential neighbors among its non-remote
478  // children, then our neighbor must be remote.
479  while (neigh != remote_elem &&
480  neigh->has_children())
481  {
482  bool found_neigh = false;
483  for (unsigned int c = 0;
484  !found_neigh &&
485  c != neigh->n_children(); ++c)
486  {
487  Elem * child = neigh->child_ptr(c);
488  if (child == remote_elem)
489  continue;
490  unsigned int n_neigh = child->n_neighbors();
491  for (unsigned int n=0; n != n_neigh; ++n)
492  {
493  Elem * ncn = child->neighbor(n);
494  if (ncn != remote_elem &&
495  ncn->is_ancestor_of(current_elem))
496  {
497  neigh = ncn;
498  found_neigh = true;
499  break;
500  }
501  }
502  }
503  if (!found_neigh)
504  neigh = const_cast<RemoteElem *>(remote_elem);
505  }
506  }
507  current_elem->set_neighbor(s, neigh);
508 #ifdef DEBUG
509  if (neigh != libmesh_nullptr && neigh != remote_elem)
510  // We ignore subactive elements here because
511  // we don't care about neighbors of subactive element.
512  if ((!neigh->active()) && (!current_elem->subactive()))
513  {
514  libMesh::err << "On processor " << this->processor_id()
515  << std::endl;
516  libMesh::err << "Bad element ID = " << current_elem->id()
517  << ", Side " << s << ", Bad neighbor ID = " << neigh->id() << std::endl;
518  libMesh::err << "Bad element proc_ID = " << current_elem->processor_id()
519  << ", Bad neighbor proc_ID = " << neigh->processor_id() << std::endl;
520  libMesh::err << "Bad element size = " << current_elem->hmin()
521  << ", Bad neighbor size = " << neigh->hmin() << std::endl;
522  libMesh::err << "Bad element center = " << current_elem->centroid()
523  << ", Bad neighbor center = " << neigh->centroid() << std::endl;
524  libMesh::err << "ERROR: "
525  << (current_elem->active()?"Active":"Ancestor")
526  << " Element at level "
527  << current_elem->level() << std::endl;
528  libMesh::err << "with "
529  << (parent->active()?"active":
530  (parent->subactive()?"subactive":"ancestor"))
531  << " parent share "
532  << (neigh->subactive()?"subactive":"ancestor")
533  << " neighbor at level " << neigh->level()
534  << std::endl;
535  NameBasedIO(*this).write ("bad_mesh.gmv");
536  libmesh_error_msg("Problematic mesh written to bad_mesh.gmv.");
537  }
538 #endif // DEBUG
539  }
540  }
541 
542  // We can skip to the next element if we're full-dimension
543  // and therefore don't have any interior parents
544  if (current_elem->dim() >= LIBMESH_DIM)
545  continue;
546 
547  // We have no interior parents unless we can find one later
548  current_elem->set_interior_parent(libmesh_nullptr);
549 
550  Elem * pip = parent->interior_parent();
551 
552  if (!pip)
553  continue;
554 
555  // If there's no interior_parent children, whether due to a
556  // remote element or a non-conformity, then there's no
557  // children to search.
558  if (pip == remote_elem || pip->active())
559  {
560  current_elem->set_interior_parent(pip);
561  continue;
562  }
563 
564  // For node comparisons we'll need a sensible tolerance
565  Real node_tolerance = current_elem->hmin() * TOLERANCE;
566 
567  // Otherwise our interior_parent should be a child of our
568  // parent's interior_parent.
569  for (unsigned int c=0; c != pip->n_children(); ++c)
570  {
571  Elem * child = pip->child_ptr(c);
572 
573  // If we have a remote_elem, that might be our
574  // interior_parent. We'll set it provisionally now and
575  // keep trying to find something better.
576  if (child == remote_elem)
577  {
578  current_elem->set_interior_parent
579  (const_cast<RemoteElem *>(remote_elem));
580  continue;
581  }
582 
583  bool child_contains_our_nodes = true;
584  for (unsigned int n=0; n != current_elem->n_nodes();
585  ++n)
586  {
587  bool child_contains_this_node = false;
588  for (unsigned int cn=0; cn != child->n_nodes();
589  ++cn)
590  if (child->point(cn).absolute_fuzzy_equals
591  (current_elem->point(n), node_tolerance))
592  {
593  child_contains_this_node = true;
594  break;
595  }
596  if (!child_contains_this_node)
597  {
598  child_contains_our_nodes = false;
599  break;
600  }
601  }
602  if (child_contains_our_nodes)
603  {
604  current_elem->set_interior_parent(child);
605  break;
606  }
607  }
608 
609  // We should have found *some* interior_parent at this
610  // point, whether semilocal or remote.
611  libmesh_assert(current_elem->interior_parent());
612  }
613  }
614 
615 #endif // AMR
616 
617 
618 #ifdef DEBUG
620  !reset_remote_elements);
622 #endif
623 }
624 
625 
626 
627 void UnstructuredMesh::read (const std::string & name,
628  void *,
629  bool skip_renumber_nodes_and_elements,
630  bool skip_find_neighbors)
631 {
632  // Set the skip_renumber_nodes_and_elements flag on all processors
633  // if necessary.
634  // This ensures that renumber_nodes_and_elements is *not* called
635  // during prepare_for_use() for certain types of mesh files.
636  // This is required in cases where there is an associated solution
637  // file which expects a certain ordering of the nodes.
638  if (name.rfind(".gmv") + 4 == name.size())
639  this->allow_renumbering(false);
640 
641  NameBasedIO(*this).read(name);
642 
643  if (skip_renumber_nodes_and_elements)
644  {
645  // Use MeshBase::allow_renumbering() yourself instead.
646  libmesh_deprecated();
647  this->allow_renumbering(false);
648  }
649 
650  // Done reading the mesh. Now prepare it for use.
651  this->prepare_for_use(/*skip_renumber (deprecated)*/ false,
652  skip_find_neighbors);
653 }
654 
655 
656 
657 void UnstructuredMesh::write (const std::string & name)
658 {
659  LOG_SCOPE("write()", "Mesh");
660 
661  NameBasedIO(*this).write(name);
662 }
663 
664 
665 
666 void UnstructuredMesh::write (const std::string & name,
667  const std::vector<Number> & v,
668  const std::vector<std::string> & vn)
669 {
670  LOG_SCOPE("write()", "Mesh");
671 
672  NameBasedIO(*this).write_nodal_data(name, v, vn);
673 }
674 
675 
676 
677 
678 
680  const processor_id_type pid) const
681 {
682 
683  // Issue a warning if the number the number of processors
684  // currently available is less that that requested for
685  // partitioning. This is not necessarily an error since
686  // you may run on one processor and still partition the
687  // mesh into several partitions.
688 #ifdef DEBUG
689  if (this->n_processors() < pid)
690  {
691  libMesh::out << "WARNING: You are creating a "
692  << "mesh for a processor id (="
693  << pid
694  << ") greater than "
695  << "the number of processors available for "
696  << "the calculation. (="
697  << this->n_processors()
698  << ")."
699  << std::endl;
700  }
701 #endif
702 
703  // Create iterators to loop over the list of elements
704  // const_active_pid_elem_iterator it(this->elements_begin(), pid);
705  // const const_active_pid_elem_iterator it_end(this->elements_end(), pid);
706 
708  const const_element_iterator it_end = this->active_pid_elements_end(pid);
709 
710  this->create_submesh (pid_mesh, it, it_end);
711 }
712 
713 
714 
715 
716 
717 
718 
721  const const_element_iterator & it_end) const
722 {
723  // Just in case the subdomain_mesh already has some information
724  // in it, get rid of it.
725  new_mesh.clear();
726 
727  // If we're not serial, our submesh isn't either.
728  // There are no remote elements to delete on an empty mesh, but
729  // calling the method to do so marks the mesh as parallel.
730  if (!this->is_serial())
731  new_mesh.delete_remote_elements();
732 
733  // Fail if (*this == new_mesh), we cannot create a submesh inside ourself!
734  // This may happen if the user accidentally passes the original mesh into
735  // this function! We will check this by making sure we did not just
736  // clear ourself.
737  libmesh_assert_not_equal_to (this->n_nodes(), 0);
738  libmesh_assert_not_equal_to (this->n_elem(), 0);
739 
740  // Container to catch boundary IDs handed back by BoundaryInfo
741  std::vector<boundary_id_type> bc_ids;
742 
743  for (; it != it_end; ++it)
744  {
745  const Elem * old_elem = *it;
746 
747  // Add an equivalent element type to the new_mesh.
748  // Copy ids for this element.
749  Elem * new_elem = Elem::build(old_elem->type()).release();
750  new_elem->set_id() = old_elem->id();
751 #ifdef LIBMESH_ENABLE_UNIQUE_ID
752  new_elem->set_unique_id() = old_elem->unique_id();
753 #endif
754  new_elem->subdomain_id() = old_elem->subdomain_id();
755  new_elem->processor_id() = old_elem->processor_id();
756 
757  new_mesh.add_elem (new_elem);
758 
759  libmesh_assert(new_elem);
760 
761  // Loop over the nodes on this element.
762  for (unsigned int n=0; n<old_elem->n_nodes(); n++)
763  {
764  const dof_id_type this_node_id = old_elem->node_id(n);
765 
766  // Add this node to the new mesh if it's not there already
767  if (!new_mesh.query_node_ptr(this_node_id))
768  {
769 #ifdef LIBMESH_ENABLE_UNIQUE_ID
770  Node *newn =
771 #endif
772  new_mesh.add_point (old_elem->point(n),
773  this_node_id,
774  old_elem->node_ptr(n)->processor_id());
775 
776 #ifdef LIBMESH_ENABLE_UNIQUE_ID
777  newn->set_unique_id() = old_elem->node_ptr(n)->unique_id();
778 #endif
779  }
780 
781  // Define this element's connectivity on the new mesh
782  new_elem->set_node(n) = new_mesh.node_ptr(this_node_id);
783  }
784 
785  // Maybe add boundary conditions for this element
786  for (unsigned short s=0; s<old_elem->n_sides(); s++)
787  {
788  this->get_boundary_info().boundary_ids(old_elem, s, bc_ids);
789  new_mesh.get_boundary_info().add_side (new_elem, s, bc_ids);
790  }
791  } // end loop over elements
792 
793  // Prepare the new_mesh for use
794  new_mesh.prepare_for_use(/*skip_renumber =*/false);
795 }
796 
797 
798 
799 #ifdef LIBMESH_ENABLE_AMR
801 {
802  LOG_SCOPE ("contract()", "Mesh");
803 
804  // Flag indicating if this call actually changes the mesh
805  bool mesh_changed = false;
806 
809 
810 #ifdef DEBUG
811  for ( ; in != end; ++in)
812  if (*in != libmesh_nullptr)
813  {
814  Elem * el = *in;
815  libmesh_assert(el->active() || el->subactive() || el->ancestor());
816  }
817  in = elements_begin();
818 #endif
819 
820  // Loop over the elements.
821  for ( ; in != end; ++in)
822  if (*in != libmesh_nullptr)
823  {
824  Elem * el = *in;
825 
826  // Delete all the subactive ones
827  if (el->subactive())
828  {
829  // No level-0 element should be subactive.
830  // Note that we CAN'T test elem->level(), as that
831  // touches elem->parent()->dim(), and elem->parent()
832  // might have already been deleted!
833  libmesh_assert(el->parent());
834 
835  // Delete the element
836  // This just sets a pointer to NULL, and doesn't
837  // invalidate any iterators
838  this->delete_elem(el);
839 
840  // the mesh has certainly changed
841  mesh_changed = true;
842  }
843  else
844  {
845  // Compress all the active ones
846  if (el->active())
847  el->contract();
848  else
849  libmesh_assert (el->ancestor());
850  }
851  }
852 
853  // Strip any newly-created NULL voids out of the element array
855 
856  // FIXME: Need to understand why deleting subactive children
857  // invalidates the point locator. For now we will clear it explicitly
858  this->clear_point_locator();
859 
860  // Allow our GhostingFunctor objects to reinit if necessary.
861  std::set<GhostingFunctor *>::iterator gf_it = this->ghosting_functors_begin();
862  const std::set<GhostingFunctor *>::iterator gf_end = this->ghosting_functors_end();
863  for (; gf_it != gf_end; ++gf_it)
864  {
865  GhostingFunctor *gf = *gf_it;
866  libmesh_assert(gf);
867  gf->mesh_reinit();
868  }
869 
870  return mesh_changed;
871 }
872 #endif // #ifdef LIBMESH_ENABLE_AMR
873 
874 } // namespace libMesh
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
bool ancestor() const
Definition: elem.C:1569
virtual void read(const std::string &mesh_file) libmesh_override
Definition: namebased_io.C:62
bool has_children() const
Definition: elem.h:2097
bool closed()
Definition: libmesh.C:281
const BoundaryInfo & get_boundary_info() const
Definition: mesh_base.h:114
unique_id_type & set_unique_id()
Definition: dof_object.h:654
virtual void reserve_nodes(const dof_id_type nn)=0
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false) libmesh_override
Used by ParallelMesh to represent an Elem owned by another processor.
Definition: remote_elem.h:44
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1795
A geometric point in (x,y,z) space associated with a DOF.
Definition: node.h:52
void create_pid_mesh(UnstructuredMesh &pid_mesh, const processor_id_type pid) const
bool subactive() const
Definition: elem.h:2077
virtual bool is_serial() const
Definition: mesh_base.h:137
void libmesh_assert_valid_amr_interior_parents(const MeshBase &mesh)
Definition: mesh_tools.C:1206
bool active() const
Definition: elem.h:2059
static UniquePtr< Elem > build(const ElemType type, Elem *p=libmesh_nullptr)
Definition: elem.C:234
virtual ElemType type() const =0
virtual element_iterator level_elements_begin(unsigned int level)=0
virtual void copy_nodes_and_elements(const UnstructuredMesh &other_mesh, const bool skip_find_neighbors=false)
void allow_renumbering(bool allow)
Definition: mesh_base.h:734
unsigned int p_level() const
Definition: elem.h:2224
bool skip_partitioning() const
Definition: mesh_base.h:760
void skip_partitioning(bool skip)
Definition: mesh_base.h:759
virtual unsigned int n_neighbors() const
Definition: elem.h:557
const Elem * parent() const
Definition: elem.h:2148
virtual void write_nodal_data(const std::string &, const std::vector< Number > &, const std::vector< std::string > &) libmesh_override
Definition: namebased_io.C:425
processor_id_type n_processors() const
const Elem * interior_parent() const
Definition: elem.C:941
The base class for all geometric element types.
Definition: elem.h:86
virtual Real hmin() const
Definition: elem.C:449
uint8_t processor_id_type
Definition: id_types.h:99
void add_child(Elem *elem)
Definition: elem.C:1604
const class libmesh_nullptr_t libmesh_nullptr
virtual const Node * node_ptr(const dof_id_type i) const =0
static const Real TOLERANCE
IterBase * end
virtual void find_neighbors(const bool reset_remote_elements=false, const bool reset_current_list=true) libmesh_override
virtual bool is_child_on_side(const unsigned int c, const unsigned int s) const =0
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
void set_interior_parent(Elem *p)
Definition: elem.C:993
Base class for Mesh.
Definition: mesh_base.h:67
dof_id_type & set_id()
Definition: dof_object.h:633
libmesh_assert(j)
std::set< GhostingFunctor * >::const_iterator ghosting_functors_end() const
Definition: mesh_base.h:789
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual unsigned int n_nodes() const =0
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:1821
std::vector< boundary_id_type > boundary_ids(const Node *node) const
void create_submesh(UnstructuredMesh &new_mesh, const_element_iterator &it, const const_element_iterator &it_end) const
virtual node_iterator nodes_begin()=0
virtual element_iterator elements_begin()=0
virtual element_iterator level_elements_end(unsigned int level)=0
unsigned int _n_parts
Definition: mesh_base.h:1315
void allow_remote_element_removal(bool allow)
Definition: mesh_base.h:743
virtual void write(const std::string &mesh_file) libmesh_override
Definition: namebased_io.C:279
virtual void delete_elem(Elem *e)=0
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1730
bool is_ancestor_of(const Elem *descendant) const
Definition: elem.h:2127
virtual element_iterator active_pid_elements_begin(processor_id_type proc_id)=0
Base class for Replicated and Distributed meshes.
virtual Elem * add_elem(Elem *e)=0
virtual element_iterator elements_end()=0
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:621
virtual bool contract() libmesh_override
virtual const Node * query_node_ptr(const dof_id_type i) const =0
virtual void write(const std::string &name) libmesh_override
void clear_point_locator()
Definition: mesh_base.C:552
const Elem * child_ptr(unsigned int i) const
Definition: elem.h:2247
void libmesh_assert_valid_amr_elem_ids(const MeshBase &mesh)
Definition: mesh_tools.C:1182
RefinementState p_refinement_flag() const
Definition: elem.h:2320
void prepare_for_use(const bool skip_renumber_nodes_and_elements=false, const bool skip_find_neighbors=false)
Definition: mesh_base.C:173
OStreamProxy err(std::cerr)
subdomain_id_type subdomain_id() const
Definition: elem.h:1805
virtual unsigned int n_sides() const =0
void set_neighbor(const unsigned int i, Elem *n)
Definition: elem.h:1852
virtual void clear()
Definition: mesh_base.C:284
virtual unsigned int n_children() const =0
UnstructuredMesh(const Parallel::Communicator &comm_in, unsigned char dim=1)
unsigned int which_child_am_i(const Elem *e) const
Definition: elem.h:2287
std::set< GhostingFunctor * >::const_iterator ghosting_functors_begin() const
Definition: mesh_base.h:783
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual node_iterator nodes_end()=0
virtual element_iterator active_pid_elements_end(processor_id_type proc_id)=0
const Point & point(const unsigned int i) const
Definition: elem.h:1667
bool absolute_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition: type_vector.h:954
virtual dof_id_type key(const unsigned int s) const =0
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
bool allow_renumbering() const
Definition: mesh_base.h:735
virtual unsigned int dim() const =0
unsigned int level() const
Definition: elem.h:2190
RefinementState refinement_flag() const
Definition: elem.h:2304
bool initialized()
Definition: libmesh.C:274
virtual Point centroid() const
Definition: elem.C:437
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1689
void libmesh_assert_valid_neighbors(const MeshBase &mesh, bool assert_valid_remote_elems=true)
Definition: mesh_tools.C:1825
virtual void delete_remote_elements()
Definition: mesh_base.h:179
dof_id_type id() const
Definition: dof_object.h:624
virtual dof_id_type n_nodes() const =0
virtual const Elem * elem_ptr(const dof_id_type i) const =0
virtual dof_id_type n_elem() const =0
unique_id_type unique_id() const
Definition: dof_object.h:641
bool allow_remote_element_removal() const
Definition: mesh_base.h:744
OStreamProxy out(std::cout)
Elem * neighbor(const unsigned int i) const
Definition: elem.h:1841
virtual void reserve_elem(const dof_id_type ne)=0
virtual UniquePtr< Elem > side_ptr(unsigned int i)=0
virtual void renumber_nodes_and_elements()=0
uint8_t dof_id_type
Definition: id_types.h:64
processor_id_type processor_id() const
processor_id_type processor_id() const
Definition: dof_object.h:686
const RemoteElem * remote_elem
Definition: remote_elem.C:57