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