libMesh
mesh_modification.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 <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 #include <cmath> // for std::acos()
23 #include <algorithm>
24 #include <limits>
25 #include <map>
26 
27 // Local includes
28 #include "libmesh/boundary_info.h"
29 #include "libmesh/function_base.h"
30 #include "libmesh/cell_tet4.h"
31 #include "libmesh/cell_tet10.h"
32 #include "libmesh/face_tri3.h"
33 #include "libmesh/face_tri6.h"
34 #include "libmesh/libmesh_logging.h"
35 #include "libmesh/mesh_communication.h"
36 #include "libmesh/mesh_modification.h"
37 #include "libmesh/mesh_tools.h"
38 #include "libmesh/parallel.h"
39 #include "libmesh/remote_elem.h"
40 #include "libmesh/string_to_enum.h"
41 #include "libmesh/unstructured_mesh.h"
42 #include "libmesh/partitioner.h"
43 
44 namespace
45 {
46 bool split_first_diagonal(const libMesh::Elem * elem,
47  unsigned int diag_1_node_1,
48  unsigned int diag_1_node_2,
49  unsigned int diag_2_node_1,
50  unsigned int diag_2_node_2)
51 {
52  return ((elem->node_id(diag_1_node_1) > elem->node_id(diag_2_node_1) &&
53  elem->node_id(diag_1_node_1) > elem->node_id(diag_2_node_2)) ||
54  (elem->node_id(diag_1_node_2) > elem->node_id(diag_2_node_1) &&
55  elem->node_id(diag_1_node_2) > elem->node_id(diag_2_node_2)));
56 }
57 
58 }
59 
60 namespace libMesh
61 {
62 
63 
64 // ------------------------------------------------------------
65 // MeshTools::Modification functions for mesh modification
67  const Real factor,
68  const bool perturb_boundary)
69 {
70  libmesh_assert (mesh.n_nodes());
71  libmesh_assert (mesh.n_elem());
72  libmesh_assert ((factor >= 0.) && (factor <= 1.));
73 
74  LOG_SCOPE("distort()", "MeshTools::Modification");
75 
76  // First find nodes on the boundary and flag them
77  // so that we don't move them
78  // on_boundary holds false (not on boundary) and true (on boundary)
79  std::vector<bool> on_boundary (mesh.max_node_id(), false);
80 
81  if (!perturb_boundary) MeshTools::find_boundary_nodes (mesh, on_boundary);
82 
83  // Now calculate the minimum distance to
84  // neighboring nodes for each node.
85  // hmin holds these distances.
86  std::vector<float> hmin (mesh.max_node_id(),
88 
89  for (const auto & elem : mesh.active_element_ptr_range())
90  for (auto & n : elem->node_ref_range())
91  hmin[n.id()] = std::min(hmin[n.id()],
92  static_cast<float>(elem->hmin()));
93 
94 
95  // Now actually move the nodes
96  {
97  const unsigned int seed = 123456;
98 
99  // seed the random number generator.
100  // We'll loop from 1 to n_nodes on every processor, even those
101  // that don't have a particular node, so that the pseudorandom
102  // numbers will be the same everywhere.
103  std::srand(seed);
104 
105  // If the node is on the boundary or
106  // the node is not used by any element (hmin[n]<1.e20)
107  // then we should not move it.
108  // [Note: Testing for (in)equality might be wrong
109  // (different types, namely float and double)]
110  for (unsigned int n=0; n<mesh.max_node_id(); n++)
111  if (!on_boundary[n] && (hmin[n] < 1.e20) )
112  {
113  // the direction, random but unit normalized
114 
115  Point dir( static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
116  (mesh.mesh_dimension() > 1) ?
117  static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX)
118  : 0.,
119  ((mesh.mesh_dimension() == 3) ?
120  static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX)
121  : 0.)
122  );
123 
124  dir(0) = (dir(0)-.5)*2.;
125  if (mesh.mesh_dimension() > 1)
126  dir(1) = (dir(1)-.5)*2.;
127  if (mesh.mesh_dimension() == 3)
128  dir(2) = (dir(2)-.5)*2.;
129 
130  dir = dir.unit();
131 
132  Node * node = mesh.query_node_ptr(n);
133  if (!node)
134  continue;
135 
136  (*node)(0) += dir(0)*factor*hmin[n];
137  if (mesh.mesh_dimension() > 1)
138  (*node)(1) += dir(1)*factor*hmin[n];
139  if (mesh.mesh_dimension() == 3)
140  (*node)(2) += dir(2)*factor*hmin[n];
141  }
142  }
143 }
144 
145 
146 
148  const FunctionBase<Real> & mapfunc)
149 {
150  libmesh_assert (mesh.n_nodes());
151  libmesh_assert (mesh.n_elem());
152 
153  LOG_SCOPE("redistribute()", "MeshTools::Modification");
154 
155  DenseVector<Real> output_vec(LIBMESH_DIM);
156 
157  // FIXME - we should thread this later.
158  UniquePtr<FunctionBase<Real>> myfunc = mapfunc.clone();
159 
160  for (auto & node : mesh.node_ptr_range())
161  {
162  (*myfunc)(*node, output_vec);
163 
164  (*node)(0) = output_vec(0);
165 #if LIBMESH_DIM > 1
166  (*node)(1) = output_vec(1);
167 #endif
168 #if LIBMESH_DIM > 2
169  (*node)(2) = output_vec(2);
170 #endif
171  }
172 }
173 
174 
175 
177  const Real xt,
178  const Real yt,
179  const Real zt)
180 {
181  const Point p(xt, yt, zt);
182 
183  for (auto & node : mesh.node_ptr_range())
184  *node += p;
185 }
186 
187 
188 // void MeshTools::Modification::rotate2D (MeshBase & mesh,
189 // const Real alpha)
190 // {
191 // libmesh_assert_not_equal_to (mesh.mesh_dimension(), 1);
192 
193 // const Real pi = std::acos(-1);
194 // const Real a = alpha/180.*pi;
195 // for (unsigned int n=0; n<mesh.n_nodes(); n++)
196 // {
197 // const Point p = mesh.node_ref(n);
198 // const Real x = p(0);
199 // const Real y = p(1);
200 // const Real z = p(2);
201 // mesh.node_ref(n) = Point(std::cos(a)*x - std::sin(a)*y,
202 // std::sin(a)*x + std::cos(a)*y,
203 // z);
204 // }
205 
206 // }
207 
208 
209 
211  const Real phi,
212  const Real theta,
213  const Real psi)
214 {
215 #if LIBMESH_DIM == 3
216  const Real p = -phi/180.*libMesh::pi;
217  const Real t = -theta/180.*libMesh::pi;
218  const Real s = -psi/180.*libMesh::pi;
219  const Real sp = std::sin(p), cp = std::cos(p);
220  const Real st = std::sin(t), ct = std::cos(t);
221  const Real ss = std::sin(s), cs = std::cos(s);
222 
223  // We follow the convention described at http://mathworld.wolfram.com/EulerAngles.html
224  // (equations 6-14 give the entries of the composite transformation matrix).
225  // The rotations are performed sequentially about the z, x, and z axes, in that order.
226  // A positive angle yields a counter-clockwise rotation about the axis in question.
227  for (auto & node : mesh.node_ptr_range())
228  {
229  const Point pt = *node;
230  const Real x = pt(0);
231  const Real y = pt(1);
232  const Real z = pt(2);
233  *node = Point(( cp*cs-sp*ct*ss)*x + ( sp*cs+cp*ct*ss)*y + (st*ss)*z,
234  (-cp*ss-sp*ct*cs)*x + (-sp*ss+cp*ct*cs)*y + (st*cs)*z,
235  ( sp*st)*x + (-cp*st)*y + (ct)*z );
236  }
237 #else
238  libmesh_error_msg("MeshTools::Modification::rotate() requires libMesh to be compiled with LIBMESH_DIM==3");
239 #endif
240 }
241 
242 
244  const Real xs,
245  const Real ys,
246  const Real zs)
247 {
248  const Real x_scale = xs;
249  Real y_scale = ys;
250  Real z_scale = zs;
251 
252  if (ys == 0.)
253  {
254  libmesh_assert_equal_to (zs, 0.);
255 
256  y_scale = z_scale = x_scale;
257  }
258 
259  // Scale the x coordinate in all dimensions
260  for (auto & node : mesh.node_ptr_range())
261  (*node)(0) *= x_scale;
262 
263  // Only scale the y coordinate in 2 and 3D
264  if (LIBMESH_DIM < 2)
265  return;
266 
267  for (auto & node : mesh.node_ptr_range())
268  (*node)(1) *= y_scale;
269 
270  // Only scale the z coordinate in 3D
271  if (LIBMESH_DIM < 3)
272  return;
273 
274  for (auto & node : mesh.node_ptr_range())
275  (*node)(2) *= z_scale;
276 }
277 
278 
279 
280 
281 // ------------------------------------------------------------
282 // UnstructuredMesh class member functions for mesh modification
284 {
285  /*
286  * when the mesh is not prepared,
287  * at least renumber the nodes and
288  * elements, so that the node ids
289  * are correct
290  */
291  if (!this->_is_prepared)
293 
294  START_LOG("all_first_order()", "Mesh");
295 
299  std::vector<bool> node_touched_by_me(this->max_node_id(), false);
300 
306  element_iterator endit = elements_end();
307  for (element_iterator it = elements_begin();
308  it != endit; ++it)
309  {
310  Elem * so_elem = *it;
311 
312  libmesh_assert(so_elem);
313 
314  /*
315  * build the first-order equivalent, add to
316  * the new_elements list.
317  */
318  Elem * lo_elem = Elem::build
320  (so_elem->type()), so_elem->parent()).release();
321 
322  const unsigned short n_sides = so_elem->n_sides();
323 
324  for (unsigned short s=0; s != n_sides; ++s)
325  if (so_elem->neighbor_ptr(s) == remote_elem)
326  lo_elem->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
327 
328 #ifdef LIBMESH_ENABLE_AMR
329  /*
330  * Reset the parent links of any child elements
331  */
332  if (so_elem->has_children())
333  for (unsigned int c = 0, nc = so_elem->n_children(); c != nc; ++c)
334  {
335  Elem * child = so_elem->child_ptr(c);
336  child->set_parent(lo_elem);
337  lo_elem->add_child(child, c);
338  }
339 
340  /*
341  * Reset the child link of any parent element
342  */
343  if (so_elem->parent())
344  {
345  unsigned int c =
346  so_elem->parent()->which_child_am_i(so_elem);
347  lo_elem->parent()->replace_child(lo_elem, c);
348  }
349 
350  /*
351  * Copy as much data to the new element as makes sense
352  */
353  lo_elem->set_p_level(so_elem->p_level());
354  lo_elem->set_refinement_flag(so_elem->refinement_flag());
355  lo_elem->set_p_refinement_flag(so_elem->p_refinement_flag());
356 #endif
357 
358  libmesh_assert_equal_to (lo_elem->n_vertices(), so_elem->n_vertices());
359 
360  /*
361  * By definition the vertices of the linear and
362  * second order element are identically numbered.
363  * transfer these.
364  */
365  for (unsigned int v=0; v < so_elem->n_vertices(); v++)
366  {
367  lo_elem->set_node(v) = so_elem->node_ptr(v);
368  node_touched_by_me[lo_elem->node_id(v)] = true;
369  }
370 
371  /*
372  * find_neighbors relies on remote_elem neighbor links being
373  * properly maintained.
374  */
375  for (unsigned short s=0; s != n_sides; s++)
376  {
377  if (so_elem->neighbor_ptr(s) == remote_elem)
378  lo_elem->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
379  }
380 
388  (this->get_boundary_info(), so_elem, lo_elem);
389 
390  /*
391  * The new first-order element is ready.
392  * Inserting it into the mesh will replace and delete
393  * the second-order element.
394  */
395  lo_elem->set_id(so_elem->id());
396 #ifdef LIBMESH_ENABLE_UNIQUE_ID
397  lo_elem->set_unique_id() = so_elem->unique_id();
398 #endif
399  lo_elem->processor_id() = so_elem->processor_id();
400  lo_elem->subdomain_id() = so_elem->subdomain_id();
401  this->insert_elem(lo_elem);
402  }
403 
404  const MeshBase::node_iterator nd_end = this->nodes_end();
406  while (nd != nd_end)
407  {
408  Node * the_node = *nd;
409  ++nd;
410  if (!node_touched_by_me[the_node->id()])
411  this->delete_node(the_node);
412  }
413 
414  // If crazy people applied boundary info to non-vertices and then
415  // deleted those non-vertices, we should make sure their boundary id
416  // caches are correct.
418 
419  STOP_LOG("all_first_order()", "Mesh");
420 
421  // On hanging nodes that used to also be second order nodes, we
422  // might now have an invalid nodal processor_id()
424 
425  // delete or renumber nodes if desired
426  this->prepare_for_use();
427 }
428 
429 
430 
431 void UnstructuredMesh::all_second_order (const bool full_ordered)
432 {
433  // This function must be run on all processors at once
434  parallel_object_only();
435 
436  /*
437  * when the mesh is not prepared,
438  * at least renumber the nodes and
439  * elements, so that the node ids
440  * are correct
441  */
442  if (!this->_is_prepared)
444 
445  /*
446  * If the mesh is empty
447  * then we have nothing to do
448  */
449  if (!this->n_elem())
450  return;
451 
452  /*
453  * If the mesh is already second order
454  * then we have nothing to do.
455  * We have to test for this in a round-about way to avoid
456  * a bug on distributed parallel meshes with more processors
457  * than elements.
458  */
459  bool already_second_order = false;
460  if (this->elements_begin() != this->elements_end() &&
461  (*(this->elements_begin()))->default_order() != FIRST)
462  already_second_order = true;
463  this->comm().max(already_second_order);
464  if (already_second_order)
465  return;
466 
467  START_LOG("all_second_order()", "Mesh");
468 
469  /*
470  * this map helps in identifying second order
471  * nodes. Namely, a second-order node:
472  * - edge node
473  * - face node
474  * - bubble node
475  * is uniquely defined through a set of adjacent
476  * vertices. This set of adjacent vertices is
477  * used to identify already added higher-order
478  * nodes. We are safe to use node id's since we
479  * make sure that these are correctly numbered.
480  */
481  std::map<std::vector<dof_id_type>, Node *> adj_vertices_to_so_nodes;
482 
483  /*
484  * for speed-up of the \p add_point() method, we
485  * can reserve memory. Guess the number of additional
486  * nodes for different dimensions
487  */
488  switch (this->mesh_dimension())
489  {
490  case 1:
491  /*
492  * in 1D, there can only be order-increase from Edge2
493  * to Edge3. Something like 1/2 of n_nodes() have
494  * to be added
495  */
496  this->reserve_nodes(static_cast<unsigned int>
497  (1.5*static_cast<double>(this->n_nodes())));
498  break;
499 
500  case 2:
501  /*
502  * in 2D, either refine from Tri3 to Tri6 (double the nodes)
503  * or from Quad4 to Quad8 (again, double) or Quad9 (2.25 that much)
504  */
505  this->reserve_nodes(static_cast<unsigned int>
506  (2*static_cast<double>(this->n_nodes())));
507  break;
508 
509 
510  case 3:
511  /*
512  * in 3D, either refine from Tet4 to Tet10 (factor = 2.5) up to
513  * Hex8 to Hex27 (something > 3). Since in 3D there _are_ already
514  * quite some nodes, and since we do not want to overburden the memory by
515  * a too conservative guess, use the lower bound
516  */
517  this->reserve_nodes(static_cast<unsigned int>
518  (2.5*static_cast<double>(this->n_nodes())));
519  break;
520 
521  default:
522  // Hm?
523  libmesh_error_msg("Unknown mesh dimension " << this->mesh_dimension());
524  }
525 
526 
527 
528  /*
529  * form a vector that will hold the node id's of
530  * the vertices that are adjacent to the son-th
531  * second-order node. Pull this outside of the
532  * loop so that silly compilers don't repeatedly
533  * create and destroy the vector.
534  */
535  std::vector<dof_id_type> adjacent_vertices_ids;
536 
544  it = elements_begin(),
545  endit = elements_end();
546 
547  for (; it != endit; ++it)
548  {
549  // the linear-order element
550  Elem * lo_elem = *it;
551 
552  libmesh_assert(lo_elem);
553 
554  // make sure it is linear order
555  if (lo_elem->default_order() != FIRST)
556  libmesh_error_msg("ERROR: This is not a linear element: type=" << lo_elem->type());
557 
558  // this does _not_ work for refined elements
559  libmesh_assert_equal_to (lo_elem->level (), 0);
560 
561  /*
562  * build the second-order equivalent, add to
563  * the new_elements list. Note that this here
564  * is the only point where \p full_ordered
565  * is necessary. The remaining code works well
566  * for either type of second-order equivalent, e.g.
567  * Hex20 or Hex27, as equivalents for Hex8
568  */
569  Elem * so_elem =
571  full_ordered) ).release();
572 
573  libmesh_assert_equal_to (lo_elem->n_vertices(), so_elem->n_vertices());
574 
575 
576  /*
577  * By definition the vertices of the linear and
578  * second order element are identically numbered.
579  * transfer these.
580  */
581  for (unsigned int v=0; v < lo_elem->n_vertices(); v++)
582  so_elem->set_node(v) = lo_elem->node_ptr(v);
583 
584  /*
585  * Now handle the additional mid-side nodes. This
586  * is simply handled through a map that remembers
587  * the already-added nodes. This map maps the global
588  * ids of the vertices (that uniquely define this
589  * higher-order node) to the new node.
590  * Notation: son = second-order node
591  */
592  const unsigned int son_begin = so_elem->n_vertices();
593  const unsigned int son_end = so_elem->n_nodes();
594 
595 
596  for (unsigned int son=son_begin; son<son_end; son++)
597  {
598  const unsigned int n_adjacent_vertices =
599  so_elem->n_second_order_adjacent_vertices(son);
600 
601  adjacent_vertices_ids.resize(n_adjacent_vertices);
602 
603  for (unsigned int v=0; v<n_adjacent_vertices; v++)
604  adjacent_vertices_ids[v] =
605  so_elem->node_id( so_elem->second_order_adjacent_vertex(son,v) );
606 
607  /*
608  * \p adjacent_vertices_ids is now in order of the current
609  * side. sort it, so that comparisons with the
610  * \p adjacent_vertices_ids created through other elements'
611  * sides can match
612  */
613  std::sort(adjacent_vertices_ids.begin(),
614  adjacent_vertices_ids.end());
615 
616 
617  // does this set of vertices already has a mid-node added?
618  std::pair<std::map<std::vector<dof_id_type>, Node *>::iterator,
619  std::map<std::vector<dof_id_type>, Node *>::iterator>
620  pos = adj_vertices_to_so_nodes.equal_range (adjacent_vertices_ids);
621 
622  // no, not added yet
623  if (pos.first == pos.second)
624  {
625  /*
626  * for this set of vertices, there is no
627  * second_order node yet. Add it.
628  *
629  * compute the location of the new node as
630  * the average over the adjacent vertices.
631  */
632  Point new_location = this->point(adjacent_vertices_ids[0]);
633  for (unsigned int v=1; v<n_adjacent_vertices; v++)
634  new_location += this->point(adjacent_vertices_ids[v]);
635 
636  new_location /= static_cast<Real>(n_adjacent_vertices);
637 
638  /* Add the new point to the mesh.
639  * If we are on a serialized mesh, then we're doing this
640  * all in sync, and the node processor_id will be
641  * consistent between processors.
642  * If we are on a distributed mesh, we can fix
643  * inconsistent processor ids later, but only if every
644  * processor gives new nodes a *locally* consistent
645  * processor id, so we'll give the new node the
646  * processor id of an adjacent element for now and then
647  * we'll update that later if appropriate.
648  */
649  Node * so_node = this->add_point
650  (new_location, DofObject::invalid_id,
651  lo_elem->processor_id());
652 
653  /*
654  * insert the new node with its defining vertex
655  * set into the map, and relocate pos to this
656  * new entry, so that the so_elem can use
657  * \p pos for inserting the node
658  */
659  adj_vertices_to_so_nodes.insert(pos.first,
660  std::make_pair(adjacent_vertices_ids,
661  so_node));
662 
663  so_elem->set_node(son) = so_node;
664  }
665  // yes, already added.
666  else
667  {
668  Node * so_node = pos.first->second;
669  libmesh_assert(so_node);
670 
671  so_elem->set_node(son) = so_node;
672 
673  // We need to ensure that the processor who should own a
674  // node *knows* they own the node.
675  if (so_node->processor_id() > lo_elem->processor_id())
676  so_node->processor_id() = lo_elem->processor_id();
677  }
678  }
679 
680  /*
681  * find_neighbors relies on remote_elem neighbor links being
682  * properly maintained.
683  */
684  for (auto s : lo_elem->side_index_range())
685  {
686  if (lo_elem->neighbor_ptr(s) == remote_elem)
687  so_elem->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
688  }
689 
702  (this->get_boundary_info(), lo_elem, so_elem);
703 
704  /*
705  * The new second-order element is ready.
706  * Inserting it into the mesh will replace and delete
707  * the first-order element.
708  */
709  so_elem->set_id(lo_elem->id());
710 #ifdef LIBMESH_ENABLE_UNIQUE_ID
711  so_elem->set_unique_id() = lo_elem->unique_id();
712 #endif
713  so_elem->processor_id() = lo_elem->processor_id();
714  so_elem->subdomain_id() = lo_elem->subdomain_id();
715  this->insert_elem(so_elem);
716  }
717 
718  // we can clear the map
719  adj_vertices_to_so_nodes.clear();
720 
721 
722  STOP_LOG("all_second_order()", "Mesh");
723 
724  // In a DistributedMesh our ghost node processor ids may be bad,
725  // the ids of nodes touching remote elements may be inconsistent,
726  // and unique_ids of newly added non-local nodes remain unset.
727  // make_nodes_parallel_consistent() will fix all this.
728  if (!this->is_serial())
730 
731  // renumber nodes, elements etc
732  this->prepare_for_use(/*skip_renumber =*/ false);
733 }
734 
735 
736 
737 
738 
739 
741 {
742  // The number of elements in the original mesh before any additions
743  // or deletions.
744  const dof_id_type n_orig_elem = mesh.n_elem();
745  const dof_id_type max_orig_id = mesh.max_elem_id();
746 
747  // We store pointers to the newly created elements in a vector
748  // until they are ready to be added to the mesh. This is because
749  // adding new elements on the fly can cause reallocation and invalidation
750  // of existing mesh element_iterators.
751  std::vector<Elem *> new_elements;
752 
753  unsigned int max_subelems = 1; // in 1D nothing needs to change
754  if (mesh.mesh_dimension() == 2) // in 2D quads can split into 2 tris
755  max_subelems = 2;
756  if (mesh.mesh_dimension() == 3) // in 3D hexes can split into 6 tets
757  max_subelems = 6;
758 
759  new_elements.reserve (max_subelems*n_orig_elem);
760 
761  // If the original mesh has *side* boundary data, we carry that over
762  // to the new mesh with triangular elements. We currently only
763  // support bringing over side-based BCs to the all-tri mesh, but
764  // that could probably be extended to node and edge-based BCs as
765  // well.
766  const bool mesh_has_boundary_data = (mesh.get_boundary_info().n_boundary_conds() > 0);
767 
768  // Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
769  std::vector<Elem *> new_bndry_elements;
770  std::vector<unsigned short int> new_bndry_sides;
771  std::vector<boundary_id_type> new_bndry_ids;
772 
773  // We may need to add new points if we run into a 1.5th order
774  // element; if we do that on a DistributedMesh in a ghost element then
775  // we will need to fix their ids / unique_ids
776  bool added_new_ghost_point = false;
777 
778  // Iterate over the elements, splitting:
779  // QUADs into pairs of conforming triangles
780  // PYRAMIDs into pairs of conforming tets,
781  // PRISMs into triplets of conforming tets, and
782  // HEXs into quintets or sextets of conforming tets.
783  // We split on the shortest diagonal to give us better
784  // triangle quality in 2D, and we split based on node ids
785  // to guarantee consistency in 3D.
786 
787  // FIXME: This algorithm does not work on refined grids!
788  {
789 #ifdef LIBMESH_ENABLE_UNIQUE_ID
790  unique_id_type max_unique_id = mesh.parallel_max_unique_id();
791 #endif
792 
793  for (auto & elem : mesh.element_ptr_range())
794  {
795  const ElemType etype = elem->type();
796 
797  // all_tri currently only works on coarse meshes
798  libmesh_assert (!elem->parent());
799 
800  // The new elements we will split the quad into.
801  // In 3D we may need as many as 6 tets per hex
802  Elem * subelem[6];
803 
804  for (unsigned int i = 0; i != max_subelems; ++i)
805  subelem[i] = libmesh_nullptr;
806 
807  switch (etype)
808  {
809  case QUAD4:
810  {
811  subelem[0] = new Tri3;
812  subelem[1] = new Tri3;
813 
814  // Check for possible edge swap
815  if ((elem->point(0) - elem->point(2)).norm() <
816  (elem->point(1) - elem->point(3)).norm())
817  {
818  subelem[0]->set_node(0) = elem->node_ptr(0);
819  subelem[0]->set_node(1) = elem->node_ptr(1);
820  subelem[0]->set_node(2) = elem->node_ptr(2);
821 
822  subelem[1]->set_node(0) = elem->node_ptr(0);
823  subelem[1]->set_node(1) = elem->node_ptr(2);
824  subelem[1]->set_node(2) = elem->node_ptr(3);
825  }
826 
827  else
828  {
829  subelem[0]->set_node(0) = elem->node_ptr(0);
830  subelem[0]->set_node(1) = elem->node_ptr(1);
831  subelem[0]->set_node(2) = elem->node_ptr(3);
832 
833  subelem[1]->set_node(0) = elem->node_ptr(1);
834  subelem[1]->set_node(1) = elem->node_ptr(2);
835  subelem[1]->set_node(2) = elem->node_ptr(3);
836  }
837 
838 
839  break;
840  }
841 
842  case QUAD8:
843  {
844  if (elem->processor_id() != mesh.processor_id())
845  added_new_ghost_point = true;
846 
847  subelem[0] = new Tri6;
848  subelem[1] = new Tri6;
849 
850  // Add a new node at the center (vertex average) of the element.
851  Node * new_node = mesh.add_point((mesh.point(elem->node_id(0)) +
852  mesh.point(elem->node_id(1)) +
853  mesh.point(elem->node_id(2)) +
854  mesh.point(elem->node_id(3)))/4,
856  elem->processor_id());
857 
858  // Check for possible edge swap
859  if ((elem->point(0) - elem->point(2)).norm() <
860  (elem->point(1) - elem->point(3)).norm())
861  {
862  subelem[0]->set_node(0) = elem->node_ptr(0);
863  subelem[0]->set_node(1) = elem->node_ptr(1);
864  subelem[0]->set_node(2) = elem->node_ptr(2);
865  subelem[0]->set_node(3) = elem->node_ptr(4);
866  subelem[0]->set_node(4) = elem->node_ptr(5);
867  subelem[0]->set_node(5) = new_node;
868 
869  subelem[1]->set_node(0) = elem->node_ptr(0);
870  subelem[1]->set_node(1) = elem->node_ptr(2);
871  subelem[1]->set_node(2) = elem->node_ptr(3);
872  subelem[1]->set_node(3) = new_node;
873  subelem[1]->set_node(4) = elem->node_ptr(6);
874  subelem[1]->set_node(5) = elem->node_ptr(7);
875 
876  }
877 
878  else
879  {
880  subelem[0]->set_node(0) = elem->node_ptr(3);
881  subelem[0]->set_node(1) = elem->node_ptr(0);
882  subelem[0]->set_node(2) = elem->node_ptr(1);
883  subelem[0]->set_node(3) = elem->node_ptr(7);
884  subelem[0]->set_node(4) = elem->node_ptr(4);
885  subelem[0]->set_node(5) = new_node;
886 
887  subelem[1]->set_node(0) = elem->node_ptr(1);
888  subelem[1]->set_node(1) = elem->node_ptr(2);
889  subelem[1]->set_node(2) = elem->node_ptr(3);
890  subelem[1]->set_node(3) = elem->node_ptr(5);
891  subelem[1]->set_node(4) = elem->node_ptr(6);
892  subelem[1]->set_node(5) = new_node;
893  }
894 
895  break;
896  }
897 
898  case QUAD9:
899  {
900  subelem[0] = new Tri6;
901  subelem[1] = new Tri6;
902 
903  // Check for possible edge swap
904  if ((elem->point(0) - elem->point(2)).norm() <
905  (elem->point(1) - elem->point(3)).norm())
906  {
907  subelem[0]->set_node(0) = elem->node_ptr(0);
908  subelem[0]->set_node(1) = elem->node_ptr(1);
909  subelem[0]->set_node(2) = elem->node_ptr(2);
910  subelem[0]->set_node(3) = elem->node_ptr(4);
911  subelem[0]->set_node(4) = elem->node_ptr(5);
912  subelem[0]->set_node(5) = elem->node_ptr(8);
913 
914  subelem[1]->set_node(0) = elem->node_ptr(0);
915  subelem[1]->set_node(1) = elem->node_ptr(2);
916  subelem[1]->set_node(2) = elem->node_ptr(3);
917  subelem[1]->set_node(3) = elem->node_ptr(8);
918  subelem[1]->set_node(4) = elem->node_ptr(6);
919  subelem[1]->set_node(5) = elem->node_ptr(7);
920  }
921 
922  else
923  {
924  subelem[0]->set_node(0) = elem->node_ptr(0);
925  subelem[0]->set_node(1) = elem->node_ptr(1);
926  subelem[0]->set_node(2) = elem->node_ptr(3);
927  subelem[0]->set_node(3) = elem->node_ptr(4);
928  subelem[0]->set_node(4) = elem->node_ptr(8);
929  subelem[0]->set_node(5) = elem->node_ptr(7);
930 
931  subelem[1]->set_node(0) = elem->node_ptr(1);
932  subelem[1]->set_node(1) = elem->node_ptr(2);
933  subelem[1]->set_node(2) = elem->node_ptr(3);
934  subelem[1]->set_node(3) = elem->node_ptr(5);
935  subelem[1]->set_node(4) = elem->node_ptr(6);
936  subelem[1]->set_node(5) = elem->node_ptr(8);
937  }
938 
939  break;
940  }
941 
942  case PRISM6:
943  {
944  // Prisms all split into three tetrahedra
945  subelem[0] = new Tet4;
946  subelem[1] = new Tet4;
947  subelem[2] = new Tet4;
948 
949  // Triangular faces are not split.
950 
951  // On quad faces, we choose the node with the highest
952  // global id, and we split on the diagonal which
953  // includes that node. This ensures that (even in
954  // parallel, even on distributed meshes) the same
955  // diagonal split will be chosen for elements on either
956  // side of the same quad face. It also ensures that we
957  // always have a mix of "clockwise" and
958  // "counterclockwise" split faces (two of one and one
959  // of the other on each prism; this is useful since the
960  // alternative all-clockwise or all-counterclockwise
961  // face splittings can't be turned into tets without
962  // adding more nodes
963 
964  // Split on 0-4 diagonal
965  if (split_first_diagonal(elem, 0,4, 1,3))
966  {
967  // Split on 0-5 diagonal
968  if (split_first_diagonal(elem, 0,5, 2,3))
969  {
970  // Split on 1-5 diagonal
971  if (split_first_diagonal(elem, 1,5, 2,4))
972  {
973  subelem[0]->set_node(0) = elem->node_ptr(0);
974  subelem[0]->set_node(1) = elem->node_ptr(4);
975  subelem[0]->set_node(2) = elem->node_ptr(5);
976  subelem[0]->set_node(3) = elem->node_ptr(3);
977 
978  subelem[1]->set_node(0) = elem->node_ptr(0);
979  subelem[1]->set_node(1) = elem->node_ptr(4);
980  subelem[1]->set_node(2) = elem->node_ptr(1);
981  subelem[1]->set_node(3) = elem->node_ptr(5);
982 
983  subelem[2]->set_node(0) = elem->node_ptr(0);
984  subelem[2]->set_node(1) = elem->node_ptr(1);
985  subelem[2]->set_node(2) = elem->node_ptr(2);
986  subelem[2]->set_node(3) = elem->node_ptr(5);
987  }
988  else // Split on 2-4 diagonal
989  {
990  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
991 
992  subelem[0]->set_node(0) = elem->node_ptr(0);
993  subelem[0]->set_node(1) = elem->node_ptr(4);
994  subelem[0]->set_node(2) = elem->node_ptr(5);
995  subelem[0]->set_node(3) = elem->node_ptr(3);
996 
997  subelem[1]->set_node(0) = elem->node_ptr(0);
998  subelem[1]->set_node(1) = elem->node_ptr(4);
999  subelem[1]->set_node(2) = elem->node_ptr(2);
1000  subelem[1]->set_node(3) = elem->node_ptr(5);
1001 
1002  subelem[2]->set_node(0) = elem->node_ptr(0);
1003  subelem[2]->set_node(1) = elem->node_ptr(1);
1004  subelem[2]->set_node(2) = elem->node_ptr(2);
1005  subelem[2]->set_node(3) = elem->node_ptr(4);
1006  }
1007  }
1008  else // Split on 2-3 diagonal
1009  {
1010  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1011 
1012  // 0-4 and 2-3 split implies 2-4 split
1013  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1014 
1015  subelem[0]->set_node(0) = elem->node_ptr(0);
1016  subelem[0]->set_node(1) = elem->node_ptr(4);
1017  subelem[0]->set_node(2) = elem->node_ptr(2);
1018  subelem[0]->set_node(3) = elem->node_ptr(3);
1019 
1020  subelem[1]->set_node(0) = elem->node_ptr(3);
1021  subelem[1]->set_node(1) = elem->node_ptr(4);
1022  subelem[1]->set_node(2) = elem->node_ptr(2);
1023  subelem[1]->set_node(3) = elem->node_ptr(5);
1024 
1025  subelem[2]->set_node(0) = elem->node_ptr(0);
1026  subelem[2]->set_node(1) = elem->node_ptr(1);
1027  subelem[2]->set_node(2) = elem->node_ptr(2);
1028  subelem[2]->set_node(3) = elem->node_ptr(4);
1029  }
1030  }
1031  else // Split on 1-3 diagonal
1032  {
1033  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
1034 
1035  // Split on 0-5 diagonal
1036  if (split_first_diagonal(elem, 0,5, 2,3))
1037  {
1038  // 1-3 and 0-5 split implies 1-5 split
1039  libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
1040 
1041  subelem[0]->set_node(0) = elem->node_ptr(1);
1042  subelem[0]->set_node(1) = elem->node_ptr(3);
1043  subelem[0]->set_node(2) = elem->node_ptr(4);
1044  subelem[0]->set_node(3) = elem->node_ptr(5);
1045 
1046  subelem[1]->set_node(0) = elem->node_ptr(1);
1047  subelem[1]->set_node(1) = elem->node_ptr(0);
1048  subelem[1]->set_node(2) = elem->node_ptr(3);
1049  subelem[1]->set_node(3) = elem->node_ptr(5);
1050 
1051  subelem[2]->set_node(0) = elem->node_ptr(0);
1052  subelem[2]->set_node(1) = elem->node_ptr(1);
1053  subelem[2]->set_node(2) = elem->node_ptr(2);
1054  subelem[2]->set_node(3) = elem->node_ptr(5);
1055  }
1056  else // Split on 2-3 diagonal
1057  {
1058  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1059 
1060  // Split on 1-5 diagonal
1061  if (split_first_diagonal(elem, 1,5, 2,4))
1062  {
1063  subelem[0]->set_node(0) = elem->node_ptr(0);
1064  subelem[0]->set_node(1) = elem->node_ptr(1);
1065  subelem[0]->set_node(2) = elem->node_ptr(2);
1066  subelem[0]->set_node(3) = elem->node_ptr(3);
1067 
1068  subelem[1]->set_node(0) = elem->node_ptr(3);
1069  subelem[1]->set_node(1) = elem->node_ptr(1);
1070  subelem[1]->set_node(2) = elem->node_ptr(2);
1071  subelem[1]->set_node(3) = elem->node_ptr(5);
1072 
1073  subelem[2]->set_node(0) = elem->node_ptr(1);
1074  subelem[2]->set_node(1) = elem->node_ptr(3);
1075  subelem[2]->set_node(2) = elem->node_ptr(4);
1076  subelem[2]->set_node(3) = elem->node_ptr(5);
1077  }
1078  else // Split on 2-4 diagonal
1079  {
1080  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1081 
1082  subelem[0]->set_node(0) = elem->node_ptr(0);
1083  subelem[0]->set_node(1) = elem->node_ptr(1);
1084  subelem[0]->set_node(2) = elem->node_ptr(2);
1085  subelem[0]->set_node(3) = elem->node_ptr(3);
1086 
1087  subelem[1]->set_node(0) = elem->node_ptr(2);
1088  subelem[1]->set_node(1) = elem->node_ptr(3);
1089  subelem[1]->set_node(2) = elem->node_ptr(4);
1090  subelem[1]->set_node(3) = elem->node_ptr(5);
1091 
1092  subelem[2]->set_node(0) = elem->node_ptr(3);
1093  subelem[2]->set_node(1) = elem->node_ptr(1);
1094  subelem[2]->set_node(2) = elem->node_ptr(2);
1095  subelem[2]->set_node(3) = elem->node_ptr(4);
1096  }
1097  }
1098  }
1099 
1100  break;
1101  }
1102 
1103  case PRISM18:
1104  {
1105  subelem[0] = new Tet10;
1106  subelem[1] = new Tet10;
1107  subelem[2] = new Tet10;
1108 
1109  // Split on 0-4 diagonal
1110  if (split_first_diagonal(elem, 0,4, 1,3))
1111  {
1112  // Split on 0-5 diagonal
1113  if (split_first_diagonal(elem, 0,5, 2,3))
1114  {
1115  // Split on 1-5 diagonal
1116  if (split_first_diagonal(elem, 1,5, 2,4))
1117  {
1118  subelem[0]->set_node(0) = elem->node_ptr(0);
1119  subelem[0]->set_node(1) = elem->node_ptr(4);
1120  subelem[0]->set_node(2) = elem->node_ptr(5);
1121  subelem[0]->set_node(3) = elem->node_ptr(3);
1122 
1123  subelem[0]->set_node(4) = elem->node_ptr(15);
1124  subelem[0]->set_node(5) = elem->node_ptr(13);
1125  subelem[0]->set_node(6) = elem->node_ptr(17);
1126  subelem[0]->set_node(7) = elem->node_ptr(9);
1127  subelem[0]->set_node(8) = elem->node_ptr(12);
1128  subelem[0]->set_node(9) = elem->node_ptr(14);
1129 
1130  subelem[1]->set_node(0) = elem->node_ptr(0);
1131  subelem[1]->set_node(1) = elem->node_ptr(4);
1132  subelem[1]->set_node(2) = elem->node_ptr(1);
1133  subelem[1]->set_node(3) = elem->node_ptr(5);
1134 
1135  subelem[1]->set_node(4) = elem->node_ptr(15);
1136  subelem[1]->set_node(5) = elem->node_ptr(10);
1137  subelem[1]->set_node(6) = elem->node_ptr(6);
1138  subelem[1]->set_node(7) = elem->node_ptr(17);
1139  subelem[1]->set_node(8) = elem->node_ptr(13);
1140  subelem[1]->set_node(9) = elem->node_ptr(16);
1141 
1142  subelem[2]->set_node(0) = elem->node_ptr(0);
1143  subelem[2]->set_node(1) = elem->node_ptr(1);
1144  subelem[2]->set_node(2) = elem->node_ptr(2);
1145  subelem[2]->set_node(3) = elem->node_ptr(5);
1146 
1147  subelem[2]->set_node(4) = elem->node_ptr(6);
1148  subelem[2]->set_node(5) = elem->node_ptr(7);
1149  subelem[2]->set_node(6) = elem->node_ptr(8);
1150  subelem[2]->set_node(7) = elem->node_ptr(17);
1151  subelem[2]->set_node(8) = elem->node_ptr(16);
1152  subelem[2]->set_node(9) = elem->node_ptr(11);
1153  }
1154  else // Split on 2-4 diagonal
1155  {
1156  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1157 
1158  subelem[0]->set_node(0) = elem->node_ptr(0);
1159  subelem[0]->set_node(1) = elem->node_ptr(4);
1160  subelem[0]->set_node(2) = elem->node_ptr(5);
1161  subelem[0]->set_node(3) = elem->node_ptr(3);
1162 
1163  subelem[0]->set_node(4) = elem->node_ptr(15);
1164  subelem[0]->set_node(5) = elem->node_ptr(13);
1165  subelem[0]->set_node(6) = elem->node_ptr(17);
1166  subelem[0]->set_node(7) = elem->node_ptr(9);
1167  subelem[0]->set_node(8) = elem->node_ptr(12);
1168  subelem[0]->set_node(9) = elem->node_ptr(14);
1169 
1170  subelem[1]->set_node(0) = elem->node_ptr(0);
1171  subelem[1]->set_node(1) = elem->node_ptr(4);
1172  subelem[1]->set_node(2) = elem->node_ptr(2);
1173  subelem[1]->set_node(3) = elem->node_ptr(5);
1174 
1175  subelem[1]->set_node(4) = elem->node_ptr(15);
1176  subelem[1]->set_node(5) = elem->node_ptr(16);
1177  subelem[1]->set_node(6) = elem->node_ptr(8);
1178  subelem[1]->set_node(7) = elem->node_ptr(17);
1179  subelem[1]->set_node(8) = elem->node_ptr(13);
1180  subelem[1]->set_node(9) = elem->node_ptr(11);
1181 
1182  subelem[2]->set_node(0) = elem->node_ptr(0);
1183  subelem[2]->set_node(1) = elem->node_ptr(1);
1184  subelem[2]->set_node(2) = elem->node_ptr(2);
1185  subelem[2]->set_node(3) = elem->node_ptr(4);
1186 
1187  subelem[2]->set_node(4) = elem->node_ptr(6);
1188  subelem[2]->set_node(5) = elem->node_ptr(7);
1189  subelem[2]->set_node(6) = elem->node_ptr(8);
1190  subelem[2]->set_node(7) = elem->node_ptr(15);
1191  subelem[2]->set_node(8) = elem->node_ptr(10);
1192  subelem[2]->set_node(9) = elem->node_ptr(16);
1193  }
1194  }
1195  else // Split on 2-3 diagonal
1196  {
1197  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1198 
1199  // 0-4 and 2-3 split implies 2-4 split
1200  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1201 
1202  subelem[0]->set_node(0) = elem->node_ptr(0);
1203  subelem[0]->set_node(1) = elem->node_ptr(4);
1204  subelem[0]->set_node(2) = elem->node_ptr(2);
1205  subelem[0]->set_node(3) = elem->node_ptr(3);
1206 
1207  subelem[0]->set_node(4) = elem->node_ptr(15);
1208  subelem[0]->set_node(5) = elem->node_ptr(16);
1209  subelem[0]->set_node(6) = elem->node_ptr(8);
1210  subelem[0]->set_node(7) = elem->node_ptr(9);
1211  subelem[0]->set_node(8) = elem->node_ptr(12);
1212  subelem[0]->set_node(9) = elem->node_ptr(17);
1213 
1214  subelem[1]->set_node(0) = elem->node_ptr(3);
1215  subelem[1]->set_node(1) = elem->node_ptr(4);
1216  subelem[1]->set_node(2) = elem->node_ptr(2);
1217  subelem[1]->set_node(3) = elem->node_ptr(5);
1218 
1219  subelem[1]->set_node(4) = elem->node_ptr(12);
1220  subelem[1]->set_node(5) = elem->node_ptr(16);
1221  subelem[1]->set_node(6) = elem->node_ptr(17);
1222  subelem[1]->set_node(7) = elem->node_ptr(14);
1223  subelem[1]->set_node(8) = elem->node_ptr(13);
1224  subelem[1]->set_node(9) = elem->node_ptr(11);
1225 
1226  subelem[2]->set_node(0) = elem->node_ptr(0);
1227  subelem[2]->set_node(1) = elem->node_ptr(1);
1228  subelem[2]->set_node(2) = elem->node_ptr(2);
1229  subelem[2]->set_node(3) = elem->node_ptr(4);
1230 
1231  subelem[2]->set_node(4) = elem->node_ptr(6);
1232  subelem[2]->set_node(5) = elem->node_ptr(7);
1233  subelem[2]->set_node(6) = elem->node_ptr(8);
1234  subelem[2]->set_node(7) = elem->node_ptr(15);
1235  subelem[2]->set_node(8) = elem->node_ptr(10);
1236  subelem[2]->set_node(9) = elem->node_ptr(16);
1237  }
1238  }
1239  else // Split on 1-3 diagonal
1240  {
1241  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
1242 
1243  // Split on 0-5 diagonal
1244  if (split_first_diagonal(elem, 0,5, 2,3))
1245  {
1246  // 1-3 and 0-5 split implies 1-5 split
1247  libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
1248 
1249  subelem[0]->set_node(0) = elem->node_ptr(1);
1250  subelem[0]->set_node(1) = elem->node_ptr(3);
1251  subelem[0]->set_node(2) = elem->node_ptr(4);
1252  subelem[0]->set_node(3) = elem->node_ptr(5);
1253 
1254  subelem[0]->set_node(4) = elem->node_ptr(15);
1255  subelem[0]->set_node(5) = elem->node_ptr(12);
1256  subelem[0]->set_node(6) = elem->node_ptr(10);
1257  subelem[0]->set_node(7) = elem->node_ptr(16);
1258  subelem[0]->set_node(8) = elem->node_ptr(14);
1259  subelem[0]->set_node(9) = elem->node_ptr(13);
1260 
1261  subelem[1]->set_node(0) = elem->node_ptr(1);
1262  subelem[1]->set_node(1) = elem->node_ptr(0);
1263  subelem[1]->set_node(2) = elem->node_ptr(3);
1264  subelem[1]->set_node(3) = elem->node_ptr(5);
1265 
1266  subelem[1]->set_node(4) = elem->node_ptr(6);
1267  subelem[1]->set_node(5) = elem->node_ptr(9);
1268  subelem[1]->set_node(6) = elem->node_ptr(15);
1269  subelem[1]->set_node(7) = elem->node_ptr(16);
1270  subelem[1]->set_node(8) = elem->node_ptr(17);
1271  subelem[1]->set_node(9) = elem->node_ptr(14);
1272 
1273  subelem[2]->set_node(0) = elem->node_ptr(0);
1274  subelem[2]->set_node(1) = elem->node_ptr(1);
1275  subelem[2]->set_node(2) = elem->node_ptr(2);
1276  subelem[2]->set_node(3) = elem->node_ptr(5);
1277 
1278  subelem[2]->set_node(4) = elem->node_ptr(6);
1279  subelem[2]->set_node(5) = elem->node_ptr(7);
1280  subelem[2]->set_node(6) = elem->node_ptr(8);
1281  subelem[2]->set_node(7) = elem->node_ptr(17);
1282  subelem[2]->set_node(8) = elem->node_ptr(16);
1283  subelem[2]->set_node(9) = elem->node_ptr(11);
1284  }
1285  else // Split on 2-3 diagonal
1286  {
1287  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1288 
1289  // Split on 1-5 diagonal
1290  if (split_first_diagonal(elem, 1,5, 2,4))
1291  {
1292  subelem[0]->set_node(0) = elem->node_ptr(0);
1293  subelem[0]->set_node(1) = elem->node_ptr(1);
1294  subelem[0]->set_node(2) = elem->node_ptr(2);
1295  subelem[0]->set_node(3) = elem->node_ptr(3);
1296 
1297  subelem[0]->set_node(4) = elem->node_ptr(6);
1298  subelem[0]->set_node(5) = elem->node_ptr(7);
1299  subelem[0]->set_node(6) = elem->node_ptr(8);
1300  subelem[0]->set_node(7) = elem->node_ptr(9);
1301  subelem[0]->set_node(8) = elem->node_ptr(15);
1302  subelem[0]->set_node(9) = elem->node_ptr(17);
1303 
1304  subelem[1]->set_node(0) = elem->node_ptr(3);
1305  subelem[1]->set_node(1) = elem->node_ptr(1);
1306  subelem[1]->set_node(2) = elem->node_ptr(2);
1307  subelem[1]->set_node(3) = elem->node_ptr(5);
1308 
1309  subelem[1]->set_node(4) = elem->node_ptr(15);
1310  subelem[1]->set_node(5) = elem->node_ptr(7);
1311  subelem[1]->set_node(6) = elem->node_ptr(17);
1312  subelem[1]->set_node(7) = elem->node_ptr(14);
1313  subelem[1]->set_node(8) = elem->node_ptr(16);
1314  subelem[1]->set_node(9) = elem->node_ptr(11);
1315 
1316  subelem[2]->set_node(0) = elem->node_ptr(1);
1317  subelem[2]->set_node(1) = elem->node_ptr(3);
1318  subelem[2]->set_node(2) = elem->node_ptr(4);
1319  subelem[2]->set_node(3) = elem->node_ptr(5);
1320 
1321  subelem[2]->set_node(4) = elem->node_ptr(15);
1322  subelem[2]->set_node(5) = elem->node_ptr(12);
1323  subelem[2]->set_node(6) = elem->node_ptr(10);
1324  subelem[2]->set_node(7) = elem->node_ptr(16);
1325  subelem[2]->set_node(8) = elem->node_ptr(14);
1326  subelem[2]->set_node(9) = elem->node_ptr(13);
1327  }
1328  else // Split on 2-4 diagonal
1329  {
1330  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1331 
1332  subelem[0]->set_node(0) = elem->node_ptr(0);
1333  subelem[0]->set_node(1) = elem->node_ptr(1);
1334  subelem[0]->set_node(2) = elem->node_ptr(2);
1335  subelem[0]->set_node(3) = elem->node_ptr(3);
1336 
1337  subelem[0]->set_node(4) = elem->node_ptr(6);
1338  subelem[0]->set_node(5) = elem->node_ptr(7);
1339  subelem[0]->set_node(6) = elem->node_ptr(8);
1340  subelem[0]->set_node(7) = elem->node_ptr(9);
1341  subelem[0]->set_node(8) = elem->node_ptr(15);
1342  subelem[0]->set_node(9) = elem->node_ptr(17);
1343 
1344  subelem[1]->set_node(0) = elem->node_ptr(2);
1345  subelem[1]->set_node(1) = elem->node_ptr(3);
1346  subelem[1]->set_node(2) = elem->node_ptr(4);
1347  subelem[1]->set_node(3) = elem->node_ptr(5);
1348 
1349  subelem[1]->set_node(4) = elem->node_ptr(17);
1350  subelem[1]->set_node(5) = elem->node_ptr(12);
1351  subelem[1]->set_node(6) = elem->node_ptr(16);
1352  subelem[1]->set_node(7) = elem->node_ptr(11);
1353  subelem[1]->set_node(8) = elem->node_ptr(14);
1354  subelem[1]->set_node(9) = elem->node_ptr(13);
1355 
1356  subelem[2]->set_node(0) = elem->node_ptr(3);
1357  subelem[2]->set_node(1) = elem->node_ptr(1);
1358  subelem[2]->set_node(2) = elem->node_ptr(2);
1359  subelem[2]->set_node(3) = elem->node_ptr(4);
1360 
1361  subelem[2]->set_node(4) = elem->node_ptr(15);
1362  subelem[2]->set_node(5) = elem->node_ptr(7);
1363  subelem[2]->set_node(6) = elem->node_ptr(17);
1364  subelem[2]->set_node(7) = elem->node_ptr(12);
1365  subelem[2]->set_node(8) = elem->node_ptr(10);
1366  subelem[2]->set_node(9) = elem->node_ptr(16);
1367  }
1368  }
1369  }
1370 
1371  break;
1372  }
1373 
1374  // No need to split elements that are already simplicial:
1375  case EDGE2:
1376  case EDGE3:
1377  case EDGE4:
1378  case TRI3:
1379  case TRI6:
1380  case TET4:
1381  case TET10:
1382  case INFEDGE2:
1383  // No way to split infinite quad/prism elements, so
1384  // hopefully no need to
1385  case INFQUAD4:
1386  case INFQUAD6:
1387  case INFPRISM6:
1388  case INFPRISM12:
1389  continue;
1390  // If we're left with an unimplemented hex we're probably
1391  // out of luck. TODO: implement hexes
1392  default:
1393  {
1394  libMesh::err << "Error, encountered unimplemented element "
1395  << Utility::enum_to_string<ElemType>(etype)
1396  << " in MeshTools::Modification::all_tri()..."
1397  << std::endl;
1398  libmesh_not_implemented();
1399  }
1400  } // end switch (etype)
1401 
1402 
1403 
1404  // Be sure the correct IDs are also set for all subelems.
1405  for (unsigned int i=0; i != max_subelems; ++i)
1406  if (subelem[i]) {
1407  subelem[i]->processor_id() = elem->processor_id();
1408  subelem[i]->subdomain_id() = elem->subdomain_id();
1409  }
1410 
1411  // On a mesh with boundary data, we need to move that data to
1412  // the new elements.
1413 
1414  // On a mesh which is distributed, we need to move
1415  // remote_elem links to the new elements.
1416  bool mesh_is_serial = mesh.is_serial();
1417 
1418  if (mesh_has_boundary_data || mesh_is_serial)
1419  {
1420  // Container to key boundary IDs handed back by the BoundaryInfo object.
1421  std::vector<boundary_id_type> bc_ids;
1422 
1423  for (auto sn : elem->side_index_range())
1424  {
1425  mesh.get_boundary_info().boundary_ids(elem, sn, bc_ids);
1426  for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
1427  {
1428  const boundary_id_type b_id = *id_it;
1429 
1430  if (mesh_is_serial && b_id == BoundaryInfo::invalid_id)
1431  continue;
1432 
1433  // Make a sorted list of node ids for elem->side(sn)
1434  UniquePtr<Elem> elem_side = elem->build_side_ptr(sn);
1435  std::vector<dof_id_type> elem_side_nodes(elem_side->n_nodes());
1436  for (std::size_t esn=0; esn<elem_side_nodes.size(); ++esn)
1437  elem_side_nodes[esn] = elem_side->node_id(esn);
1438  std::sort(elem_side_nodes.begin(), elem_side_nodes.end());
1439 
1440  for (unsigned int i=0; i != max_subelems; ++i)
1441  if (subelem[i])
1442  {
1443  for (auto subside : subelem[i]->side_index_range())
1444  {
1445  UniquePtr<Elem> subside_elem = subelem[i]->build_side_ptr(subside);
1446 
1447  // Make a list of *vertex* node ids for this subside, see if they are all present
1448  // in elem->side(sn). Note 1: we can't just compare elem->key(sn) to
1449  // subelem[i]->key(subside) in the Prism cases, since the new side is
1450  // a different type. Note 2: we only use vertex nodes since, in the future,
1451  // a Hex20 or Prism15's QUAD8 face may be split into two Tri6 faces, and the
1452  // original face will not contain the mid-edge node.
1453  std::vector<dof_id_type> subside_nodes(subside_elem->n_vertices());
1454  for (std::size_t ssn=0; ssn<subside_nodes.size(); ++ssn)
1455  subside_nodes[ssn] = subside_elem->node_id(ssn);
1456  std::sort(subside_nodes.begin(), subside_nodes.end());
1457 
1458  // std::includes returns true if every element of the second sorted range is
1459  // contained in the first sorted range.
1460  if (std::includes(elem_side_nodes.begin(), elem_side_nodes.end(),
1461  subside_nodes.begin(), subside_nodes.end()))
1462  {
1463  if (b_id != BoundaryInfo::invalid_id)
1464  {
1465  new_bndry_ids.push_back(b_id);
1466  new_bndry_elements.push_back(subelem[i]);
1467  new_bndry_sides.push_back(subside);
1468  }
1469 
1470  // If the original element had a RemoteElem neighbor on side 'sn',
1471  // then the subelem has one on side 'subside'.
1472  if (elem->neighbor_ptr(sn) == remote_elem)
1473  subelem[i]->set_neighbor(subside, const_cast<RemoteElem*>(remote_elem));
1474  }
1475  }
1476  }
1477  } // end for loop over boundary IDs
1478  } // end for loop over sides
1479 
1480  // Remove the original element from the BoundaryInfo structure.
1481  mesh.get_boundary_info().remove(elem);
1482 
1483  } // end if (mesh_has_boundary_data)
1484 
1485  // Determine new IDs for the split elements which will be
1486  // the same on all processors, therefore keeping the Mesh
1487  // in sync. Note: we offset the new IDs by max_orig_id to
1488  // avoid overwriting any of the original IDs.
1489  for (unsigned int i=0; i != max_subelems; ++i)
1490  if (subelem[i])
1491  {
1492  // Determine new IDs for the split elements which will be
1493  // the same on all processors, therefore keeping the Mesh
1494  // in sync. Note: we offset the new IDs by the max of the
1495  // pre-existing ids to avoid conflicting with originals.
1496  subelem[i]->set_id( max_orig_id + 6*elem->id() + i );
1497 
1498 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1499  subelem[i]->set_unique_id() = max_unique_id + 2*elem->unique_id() + i;
1500 #endif
1501 
1502  // Prepare to add the newly-created simplices
1503  new_elements.push_back(subelem[i]);
1504  }
1505 
1506  // Delete the original element
1507  mesh.delete_elem(elem);
1508  } // End for loop over elements
1509  } // end scope
1510 
1511 
1512  // Now, iterate over the new elements vector, and add them each to
1513  // the Mesh.
1514  {
1515  std::vector<Elem *>::iterator el = new_elements.begin();
1516  const std::vector<Elem *>::iterator end = new_elements.end();
1517  for (; el != end; ++el)
1518  mesh.add_elem(*el);
1519  }
1520 
1521  if (mesh_has_boundary_data)
1522  {
1523  // If the old mesh had boundary data, the new mesh better have
1524  // some. However, we can't assert that the size of
1525  // new_bndry_elements vector is > 0, since we may not have split
1526  // any elements actually on the boundary. We also can't assert
1527  // that the original number of boundary sides is equal to the
1528  // sum of the boundary sides currently in the mesh and the
1529  // newly-added boundary sides, since in 3D, we may have split a
1530  // boundary QUAD into two boundary TRIs. Therefore, we won't be
1531  // too picky about the actual number of BCs, and just assert that
1532  // there are some, somewhere.
1533 #ifndef NDEBUG
1534  bool nbe_nonempty = new_bndry_elements.size();
1535  mesh.comm().max(nbe_nonempty);
1536  libmesh_assert(nbe_nonempty ||
1537  mesh.get_boundary_info().n_boundary_conds()>0);
1538 #endif
1539 
1540  // We should also be sure that the lengths of the new boundary data vectors
1541  // are all the same.
1542  libmesh_assert_equal_to (new_bndry_elements.size(), new_bndry_sides.size());
1543  libmesh_assert_equal_to (new_bndry_sides.size(), new_bndry_ids.size());
1544 
1545  // Add the new boundary info to the mesh
1546  for (std::size_t s=0; s<new_bndry_elements.size(); ++s)
1547  mesh.get_boundary_info().add_side(new_bndry_elements[s],
1548  new_bndry_sides[s],
1549  new_bndry_ids[s]);
1550  }
1551 
1552  // In a DistributedMesh any newly added ghost node ids may be
1553  // inconsistent, and unique_ids of newly added ghost nodes remain
1554  // unset.
1555  // make_nodes_parallel_consistent() will fix all this.
1556  if (!mesh.is_serial())
1557  {
1558  mesh.comm().max(added_new_ghost_point);
1559 
1560  if (added_new_ghost_point)
1562  }
1563 
1564 
1565 
1566  // Prepare the newly created mesh for use.
1567  mesh.prepare_for_use(/*skip_renumber =*/ false);
1568 
1569  // Let the new_elements and new_bndry_elements vectors go out of scope.
1570 }
1571 
1572 
1574  const unsigned int n_iterations,
1575  const Real power)
1576 {
1580  libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1581 
1582  /*
1583  * find the boundary nodes
1584  */
1585  std::vector<bool> on_boundary;
1586  MeshTools::find_boundary_nodes(mesh, on_boundary);
1587 
1588  for (unsigned int iter=0; iter<n_iterations; iter++)
1589 
1590  {
1591  /*
1592  * loop over the mesh refinement level
1593  */
1594  unsigned int n_levels = MeshTools::n_levels(mesh);
1595  for (unsigned int refinement_level=0; refinement_level != n_levels;
1596  refinement_level++)
1597  {
1598  // initialize the storage (have to do it on every level to get empty vectors
1599  std::vector<Point> new_positions;
1600  std::vector<Real> weight;
1601  new_positions.resize(mesh.n_nodes());
1602  weight.resize(mesh.n_nodes());
1603 
1604  {
1605  /*
1606  * Loop over the elements to calculate new node positions
1607  */
1608  MeshBase::element_iterator el = mesh.level_elements_begin(refinement_level);
1609  const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level);
1610 
1611  for (; el != end; ++el)
1612  {
1613  /*
1614  * Constant handle for the element
1615  */
1616  const Elem * elem = *el;
1617 
1618  /*
1619  * We relax all nodes on level 0 first
1620  * If the element is refined (level > 0), we interpolate the
1621  * parents nodes with help of the embedding matrix
1622  */
1623  if (refinement_level == 0)
1624  {
1625  for (auto s : elem->side_index_range())
1626  {
1627  /*
1628  * Only operate on sides which are on the
1629  * boundary or for which the current element's
1630  * id is greater than its neighbor's.
1631  * Sides get only built once.
1632  */
1633  if ((elem->neighbor_ptr(s) != libmesh_nullptr) &&
1634  (elem->id() > elem->neighbor_ptr(s)->id()))
1635  {
1637 
1638  const Node & node0 = side->node_ref(0);
1639  const Node & node1 = side->node_ref(1);
1640 
1641  Real node_weight = 1.;
1642  // calculate the weight of the nodes
1643  if (power > 0)
1644  {
1645  Point diff = node0-node1;
1646  node_weight = std::pow(diff.norm(), power);
1647  }
1648 
1649  const dof_id_type id0 = node0.id(), id1 = node1.id();
1650  new_positions[id0].add_scaled( node1, node_weight );
1651  new_positions[id1].add_scaled( node0, node_weight );
1652  weight[id0] += node_weight;
1653  weight[id1] += node_weight;
1654  }
1655  } // element neighbor loop
1656  }
1657 #ifdef LIBMESH_ENABLE_AMR
1658  else // refinement_level > 0
1659  {
1660  /*
1661  * Find the positions of the hanging nodes of refined elements.
1662  * We do this by calculating their position based on the parent
1663  * (one level less refined) element, and the embedding matrix
1664  */
1665 
1666  const Elem * parent = elem->parent();
1667 
1668  /*
1669  * find out which child I am
1670  */
1671  unsigned int c = parent->which_child_am_i(elem);
1672  /*
1673  *loop over the childs (that is, the current elements) nodes
1674  */
1675  for (auto nc : elem->node_index_range())
1676  {
1677  /*
1678  * the new position of the node
1679  */
1680  Point point;
1681  for (auto n : parent->node_index_range())
1682  {
1683  /*
1684  * The value from the embedding matrix
1685  */
1686  const float em_val = parent->embedding_matrix(c,nc,n);
1687 
1688  if (em_val != 0.)
1689  point.add_scaled (parent->point(n), em_val);
1690  }
1691 
1692  const dof_id_type id = elem->node_ptr(nc)->id();
1693  new_positions[id] = point;
1694  weight[id] = 1.;
1695  }
1696  } // if element refinement_level
1697 #endif // #ifdef LIBMESH_ENABLE_AMR
1698 
1699  } // element loop
1700 
1701  /*
1702  * finally reposition the vertex nodes
1703  */
1704  for (unsigned int nid=0; nid<mesh.n_nodes(); ++nid)
1705  if (!on_boundary[nid] && weight[nid] > 0.)
1706  mesh.node_ref(nid) = new_positions[nid]/weight[nid];
1707  }
1708 
1709  {
1710  /*
1711  * Now handle the additional second_order nodes by calculating
1712  * their position based on the vertex positions
1713  * we do a second loop over the level elements
1714  */
1715  MeshBase::element_iterator el = mesh.level_elements_begin(refinement_level);
1716  const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level);
1717 
1718  for (; el != end; ++el)
1719  {
1720  /*
1721  * Constant handle for the element
1722  */
1723  const Elem * elem = *el;
1724  const unsigned int son_begin = elem->n_vertices();
1725  const unsigned int son_end = elem->n_nodes();
1726  for (unsigned int n=son_begin; n<son_end; n++)
1727  {
1728  const unsigned int n_adjacent_vertices =
1730 
1731  Point point;
1732  for (unsigned int v=0; v<n_adjacent_vertices; v++)
1733  point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
1734 
1735  const dof_id_type id = elem->node_ptr(n)->id();
1736  mesh.node_ref(id) = point/n_adjacent_vertices;
1737  }
1738  }
1739  }
1740 
1741  } // refinement_level loop
1742 
1743  } // end iteration
1744 }
1745 
1746 
1747 
1748 #ifdef LIBMESH_ENABLE_AMR
1750 {
1751  // Algorithm:
1752  // .) For each active element in the mesh: construct a
1753  // copy which is the same in every way *except* it is
1754  // a level 0 element. Store the pointers to these in
1755  // a separate vector. Save any boundary information as well.
1756  // Delete the active element from the mesh.
1757  // .) Loop over all (remaining) elements in the mesh, delete them.
1758  // .) Add the level-0 copies back to the mesh
1759 
1760  // Temporary storage for new element pointers
1761  std::vector<Elem *> new_elements;
1762 
1763  // BoundaryInfo Storage for element ids, sides, and BC ids
1764  std::vector<Elem *> saved_boundary_elements;
1765  std::vector<boundary_id_type> saved_bc_ids;
1766  std::vector<unsigned short int> saved_bc_sides;
1767 
1768  // Container to catch boundary ids passed back by BoundaryInfo
1769  std::vector<boundary_id_type> bc_ids;
1770 
1771  // Reserve a reasonable amt. of space for each
1772  new_elements.reserve(mesh.n_active_elem());
1773  saved_boundary_elements.reserve(mesh.get_boundary_info().n_boundary_conds());
1774  saved_bc_ids.reserve(mesh.get_boundary_info().n_boundary_conds());
1775  saved_bc_sides.reserve(mesh.get_boundary_info().n_boundary_conds());
1776 
1777  for (auto & elem : mesh.active_element_ptr_range())
1778  {
1779  // Make a new element of the same type
1780  Elem * copy = Elem::build(elem->type()).release();
1781 
1782  // Set node pointers (they still point to nodes in the original mesh)
1783  for (auto n : elem->node_index_range())
1784  copy->set_node(n) = elem->node_ptr(n);
1785 
1786  // Copy over ids
1787  copy->processor_id() = elem->processor_id();
1788  copy->subdomain_id() = elem->subdomain_id();
1789 
1790  // Retain the original element's ID(s) as well, otherwise
1791  // the Mesh may try to create them for you...
1792  copy->set_id( elem->id() );
1793 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1794  copy->set_unique_id() = elem->unique_id();
1795 #endif
1796 
1797  // This element could have boundary info or DistributedMesh
1798  // remote_elem links as well. We need to save the (elem,
1799  // side, bc_id) triples and those links
1800  for (auto s : elem->side_index_range())
1801  {
1802  if (elem->neighbor_ptr(s) == remote_elem)
1803  copy->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
1804 
1805  mesh.get_boundary_info().boundary_ids(elem, s, bc_ids);
1806  for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
1807  {
1808  const boundary_id_type bc_id = *id_it;
1809 
1810  if (bc_id != BoundaryInfo::invalid_id)
1811  {
1812  saved_boundary_elements.push_back(copy);
1813  saved_bc_ids.push_back(bc_id);
1814  saved_bc_sides.push_back(s);
1815  }
1816  }
1817  }
1818 
1819 
1820  // We're done with this element
1821  mesh.delete_elem(elem);
1822 
1823  // But save the copy
1824  new_elements.push_back(copy);
1825  }
1826 
1827  // Make sure we saved the same number of boundary conditions
1828  // in each vector.
1829  libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
1830  libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
1831 
1832  // Loop again, delete any remaining elements
1833  for (auto & elem : mesh.element_ptr_range())
1834  mesh.delete_elem(elem);
1835 
1836  // Add the copied (now level-0) elements back to the mesh
1837  {
1838  for (std::vector<Elem *>::iterator it = new_elements.begin();
1839  it != new_elements.end();
1840  ++it)
1841  {
1842 #ifndef NDEBUG
1843  dof_id_type orig_id = (*it)->id();
1844 
1845  // ugly mid-statement endif to avoid unused variable warnings
1846  Elem * added_elem =
1847 #endif
1848  mesh.add_elem(*it);
1849 
1850 #ifndef NDEBUG
1851  dof_id_type added_id = added_elem->id();
1852 #endif
1853 
1854  // If the Elem, as it was re-added to the mesh, now has a
1855  // different ID (this is unlikely, so it's just an assert)
1856  // the boundary information will no longer be correct.
1857  libmesh_assert_equal_to (orig_id, added_id);
1858  }
1859  }
1860 
1861  // Finally, also add back the saved boundary information
1862  for (std::size_t e=0; e<saved_boundary_elements.size(); ++e)
1863  mesh.get_boundary_info().add_side(saved_boundary_elements[e],
1864  saved_bc_sides[e],
1865  saved_bc_ids[e]);
1866 
1867  // Trim unused and renumber nodes and elements
1868  mesh.prepare_for_use(/*skip_renumber =*/ false);
1869 }
1870 #endif // #ifdef LIBMESH_ENABLE_AMR
1871 
1872 
1873 
1875  const boundary_id_type old_id,
1876  const boundary_id_type new_id)
1877 {
1878  if (old_id == new_id)
1879  {
1880  // If the IDs are the same, this is a no-op.
1881  return;
1882  }
1883 
1884  // A reference to the Mesh's BoundaryInfo object, for convenience.
1885  BoundaryInfo & bi = mesh.get_boundary_info();
1886 
1887  {
1888  // Build a list of all nodes that have boundary IDs
1889  std::vector<dof_id_type> node_list;
1890  std::vector<boundary_id_type> bc_id_list;
1891  bi.build_node_list (node_list, bc_id_list);
1892 
1893  // Temporary vector to hold ids
1894  std::vector<boundary_id_type> bndry_ids;
1895 
1896  // For each node with the old_id...
1897  for (std::size_t idx=0; idx<node_list.size(); ++idx)
1898  if (bc_id_list[idx] == old_id)
1899  {
1900  // Get the node in question
1901  const Node * node = mesh.node_ptr(node_list[idx]);
1902 
1903  // Get all the current IDs for this node.
1904  bi.boundary_ids(node, bndry_ids);
1905 
1906  // Update the IDs accordingly
1907  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
1908 
1909  // Remove all traces of that node from the BoundaryInfo object
1910  bi.remove(node);
1911 
1912  // Add it back with the updated IDs
1913  bi.add_node(node, bndry_ids);
1914  }
1915  }
1916 
1917  {
1918  // Build a list of all edges that have boundary IDs
1919  std::vector<dof_id_type> elem_list;
1920  std::vector<unsigned short int> edge_list;
1921  std::vector<boundary_id_type> bc_id_list;
1922  bi.build_edge_list (elem_list, edge_list, bc_id_list);
1923 
1924  // Temporary vector to hold ids
1925  std::vector<boundary_id_type> bndry_ids;
1926 
1927  // For each edge with the old_id...
1928  for (std::size_t idx=0; idx<elem_list.size(); ++idx)
1929  if (bc_id_list[idx] == old_id)
1930  {
1931  // Get the elem in question
1932  const Elem * elem = mesh.elem_ptr(elem_list[idx]);
1933 
1934  // The edge of the elem in question
1935  unsigned short int edge = edge_list[idx];
1936 
1937  // Get all the current IDs for the edge in question.
1938  bi.edge_boundary_ids(elem, edge, bndry_ids);
1939 
1940  // Update the IDs accordingly
1941  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
1942 
1943  // Remove all traces of that edge from the BoundaryInfo object
1944  bi.remove_edge(elem, edge);
1945 
1946  // Add it back with the updated IDs
1947  bi.add_edge(elem, edge, bndry_ids);
1948  }
1949  }
1950 
1951  {
1952  // Build a list of all shell-faces that have boundary IDs
1953  std::vector<dof_id_type> elem_list;
1954  std::vector<unsigned short int> shellface_list;
1955  std::vector<boundary_id_type> bc_id_list;
1956  bi.build_shellface_list (elem_list, shellface_list, bc_id_list);
1957 
1958  // Temporary vector to hold ids
1959  std::vector<boundary_id_type> bndry_ids;
1960 
1961  // For each shellface with the old_id...
1962  for (std::size_t idx=0; idx<elem_list.size(); ++idx)
1963  if (bc_id_list[idx] == old_id)
1964  {
1965  // Get the elem in question
1966  const Elem * elem = mesh.elem_ptr(elem_list[idx]);
1967 
1968  // The shellface of the elem in question
1969  unsigned short int shellface = shellface_list[idx];
1970 
1971  // Get all the current IDs for the shellface in question.
1972  bi.shellface_boundary_ids(elem, shellface, bndry_ids);
1973 
1974  // Update the IDs accordingly
1975  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
1976 
1977  // Remove all traces of that shellface from the BoundaryInfo object
1978  bi.remove_shellface(elem, shellface);
1979 
1980  // Add it back with the updated IDs
1981  bi.add_shellface(elem, shellface, bndry_ids);
1982  }
1983  }
1984 
1985  {
1986  // Build a list of all sides that have boundary IDs
1987  std::vector<dof_id_type> elem_list;
1988  std::vector<unsigned short int> side_list;
1989  std::vector<boundary_id_type> bc_id_list;
1990  bi.build_side_list (elem_list, side_list, bc_id_list);
1991 
1992  // Temporary vector to hold ids
1993  std::vector<boundary_id_type> bndry_ids;
1994 
1995  // For each side with the old_id...
1996  for (std::size_t idx=0; idx<elem_list.size(); ++idx)
1997  if (bc_id_list[idx] == old_id)
1998  {
1999  // Get the elem in question
2000  const Elem * elem = mesh.elem_ptr(elem_list[idx]);
2001 
2002  // The side of the elem in question
2003  unsigned short int side = side_list[idx];
2004 
2005  // Get all the current IDs for the side in question.
2006  bi.boundary_ids(elem, side, bndry_ids);
2007 
2008  // Update the IDs accordingly
2009  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
2010 
2011  // Remove all traces of that side from the BoundaryInfo object
2012  bi.remove_side(elem, side);
2013 
2014  // Add it back with the updated IDs
2015  bi.add_side(elem, side, bndry_ids);
2016  }
2017  }
2018 
2019  // Remove any remaining references to the old_id so it does not show
2020  // up in lists of boundary ids, etc.
2021  bi.remove_id(old_id);
2022 }
2023 
2024 
2025 
2027  const subdomain_id_type old_id,
2028  const subdomain_id_type new_id)
2029 {
2030  if (old_id == new_id)
2031  {
2032  // If the IDs are the same, this is a no-op.
2033  return;
2034  }
2035 
2036  for (auto & elem : mesh.element_ptr_range())
2037  {
2038  if (elem->subdomain_id() == old_id)
2039  elem->subdomain_id() = new_id;
2040  }
2041 }
2042 
2043 
2044 } // namespace libMesh
The definition of the element_iterator struct.
Definition: mesh_base.h:1476
OStreamProxy err
void set_p_level(const unsigned int p)
Sets the value of the p-refinement level for the element.
Definition: elem.h:2559
virtual unique_id_type parallel_max_unique_id() const =0
bool has_children() const
Definition: elem.h:2295
The Tri3 is an element in 2D composed of 3 nodes.
Definition: face_tri3.h:54
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
unique_id_type & set_unique_id()
Definition: dof_object.h:662
virtual void reserve_nodes(const dof_id_type nn)=0
Reserves space for a known number of nodes.
virtual UniquePtr< Elem > build_side_ptr(const unsigned int i, bool proxy=true)=0
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int n) const
Definition: elem.C:2698
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
void remove_edge(const Elem *elem, const unsigned short int edge)
Removes all boundary conditions associated with edge edge of element elem, if any exist...
A Node is like a Point, but with more information.
Definition: node.h:52
virtual dof_id_type n_active_elem() const =0
void rotate(MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
Rotates the mesh in the xy plane.
virtual bool is_serial() const
Definition: mesh_base.h:140
virtual const Point & point(const dof_id_type i) const =0
static UniquePtr< Elem > build(const ElemType type, Elem *p=libmesh_nullptr)
Definition: elem.C:238
void set_parent(Elem *p)
Sets the pointer to the element&#39;s parent.
Definition: elem.h:2362
virtual void all_first_order() libmesh_override
Converts a mesh with higher-order elements into a mesh with linear elements.
virtual ElemType type() const =0
virtual element_iterator level_elements_begin(unsigned int level)=0
Iterate over elements of a given level.
unsigned int p_level() const
Definition: elem.h:2422
void max(T &r) const
Take a local variable and replace it with the maximum of it&#39;s values on all processors.
void add_scaled(const TypeVector< T2 > &, const T)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:624
void remove_id(boundary_id_type id)
Removes all entities (nodes, sides, edges, shellfaces) with boundary id id from their respective cont...
void scale(MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
Scales the mesh.
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2083
void remove_shellface(const Elem *elem, const unsigned short int shellface)
Removes all boundary conditions associated with shell face shellface of element elem, if any exist.
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:416
void distort(MeshBase &mesh, const Real factor, const bool perturb_boundary=false)
Randomly perturb the nodal locations.
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:420
ElemType
Defines an enum for geometric element types.
IntRange< unsigned short > node_index_range() const
Definition: elem.h:2065
const Elem * parent() const
Definition: elem.h:2346
virtual void all_second_order(const bool full_ordered=true) libmesh_override
Converts a (conforming, non-refined) mesh with linear elements into a mesh with second-order elements...
void remove(const Node *node)
Removes the boundary conditions associated with node node, if any exist.
unsigned short int side
Definition: xdr_io.C:49
virtual dof_id_type max_node_id() const =0
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
virtual Real hmin() const
Definition: elem.C:458
Real norm() const
Definition: type_vector.h:909
MeshBase & mesh
void add_child(Elem *elem)
Adds a child pointer to the array of children of this element.
Definition: elem.C:1608
const class libmesh_nullptr_t libmesh_nullptr
void set_refinement_flag(const RefinementState rflag)
Sets the value of the refinement flag for the element.
Definition: elem.h:2513
virtual const Node * node_ptr(const dof_id_type i) const =0
std::size_t n_boundary_conds() const
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
The libMesh namespace provides an interface to certain functionality in the library.
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.
long double max(long double a, double b)
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.
boundary_id_type bc_id
Definition: xdr_io.C:50
This is the MeshBase class.
Definition: mesh_base.h:68
void add(const TypeVector< T2 > &)
Add to this vector without creating a temporary.
Definition: type_vector.h:600
dof_id_type & set_id()
Definition: dof_object.h:641
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
The Tri6 is an element in 2D composed of 6 nodes.
Definition: face_tri6.h:49
virtual unsigned int n_nodes() const =0
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:1967
std::vector< boundary_id_type > boundary_ids(const Node *node) const
void change_subdomain_id(MeshBase &mesh, const subdomain_id_type old_id, const subdomain_id_type new_id)
Finds any subdomain ids that are currently old_id, changes them to new_id.
virtual node_iterator nodes_begin()=0
Iterate over all the nodes in the Mesh.
virtual element_iterator elements_begin()=0
Iterate over all the elements in the Mesh.
virtual element_iterator level_elements_end(unsigned int level)=0
void all_tri(MeshBase &mesh)
Converts the 2D quadrilateral elements of a Mesh into triangular elements.
void replace_child(Elem *elem, unsigned int c)
Replaces the child pointer at the specified index in the child array.
Definition: elem.C:1654
virtual dof_id_type max_elem_id() const =0
void smooth(MeshBase &, unsigned int, Real)
Smooth the mesh with a simple Laplace smoothing algorithm.
void add_node(const Node *node, const boundary_id_type id)
Add Node node with boundary id id to the boundary information data structures.
PetscErrorCode Vec x
int8_t boundary_id_type
Definition: id_types.h:51
void shellface_boundary_ids(const Elem *const elem, const unsigned short int shellface, std::vector< boundary_id_type > &vec_to_fill) const
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.
virtual SimpleRange< element_iterator > element_ptr_range()=0
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1874
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:245
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
virtual element_iterator elements_end()=0
void change_boundary_id(MeshBase &mesh, const boundary_id_type old_id, const boundary_id_type new_id)
Finds any boundary ids that are currently old_id, changes them to new_id.
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:603
virtual SimpleRange< node_iterator > node_ptr_range()=0
virtual const Node * query_node_ptr(const dof_id_type i) const =0
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:56
const Elem * child_ptr(unsigned int i) const
Definition: elem.h:2445
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:324
virtual void delete_node(Node *n)=0
Removes the Node n from the mesh.
RefinementState p_refinement_flag() const
Definition: elem.h:2521
virtual Order default_order() const =0
void prepare_for_use(const bool skip_renumber_nodes_and_elements=false, const bool skip_find_neighbors=false)
Prepare a newly created (or read) mesh for use.
Definition: mesh_base.C:174
subdomain_id_type subdomain_id() const
Definition: elem.h:1951
virtual unsigned int n_sides() const =0
void set_neighbor(const unsigned int i, Elem *n)
Assigns n as the neighbor.
Definition: elem.h:2000
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 regenerate_id_sets()
Clears and regenerates the cached sets of ids.
void translate(MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.)
Translates the mesh.
virtual unsigned int n_children() const =0
static ElemType second_order_equivalent_type(const ElemType et, const bool full_ordered=true)
Definition: elem.C:2783
double pow(double a, int b)
TypeVector< T > unit() const
Definition: type_vector.C:37
SimpleRange< NodeRefIter > node_ref_range()
Returns a range with all nodes of an element, usable in range-based for loops.
Definition: elem.h:2047
The definition of the node_iterator struct.
Definition: mesh_base.h:1528
unsigned int which_child_am_i(const Elem *e) const
Definition: elem.h:2487
virtual Elem * insert_elem(Elem *e)=0
Insert elem e to the element array, preserving its id and replacing/deleting any existing element wit...
void remove_side(const Elem *elem, const unsigned short int side)
Removes all boundary conditions associated with side side of element elem, if any exist...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual node_iterator nodes_end()=0
const Point & point(const unsigned int i) const
Definition: elem.h:1809
void copy_boundary_ids(const BoundaryInfo &old_boundary_info, const Elem *const old_elem, const Elem *const new_elem)
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 make_nodes_parallel_consistent(MeshBase &)
Copy processor_ids and ids on ghost nodes from their local processors.
const Parallel::Communicator & comm() const
std::vector< boundary_id_type > edge_boundary_ids(const Elem *const elem, const unsigned short int edge) const
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
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...
bool _is_prepared
Flag indicating if the mesh has been prepared for use.
Definition: mesh_base.h:1347
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 ...
virtual unsigned int n_vertices() const =0
unsigned int level() const
Definition: elem.h:2388
RefinementState refinement_flag() const
Definition: elem.h:2505
virtual UniquePtr< FunctionBase< Output > > clone() const =0
virtual const Node & node(const dof_id_type i) const
Definition: mesh_base.h:440
void redistribute(MeshBase &mesh, const FunctionBase< Real > &mapfunc)
Deterministically perturb the nodal locations.
Defines a dense vector for use in Finite Element-type computations.
void set_p_refinement_flag(const RefinementState pflag)
Sets the value of the p-refinement flag for the element.
Definition: elem.h:2529
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1831
dof_id_type id() const
Definition: dof_object.h:632
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:649
static ElemType first_order_equivalent_type(const ElemType et)
Definition: elem.C:2724
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const
Definition: elem.C:2706
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.
long double min(long double a, double b)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
uint8_t unique_id_type
Definition: id_types.h:79
void flatten(MeshBase &mesh)
Removes all the refinement tree structure of Mesh, leaving only the highest-level (most-refined) elem...
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...
void find_boundary_nodes(const MeshBase &mesh, std::vector< bool > &on_boundary)
Calling this function on a 2D mesh will convert all the elements to triangles.
Definition: mesh_tools.C:290
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.
processor_id_type processor_id() const
uint8_t dof_id_type
Definition: id_types.h:64
processor_id_type processor_id() const
Definition: dof_object.h:694
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 Real pi
.
Definition: libmesh.h:172
virtual float embedding_matrix(const unsigned int child_num, const unsigned int child_node_num, const unsigned int parent_node_num) const =0
const RemoteElem * remote_elem
Definition: remote_elem.C:57