libMesh
mesh_modification.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // 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 #include <array>
27 
28 // Local includes
29 #include "libmesh/boundary_info.h"
30 #include "libmesh/function_base.h"
31 #include "libmesh/cell_tet4.h"
32 #include "libmesh/cell_tet10.h"
33 #include "libmesh/face_tri3.h"
34 #include "libmesh/face_tri6.h"
35 #include "libmesh/libmesh_logging.h"
36 #include "libmesh/mesh_communication.h"
37 #include "libmesh/mesh_modification.h"
38 #include "libmesh/mesh_tools.h"
39 #include "libmesh/parallel.h"
40 #include "libmesh/remote_elem.h"
41 #include "libmesh/enum_to_string.h"
42 #include "libmesh/unstructured_mesh.h"
43 #include "libmesh/elem_side_builder.h"
44 #include "libmesh/tensor_value.h"
45 
46 namespace
47 {
48 using namespace libMesh;
49 
50 bool split_first_diagonal(const Elem * elem,
51  unsigned int diag_1_node_1,
52  unsigned int diag_1_node_2,
53  unsigned int diag_2_node_1,
54  unsigned int diag_2_node_2)
55 {
56  return ((elem->node_id(diag_1_node_1) > elem->node_id(diag_2_node_1) &&
57  elem->node_id(diag_1_node_1) > elem->node_id(diag_2_node_2)) ||
58  (elem->node_id(diag_1_node_2) > elem->node_id(diag_2_node_1) &&
59  elem->node_id(diag_1_node_2) > elem->node_id(diag_2_node_2)));
60 }
61 
62 
63 // Return the local index of the vertex on \p elem with the highest
64 // node id.
65 unsigned int highest_vertex_on(const Elem * elem)
66 {
67  unsigned int highest_n = 0;
68  dof_id_type highest_n_id = elem->node_id(0);
69  for (auto n : make_range(1u, elem->n_vertices()))
70  {
71  const dof_id_type n_id = elem->node_id(n);
72  if (n_id > highest_n_id)
73  {
74  highest_n = n;
75  highest_n_id = n_id;
76  }
77  }
78 
79  return highest_n;
80 }
81 
82 
83 static const std::array<std::array<unsigned int, 3>, 8> opposing_nodes =
84 {{ {2,5,7},{3,4,6},{0,5,7},{1,4,6},{1,3,6},{0,2,7},{1,3,4},{0,2,5} }};
85 
86 
87 // Find the highest id on these side nodes of this element
88 std::pair<unsigned int, unsigned int>
89 split_diagonal(const Elem * elem,
90  const std::vector<unsigned int> & nodes_on_side)
91 {
92  libmesh_assert_equal_to(elem->type(), HEX8);
93 
94  unsigned int highest_n = nodes_on_side.front();
95  dof_id_type highest_n_id = elem->node_id(nodes_on_side.front());
96  for (auto n : nodes_on_side)
97  {
98  const dof_id_type n_id = elem->node_id(n);
99  if (n_id > highest_n_id)
100  {
101  highest_n = n;
102  highest_n_id = n_id;
103  }
104  }
105 
106  for (auto n : nodes_on_side)
107  {
108  for (auto n2 : opposing_nodes[highest_n])
109  if (n2 == n)
110  return std::make_pair(highest_n, n2);
111  }
112 
113  libmesh_error();
114 
115  return std::make_pair(libMesh::invalid_uint, libMesh::invalid_uint);
116 }
117 
118 
119 // Reconstruct a C++20 feature in C++14
120 template <typename T>
121 struct reversion_wrapper { T& iterable; };
122 
123 template <typename T>
124 auto begin (reversion_wrapper<T> w) {return std::rbegin(w.iterable);}
125 
126 template <typename T>
127 auto end (reversion_wrapper<T> w) {return std::rend(w.iterable);}
128 
129 template <typename T>
130 reversion_wrapper<T> reverse(T&& iterable) {return {iterable};}
131 
132 }
133 
134 
135 namespace libMesh
136 {
137 
138 
139 // ------------------------------------------------------------
140 // MeshTools::Modification functions for mesh modification
142  const Real factor,
143  const bool perturb_boundary)
144 {
147  libmesh_assert ((factor >= 0.) && (factor <= 1.));
148 
149  LOG_SCOPE("distort()", "MeshTools::Modification");
150 
151  // If we are not perturbing boundary nodes, make a
152  // quickly-searchable list of node ids we can check against.
153  std::unordered_set<dof_id_type> boundary_node_ids;
154  if (!perturb_boundary)
155  boundary_node_ids = MeshTools::find_boundary_nodes (mesh);
156 
157  // Now calculate the minimum distance to
158  // neighboring nodes for each node.
159  // hmin holds these distances.
160  std::vector<float> hmin (mesh.max_node_id(),
161  std::numeric_limits<float>::max());
162 
163  for (const auto & elem : mesh.active_element_ptr_range())
164  for (auto & n : elem->node_ref_range())
165  hmin[n.id()] = std::min(hmin[n.id()],
166  static_cast<float>(elem->hmin()));
167 
168  // Now actually move the nodes
169  {
170  const unsigned int seed = 123456;
171 
172  // seed the random number generator.
173  // We'll loop from 1 to n_nodes on every processor, even those
174  // that don't have a particular node, so that the pseudorandom
175  // numbers will be the same everywhere.
176  std::srand(seed);
177 
178  // If the node is on the boundary or
179  // the node is not used by any element (hmin[n]<1.e20)
180  // then we should not move it.
181  // [Note: Testing for (in)equality might be wrong
182  // (different types, namely float and double)]
183  for (auto n : make_range(mesh.max_node_id()))
184  if ((perturb_boundary || !boundary_node_ids.count(n)) && hmin[n] < 1.e20)
185  {
186  // the direction, random but unit normalized
187  Point dir (static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
188  (mesh.mesh_dimension() > 1) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.,
189  ((mesh.mesh_dimension() == 3) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.));
190 
191  dir(0) = (dir(0)-.5)*2.;
192 #if LIBMESH_DIM > 1
193  if (mesh.mesh_dimension() > 1)
194  dir(1) = (dir(1)-.5)*2.;
195 #endif
196 #if LIBMESH_DIM > 2
197  if (mesh.mesh_dimension() == 3)
198  dir(2) = (dir(2)-.5)*2.;
199 #endif
200 
201  dir = dir.unit();
202 
203  Node * node = mesh.query_node_ptr(n);
204  if (!node)
205  continue;
206 
207  (*node)(0) += dir(0)*factor*hmin[n];
208 #if LIBMESH_DIM > 1
209  if (mesh.mesh_dimension() > 1)
210  (*node)(1) += dir(1)*factor*hmin[n];
211 #endif
212 #if LIBMESH_DIM > 2
213  if (mesh.mesh_dimension() == 3)
214  (*node)(2) += dir(2)*factor*hmin[n];
215 #endif
216  }
217  }
218 }
219 
220 
221 
223 {
224  LOG_SCOPE("permute_elements()", "MeshTools::Modification");
225 
226  // We don't yet support doing permute() on a parent element, which
227  // would require us to consistently permute all its children and
228  // give them different local child numbers.
229  unsigned int n_levels = MeshTools::n_levels(mesh);
230  if (n_levels > 1)
231  libmesh_error();
232 
233  const unsigned int seed = 123456;
234 
235  // seed the random number generator.
236  // We'll loop from 1 to max_elem_id on every processor, even those
237  // that don't have a particular element, so that the pseudorandom
238  // numbers will be the same everywhere.
239  std::srand(seed);
240 
241 
242  for (auto e_id : make_range(mesh.max_elem_id()))
243  {
244  int my_rand = std::rand();
245 
246  Elem * elem = mesh.query_elem_ptr(e_id);
247 
248  if (!elem)
249  continue;
250 
251  const unsigned int max_permutation = elem->n_permutations();
252  if (!max_permutation)
253  continue;
254 
255  const unsigned int perm = my_rand % max_permutation;
256 
257  elem->permute(perm);
258  }
259 }
260 
261 
263 {
264  LOG_SCOPE("orient_elements()", "MeshTools::Modification");
265 
266  // We don't yet support doing orient() on a parent element, which
267  // would require us to consistently orient all its children and
268  // give them different local child numbers.
269  unsigned int n_levels = MeshTools::n_levels(mesh);
270  if (n_levels > 1)
271  libmesh_not_implemented_msg("orient_elements() does not support refined meshes");
272 
273  BoundaryInfo & boundary_info = mesh.get_boundary_info();
274  for (auto elem : mesh.element_ptr_range())
275  elem->orient(&boundary_info);
276 }
277 
278 
279 
281  const FunctionBase<Real> & mapfunc)
282 {
285 
286  LOG_SCOPE("redistribute()", "MeshTools::Modification");
287 
288  DenseVector<Real> output_vec(LIBMESH_DIM);
289 
290  // FIXME - we should thread this later.
291  std::unique_ptr<FunctionBase<Real>> myfunc = mapfunc.clone();
292 
293  for (auto & node : mesh.node_ptr_range())
294  {
295  (*myfunc)(*node, output_vec);
296 
297  (*node)(0) = output_vec(0);
298 #if LIBMESH_DIM > 1
299  (*node)(1) = output_vec(1);
300 #endif
301 #if LIBMESH_DIM > 2
302  (*node)(2) = output_vec(2);
303 #endif
304  }
305 }
306 
307 
308 
310  const Real xt,
311  const Real yt,
312  const Real zt)
313 {
314  const Point p(xt, yt, zt);
315 
316  for (auto & node : mesh.node_ptr_range())
317  *node += p;
318 }
319 
320 
321 // void MeshTools::Modification::rotate2D (MeshBase & mesh,
322 // const Real alpha)
323 // {
324 // libmesh_assert_not_equal_to (mesh.mesh_dimension(), 1);
325 
326 // const Real pi = std::acos(-1);
327 // const Real a = alpha/180.*pi;
328 // for (unsigned int n=0; n<mesh.n_nodes(); n++)
329 // {
330 // const Point p = mesh.node_ref(n);
331 // const Real x = p(0);
332 // const Real y = p(1);
333 // const Real z = p(2);
334 // mesh.node_ref(n) = Point(std::cos(a)*x - std::sin(a)*y,
335 // std::sin(a)*x + std::cos(a)*y,
336 // z);
337 // }
338 
339 // }
340 
341 
342 
345  const Real phi,
346  const Real theta,
347  const Real psi)
348 {
349 #if LIBMESH_DIM == 3
350  const auto R = RealTensorValue::intrinsic_rotation_matrix(phi, theta, psi);
351 
352  if (theta)
354 
355  for (auto & node : mesh.node_ptr_range())
356  {
357  Point & pt = *node;
358  pt = R * pt;
359  }
360 
361  return R;
362 
363 #else
364  libmesh_ignore(mesh, phi, theta, psi);
365  libmesh_error_msg("MeshTools::Modification::rotate() requires libMesh to be compiled with LIBMESH_DIM==3");
366  // We'll never get here
367  return RealTensorValue();
368 #endif
369 }
370 
371 
373  const Real xs,
374  const Real ys,
375  const Real zs)
376 {
377  const Real x_scale = xs;
378  Real y_scale = ys;
379  Real z_scale = zs;
380 
381  if (ys == 0.)
382  {
383  libmesh_assert_equal_to (zs, 0.);
384 
385  y_scale = z_scale = x_scale;
386  }
387 
388  // Scale the x coordinate in all dimensions
389  for (auto & node : mesh.node_ptr_range())
390  (*node)(0) *= x_scale;
391 
392  // Only scale the y coordinate in 2 and 3D
393  if (LIBMESH_DIM < 2)
394  return;
395 
396  for (auto & node : mesh.node_ptr_range())
397  (*node)(1) *= y_scale;
398 
399  // Only scale the z coordinate in 3D
400  if (LIBMESH_DIM < 3)
401  return;
402 
403  for (auto & node : mesh.node_ptr_range())
404  (*node)(2) *= z_scale;
405 }
406 
407 
408 
410 {
411  if (!mesh.is_replicated() && !mesh.is_prepared())
413 
414  // The number of elements in the original mesh before any additions
415  // or deletions.
416  const dof_id_type n_orig_elem = mesh.n_elem();
417  const dof_id_type max_orig_id = mesh.max_elem_id();
418 
419  // We store pointers to the newly created elements in a vector
420  // until they are ready to be added to the mesh. This is because
421  // adding new elements on the fly can cause reallocation and invalidation
422  // of existing mesh element_iterators.
423  std::vector<std::unique_ptr<Elem>> new_elements;
424 
425  unsigned int max_subelems = 1; // in 1D nothing needs to change
426  if (mesh.mesh_dimension() == 2) // in 2D quads can split into 2 tris
427  max_subelems = 2;
428  if (mesh.mesh_dimension() == 3) // in 3D hexes can split into 6 tets
429  max_subelems = 6;
430 
431  new_elements.reserve (max_subelems*n_orig_elem);
432 
433  // If the original mesh has *side* boundary data, we carry that over
434  // to the new mesh with triangular elements. We currently only
435  // support bringing over side-based BCs to the all-tri mesh, but
436  // that could probably be extended to node and edge-based BCs as
437  // well.
438  const bool mesh_has_boundary_data = (mesh.get_boundary_info().n_boundary_conds() > 0);
439 
440  // Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
441  std::vector<Elem *> new_bndry_elements;
442  std::vector<unsigned short int> new_bndry_sides;
443  std::vector<boundary_id_type> new_bndry_ids;
444 
445  // We may need to add new points if we run into a 1.5th order
446  // element; if we do that on a DistributedMesh in a ghost element then
447  // we will need to fix their ids / unique_ids
448  bool added_new_ghost_point = false;
449 
450  // Iterate over the elements, splitting:
451  // QUADs into pairs of conforming triangles
452  // PYRAMIDs into pairs of conforming tets,
453  // PRISMs into triplets of conforming tets, and
454  // HEXs into quintets or sextets of conforming tets.
455  // We split on the shortest diagonal to give us better
456  // triangle quality in 2D, and we split based on node ids
457  // to guarantee consistency in 3D.
458 
459  // FIXME: This algorithm does not work on refined grids!
460  {
461 #ifdef LIBMESH_ENABLE_UNIQUE_ID
462  unique_id_type max_unique_id = mesh.parallel_max_unique_id();
463 #endif
464 
465  // For avoiding extraneous allocation when building side elements
466  std::unique_ptr<const Elem> elem_side, subside_elem;
467 
468  for (auto & elem : mesh.element_ptr_range())
469  {
470  const ElemType etype = elem->type();
471 
472  // all_tri currently only works on coarse meshes
473  if (elem->parent())
474  libmesh_not_implemented_msg("Cannot convert a refined element into simplices\n");
475 
476  // The new elements we will split the quad into.
477  // In 3D we may need as many as 6 tets per hex
478  std::array<std::unique_ptr<Elem>, 6> subelem {};
479 
480  switch (etype)
481  {
482  case QUAD4:
483  {
484  subelem[0] = Elem::build(TRI3);
485  subelem[1] = Elem::build(TRI3);
486 
487  // Check for possible edge swap
488  if ((elem->point(0) - elem->point(2)).norm() <
489  (elem->point(1) - elem->point(3)).norm())
490  {
491  subelem[0]->set_node(0, elem->node_ptr(0));
492  subelem[0]->set_node(1, elem->node_ptr(1));
493  subelem[0]->set_node(2, elem->node_ptr(2));
494 
495  subelem[1]->set_node(0, elem->node_ptr(0));
496  subelem[1]->set_node(1, elem->node_ptr(2));
497  subelem[1]->set_node(2, elem->node_ptr(3));
498  }
499 
500  else
501  {
502  subelem[0]->set_node(0, elem->node_ptr(0));
503  subelem[0]->set_node(1, elem->node_ptr(1));
504  subelem[0]->set_node(2, elem->node_ptr(3));
505 
506  subelem[1]->set_node(0, elem->node_ptr(1));
507  subelem[1]->set_node(1, elem->node_ptr(2));
508  subelem[1]->set_node(2, elem->node_ptr(3));
509  }
510 
511 
512  break;
513  }
514 
515  case QUAD8:
516  {
517  if (elem->processor_id() != mesh.processor_id())
518  added_new_ghost_point = true;
519 
520  subelem[0] = Elem::build(TRI6);
521  subelem[1] = Elem::build(TRI6);
522 
523  // Add a new node at the center (vertex average) of the element.
524  Node * new_node = mesh.add_point((mesh.point(elem->node_id(0)) +
525  mesh.point(elem->node_id(1)) +
526  mesh.point(elem->node_id(2)) +
527  mesh.point(elem->node_id(3)))/4,
529  elem->processor_id());
530 
531  // Check for possible edge swap
532  if ((elem->point(0) - elem->point(2)).norm() <
533  (elem->point(1) - elem->point(3)).norm())
534  {
535  subelem[0]->set_node(0, elem->node_ptr(0));
536  subelem[0]->set_node(1, elem->node_ptr(1));
537  subelem[0]->set_node(2, elem->node_ptr(2));
538  subelem[0]->set_node(3, elem->node_ptr(4));
539  subelem[0]->set_node(4, elem->node_ptr(5));
540  subelem[0]->set_node(5, new_node);
541 
542  subelem[1]->set_node(0, elem->node_ptr(0));
543  subelem[1]->set_node(1, elem->node_ptr(2));
544  subelem[1]->set_node(2, elem->node_ptr(3));
545  subelem[1]->set_node(3, new_node);
546  subelem[1]->set_node(4, elem->node_ptr(6));
547  subelem[1]->set_node(5, elem->node_ptr(7));
548 
549  }
550 
551  else
552  {
553  subelem[0]->set_node(0, elem->node_ptr(3));
554  subelem[0]->set_node(1, elem->node_ptr(0));
555  subelem[0]->set_node(2, elem->node_ptr(1));
556  subelem[0]->set_node(3, elem->node_ptr(7));
557  subelem[0]->set_node(4, elem->node_ptr(4));
558  subelem[0]->set_node(5, new_node);
559 
560  subelem[1]->set_node(0, elem->node_ptr(1));
561  subelem[1]->set_node(1, elem->node_ptr(2));
562  subelem[1]->set_node(2, elem->node_ptr(3));
563  subelem[1]->set_node(3, elem->node_ptr(5));
564  subelem[1]->set_node(4, elem->node_ptr(6));
565  subelem[1]->set_node(5, new_node);
566  }
567 
568  break;
569  }
570 
571  case QUAD9:
572  {
573  subelem[0] = Elem::build(TRI6);
574  subelem[1] = Elem::build(TRI6);
575 
576  // Check for possible edge swap
577  if ((elem->point(0) - elem->point(2)).norm() <
578  (elem->point(1) - elem->point(3)).norm())
579  {
580  subelem[0]->set_node(0, elem->node_ptr(0));
581  subelem[0]->set_node(1, elem->node_ptr(1));
582  subelem[0]->set_node(2, elem->node_ptr(2));
583  subelem[0]->set_node(3, elem->node_ptr(4));
584  subelem[0]->set_node(4, elem->node_ptr(5));
585  subelem[0]->set_node(5, elem->node_ptr(8));
586 
587  subelem[1]->set_node(0, elem->node_ptr(0));
588  subelem[1]->set_node(1, elem->node_ptr(2));
589  subelem[1]->set_node(2, elem->node_ptr(3));
590  subelem[1]->set_node(3, elem->node_ptr(8));
591  subelem[1]->set_node(4, elem->node_ptr(6));
592  subelem[1]->set_node(5, elem->node_ptr(7));
593  }
594 
595  else
596  {
597  subelem[0]->set_node(0, elem->node_ptr(0));
598  subelem[0]->set_node(1, elem->node_ptr(1));
599  subelem[0]->set_node(2, elem->node_ptr(3));
600  subelem[0]->set_node(3, elem->node_ptr(4));
601  subelem[0]->set_node(4, elem->node_ptr(8));
602  subelem[0]->set_node(5, elem->node_ptr(7));
603 
604  subelem[1]->set_node(0, elem->node_ptr(1));
605  subelem[1]->set_node(1, elem->node_ptr(2));
606  subelem[1]->set_node(2, elem->node_ptr(3));
607  subelem[1]->set_node(3, elem->node_ptr(5));
608  subelem[1]->set_node(4, elem->node_ptr(6));
609  subelem[1]->set_node(5, elem->node_ptr(8));
610  }
611 
612  break;
613  }
614 
615  case HEX8:
616  {
617  BoundaryInfo & boundary_info = mesh.get_boundary_info();
618 
619  // Hexes all split into six tetrahedra
620  subelem[0] = Elem::build(TET4);
621  subelem[1] = Elem::build(TET4);
622  subelem[2] = Elem::build(TET4);
623  subelem[3] = Elem::build(TET4);
624  subelem[4] = Elem::build(TET4);
625  subelem[5] = Elem::build(TET4);
626 
627  // On faces, we choose the node with the highest
628  // global id, and we split on the diagonal which
629  // includes that node. This ensures that (even in
630  // parallel, even on distributed meshes) the same
631  // diagonal split will be chosen for elements on either
632  // side of the same quad face.
633  const unsigned int highest_n = highest_vertex_on(elem);
634 
635  // opposing_node[n] is the local node number of the node
636  // on the farthest corner of a hex8 from local node n
637  static const std::array<unsigned int, 8> opposing_node =
638  {6, 7, 4, 5, 2, 3, 0, 1};
639 
640  static const std::vector<std::vector<unsigned int>> sides_opposing_highest =
641  {{2,3,5},{3,4,5},{1,4,5},{1,2,5},{0,2,3},{0,3,4},{0,1,4},{0,1,2}};
642  static const std::vector<std::vector<unsigned int>> nodes_neighboring_highest =
643  {{1,3,4},{0,2,5},{1,3,6},{0,2,7},{0,5,7},{1,4,6},{2,5,7},{3,4,6}};
644 
645  // Start by looking in three directions away from the
646  // highest-id node. In each direction there will be two
647  // different possibilities for the split depending on
648  // how the opposing face nodes are numbered.
649  //
650  // This is tricky enough that I'm not going to worry
651  // about manually keeping tets oriented; we'll just call
652  // orient() on each as we go.
653 
654  unsigned int next_subelem = 0;
655  for (auto side : sides_opposing_highest[highest_n])
656  {
657  const std::vector<unsigned int> nodes_on_side =
658  elem->nodes_on_side(side);
659 
660  auto [dn, dn2] = split_diagonal(elem, nodes_on_side);
661 
662  unsigned int split_on_neighbor = false;
663  for (auto n : nodes_neighboring_highest[highest_n])
664  if (dn == n || dn2 == n)
665  {
666  split_on_neighbor = true;
667  break;
668  }
669 
670  // Add one or two elements for each opposing side,
671  // depending on whether the diagonal split there
672  // connects to the neighboring diagonal split or
673  // not.
674  if (split_on_neighbor)
675  {
676  subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
677  subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
678  subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
679  for (auto n : nodes_on_side)
680  if (n != dn && n != dn2)
681  {
682  subelem[next_subelem]->set_node(3, elem->node_ptr(n));
683  break;
684  }
685  subelem[next_subelem]->orient(&boundary_info);
686  ++next_subelem;
687 
688  subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
689  subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
690  subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
691  for (auto n : reverse(nodes_on_side))
692  if (n != dn && n != dn2)
693  {
694  subelem[next_subelem]->set_node(3, elem->node_ptr(n));
695  break;
696  }
697  subelem[next_subelem]->orient(&boundary_info);
698  ++next_subelem;
699  }
700  else
701  {
702  subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
703  subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
704  subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
705  for (auto n : nodes_on_side)
706  for (auto n2 : nodes_neighboring_highest[highest_n])
707  if (n == n2)
708  {
709  subelem[next_subelem]->set_node(3, elem->node_ptr(n));
710  goto break_both_loops;
711  }
712 
713  break_both_loops:
714  subelem[next_subelem]->orient(&boundary_info);
715  ++next_subelem;
716  }
717  }
718 
719  // At this point we've created between 3 and 6 tets.
720  // What's left to do depends on how many.
721 
722  // If we just chopped off three vertices into three
723  // tets, then the best way to split this hex would be
724  // the symmetric five-split. Chop off the opposing
725  // vertex too, and then the remaining interior is our
726  // final tet.
727  if (next_subelem == 3)
728  {
729  subelem[next_subelem]->set_node(0, elem->node_ptr(opposing_nodes[highest_n][0]));
730  subelem[next_subelem]->set_node(1, elem->node_ptr(opposing_nodes[highest_n][1]));
731  subelem[next_subelem]->set_node(2, elem->node_ptr(opposing_nodes[highest_n][2]));
732  subelem[next_subelem]->set_node(3, elem->node_ptr(opposing_node[highest_n]));
733  subelem[next_subelem]->orient(&boundary_info);
734  ++next_subelem;
735 
736  subelem[next_subelem]->set_node(0, elem->node_ptr(opposing_nodes[highest_n][0]));
737  subelem[next_subelem]->set_node(1, elem->node_ptr(opposing_nodes[highest_n][1]));
738  subelem[next_subelem]->set_node(2, elem->node_ptr(opposing_nodes[highest_n][2]));
739  subelem[next_subelem]->set_node(3, elem->node_ptr(highest_n));
740  subelem[next_subelem]->orient(&boundary_info);
741  ++next_subelem;
742 
743  // We don't need the 6th tet after all
744  subelem[next_subelem].reset();
745  ++next_subelem;
746  }
747 
748  // If we just chopped off one (or two) vertices into
749  // tets, then the remaining gap is best (or only) filled
750  // by pairing another tet with each.
751  if (next_subelem == 4 ||
752  next_subelem == 5)
753  {
754  for (auto side : sides_opposing_highest[highest_n])
755  {
756  const std::vector<unsigned int> nodes_on_side =
757  elem->nodes_on_side(side);
758 
759  auto [dn, dn2] = split_diagonal(elem, nodes_on_side);
760 
761  unsigned int split_on_neighbor = false;
762  for (auto n : nodes_neighboring_highest[highest_n])
763  if (dn == n || dn2 == n)
764  {
765  split_on_neighbor = true;
766  break;
767  }
768 
769  // The two !split_on_neighbor sides are where we
770  // need the two remaining tets
771  if (!split_on_neighbor)
772  {
773  subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
774  subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
775  subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
776  subelem[next_subelem]->set_node(3, elem->node_ptr(opposing_node[highest_n]));
777  subelem[next_subelem]->orient(&boundary_info);
778  ++next_subelem;
779  }
780  }
781  }
782 
783  // Whether we got there by creating six tets from the
784  // first for loop or by patching up the split afterward,
785  // we should have considered six tets (possibly
786  // including one deleted one...) at this point.
787  libmesh_assert(next_subelem == 6);
788 
789  break;
790  }
791 
792  case PRISM6:
793  {
794  // Prisms all split into three tetrahedra
795  subelem[0] = Elem::build(TET4);
796  subelem[1] = Elem::build(TET4);
797  subelem[2] = Elem::build(TET4);
798 
799  // Triangular faces are not split.
800 
801  // On quad faces, we choose the node with the highest
802  // global id, and we split on the diagonal which
803  // includes that node. This ensures that (even in
804  // parallel, even on distributed meshes) the same
805  // diagonal split will be chosen for elements on either
806  // side of the same quad face. It also ensures that we
807  // always have a mix of "clockwise" and
808  // "counterclockwise" split faces (two of one and one
809  // of the other on each prism; this is useful since the
810  // alternative all-clockwise or all-counterclockwise
811  // face splittings can't be turned into tets without
812  // adding more nodes
813 
814  // Split on 0-4 diagonal
815  if (split_first_diagonal(elem, 0,4, 1,3))
816  {
817  // Split on 0-5 diagonal
818  if (split_first_diagonal(elem, 0,5, 2,3))
819  {
820  // Split on 1-5 diagonal
821  if (split_first_diagonal(elem, 1,5, 2,4))
822  {
823  subelem[0]->set_node(0, elem->node_ptr(0));
824  subelem[0]->set_node(1, elem->node_ptr(4));
825  subelem[0]->set_node(2, elem->node_ptr(5));
826  subelem[0]->set_node(3, elem->node_ptr(3));
827 
828  subelem[1]->set_node(0, elem->node_ptr(0));
829  subelem[1]->set_node(1, elem->node_ptr(4));
830  subelem[1]->set_node(2, elem->node_ptr(1));
831  subelem[1]->set_node(3, elem->node_ptr(5));
832 
833  subelem[2]->set_node(0, elem->node_ptr(0));
834  subelem[2]->set_node(1, elem->node_ptr(1));
835  subelem[2]->set_node(2, elem->node_ptr(2));
836  subelem[2]->set_node(3, elem->node_ptr(5));
837  }
838  else // Split on 2-4 diagonal
839  {
840  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
841 
842  subelem[0]->set_node(0, elem->node_ptr(0));
843  subelem[0]->set_node(1, elem->node_ptr(4));
844  subelem[0]->set_node(2, elem->node_ptr(5));
845  subelem[0]->set_node(3, elem->node_ptr(3));
846 
847  subelem[1]->set_node(0, elem->node_ptr(0));
848  subelem[1]->set_node(1, elem->node_ptr(4));
849  subelem[1]->set_node(2, elem->node_ptr(2));
850  subelem[1]->set_node(3, elem->node_ptr(5));
851 
852  subelem[2]->set_node(0, elem->node_ptr(0));
853  subelem[2]->set_node(1, elem->node_ptr(1));
854  subelem[2]->set_node(2, elem->node_ptr(2));
855  subelem[2]->set_node(3, elem->node_ptr(4));
856  }
857  }
858  else // Split on 2-3 diagonal
859  {
860  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
861 
862  // 0-4 and 2-3 split implies 2-4 split
863  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
864 
865  subelem[0]->set_node(0, elem->node_ptr(0));
866  subelem[0]->set_node(1, elem->node_ptr(4));
867  subelem[0]->set_node(2, elem->node_ptr(2));
868  subelem[0]->set_node(3, elem->node_ptr(3));
869 
870  subelem[1]->set_node(0, elem->node_ptr(3));
871  subelem[1]->set_node(1, elem->node_ptr(4));
872  subelem[1]->set_node(2, elem->node_ptr(2));
873  subelem[1]->set_node(3, elem->node_ptr(5));
874 
875  subelem[2]->set_node(0, elem->node_ptr(0));
876  subelem[2]->set_node(1, elem->node_ptr(1));
877  subelem[2]->set_node(2, elem->node_ptr(2));
878  subelem[2]->set_node(3, elem->node_ptr(4));
879  }
880  }
881  else // Split on 1-3 diagonal
882  {
883  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
884 
885  // Split on 0-5 diagonal
886  if (split_first_diagonal(elem, 0,5, 2,3))
887  {
888  // 1-3 and 0-5 split implies 1-5 split
889  libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
890 
891  subelem[0]->set_node(0, elem->node_ptr(1));
892  subelem[0]->set_node(1, elem->node_ptr(3));
893  subelem[0]->set_node(2, elem->node_ptr(4));
894  subelem[0]->set_node(3, elem->node_ptr(5));
895 
896  subelem[1]->set_node(0, elem->node_ptr(1));
897  subelem[1]->set_node(1, elem->node_ptr(0));
898  subelem[1]->set_node(2, elem->node_ptr(3));
899  subelem[1]->set_node(3, elem->node_ptr(5));
900 
901  subelem[2]->set_node(0, elem->node_ptr(0));
902  subelem[2]->set_node(1, elem->node_ptr(1));
903  subelem[2]->set_node(2, elem->node_ptr(2));
904  subelem[2]->set_node(3, elem->node_ptr(5));
905  }
906  else // Split on 2-3 diagonal
907  {
908  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
909 
910  // Split on 1-5 diagonal
911  if (split_first_diagonal(elem, 1,5, 2,4))
912  {
913  subelem[0]->set_node(0, elem->node_ptr(0));
914  subelem[0]->set_node(1, elem->node_ptr(1));
915  subelem[0]->set_node(2, elem->node_ptr(2));
916  subelem[0]->set_node(3, elem->node_ptr(3));
917 
918  subelem[1]->set_node(0, elem->node_ptr(3));
919  subelem[1]->set_node(1, elem->node_ptr(1));
920  subelem[1]->set_node(2, elem->node_ptr(2));
921  subelem[1]->set_node(3, elem->node_ptr(5));
922 
923  subelem[2]->set_node(0, elem->node_ptr(1));
924  subelem[2]->set_node(1, elem->node_ptr(3));
925  subelem[2]->set_node(2, elem->node_ptr(4));
926  subelem[2]->set_node(3, elem->node_ptr(5));
927  }
928  else // Split on 2-4 diagonal
929  {
930  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
931 
932  subelem[0]->set_node(0, elem->node_ptr(0));
933  subelem[0]->set_node(1, elem->node_ptr(1));
934  subelem[0]->set_node(2, elem->node_ptr(2));
935  subelem[0]->set_node(3, elem->node_ptr(3));
936 
937  subelem[1]->set_node(0, elem->node_ptr(2));
938  subelem[1]->set_node(1, elem->node_ptr(3));
939  subelem[1]->set_node(2, elem->node_ptr(4));
940  subelem[1]->set_node(3, elem->node_ptr(5));
941 
942  subelem[2]->set_node(0, elem->node_ptr(3));
943  subelem[2]->set_node(1, elem->node_ptr(1));
944  subelem[2]->set_node(2, elem->node_ptr(2));
945  subelem[2]->set_node(3, elem->node_ptr(4));
946  }
947  }
948  }
949 
950  break;
951  }
952 
953  case PRISM20:
954  case PRISM21:
955  libmesh_experimental(); // We should upgrade this to TET14...
956  libmesh_fallthrough();
957  case PRISM18:
958  {
959  subelem[0] = Elem::build(TET10);
960  subelem[1] = Elem::build(TET10);
961  subelem[2] = Elem::build(TET10);
962 
963  // Split on 0-4 diagonal
964  if (split_first_diagonal(elem, 0,4, 1,3))
965  {
966  // Split on 0-5 diagonal
967  if (split_first_diagonal(elem, 0,5, 2,3))
968  {
969  // Split on 1-5 diagonal
970  if (split_first_diagonal(elem, 1,5, 2,4))
971  {
972  subelem[0]->set_node(0, elem->node_ptr(0));
973  subelem[0]->set_node(1, elem->node_ptr(4));
974  subelem[0]->set_node(2, elem->node_ptr(5));
975  subelem[0]->set_node(3, elem->node_ptr(3));
976 
977  subelem[0]->set_node(4, elem->node_ptr(15));
978  subelem[0]->set_node(5, elem->node_ptr(13));
979  subelem[0]->set_node(6, elem->node_ptr(17));
980  subelem[0]->set_node(7, elem->node_ptr(9));
981  subelem[0]->set_node(8, elem->node_ptr(12));
982  subelem[0]->set_node(9, elem->node_ptr(14));
983 
984  subelem[1]->set_node(0, elem->node_ptr(0));
985  subelem[1]->set_node(1, elem->node_ptr(4));
986  subelem[1]->set_node(2, elem->node_ptr(1));
987  subelem[1]->set_node(3, elem->node_ptr(5));
988 
989  subelem[1]->set_node(4, elem->node_ptr(15));
990  subelem[1]->set_node(5, elem->node_ptr(10));
991  subelem[1]->set_node(6, elem->node_ptr(6));
992  subelem[1]->set_node(7, elem->node_ptr(17));
993  subelem[1]->set_node(8, elem->node_ptr(13));
994  subelem[1]->set_node(9, elem->node_ptr(16));
995 
996  subelem[2]->set_node(0, elem->node_ptr(0));
997  subelem[2]->set_node(1, elem->node_ptr(1));
998  subelem[2]->set_node(2, elem->node_ptr(2));
999  subelem[2]->set_node(3, elem->node_ptr(5));
1000 
1001  subelem[2]->set_node(4, elem->node_ptr(6));
1002  subelem[2]->set_node(5, elem->node_ptr(7));
1003  subelem[2]->set_node(6, elem->node_ptr(8));
1004  subelem[2]->set_node(7, elem->node_ptr(17));
1005  subelem[2]->set_node(8, elem->node_ptr(16));
1006  subelem[2]->set_node(9, elem->node_ptr(11));
1007  }
1008  else // Split on 2-4 diagonal
1009  {
1010  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1011 
1012  subelem[0]->set_node(0, elem->node_ptr(0));
1013  subelem[0]->set_node(1, elem->node_ptr(4));
1014  subelem[0]->set_node(2, elem->node_ptr(5));
1015  subelem[0]->set_node(3, elem->node_ptr(3));
1016 
1017  subelem[0]->set_node(4, elem->node_ptr(15));
1018  subelem[0]->set_node(5, elem->node_ptr(13));
1019  subelem[0]->set_node(6, elem->node_ptr(17));
1020  subelem[0]->set_node(7, elem->node_ptr(9));
1021  subelem[0]->set_node(8, elem->node_ptr(12));
1022  subelem[0]->set_node(9, elem->node_ptr(14));
1023 
1024  subelem[1]->set_node(0, elem->node_ptr(0));
1025  subelem[1]->set_node(1, elem->node_ptr(4));
1026  subelem[1]->set_node(2, elem->node_ptr(2));
1027  subelem[1]->set_node(3, elem->node_ptr(5));
1028 
1029  subelem[1]->set_node(4, elem->node_ptr(15));
1030  subelem[1]->set_node(5, elem->node_ptr(16));
1031  subelem[1]->set_node(6, elem->node_ptr(8));
1032  subelem[1]->set_node(7, elem->node_ptr(17));
1033  subelem[1]->set_node(8, elem->node_ptr(13));
1034  subelem[1]->set_node(9, elem->node_ptr(11));
1035 
1036  subelem[2]->set_node(0, elem->node_ptr(0));
1037  subelem[2]->set_node(1, elem->node_ptr(1));
1038  subelem[2]->set_node(2, elem->node_ptr(2));
1039  subelem[2]->set_node(3, elem->node_ptr(4));
1040 
1041  subelem[2]->set_node(4, elem->node_ptr(6));
1042  subelem[2]->set_node(5, elem->node_ptr(7));
1043  subelem[2]->set_node(6, elem->node_ptr(8));
1044  subelem[2]->set_node(7, elem->node_ptr(15));
1045  subelem[2]->set_node(8, elem->node_ptr(10));
1046  subelem[2]->set_node(9, elem->node_ptr(16));
1047  }
1048  }
1049  else // Split on 2-3 diagonal
1050  {
1051  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1052 
1053  // 0-4 and 2-3 split implies 2-4 split
1054  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1055 
1056  subelem[0]->set_node(0, elem->node_ptr(0));
1057  subelem[0]->set_node(1, elem->node_ptr(4));
1058  subelem[0]->set_node(2, elem->node_ptr(2));
1059  subelem[0]->set_node(3, elem->node_ptr(3));
1060 
1061  subelem[0]->set_node(4, elem->node_ptr(15));
1062  subelem[0]->set_node(5, elem->node_ptr(16));
1063  subelem[0]->set_node(6, elem->node_ptr(8));
1064  subelem[0]->set_node(7, elem->node_ptr(9));
1065  subelem[0]->set_node(8, elem->node_ptr(12));
1066  subelem[0]->set_node(9, elem->node_ptr(17));
1067 
1068  subelem[1]->set_node(0, elem->node_ptr(3));
1069  subelem[1]->set_node(1, elem->node_ptr(4));
1070  subelem[1]->set_node(2, elem->node_ptr(2));
1071  subelem[1]->set_node(3, elem->node_ptr(5));
1072 
1073  subelem[1]->set_node(4, elem->node_ptr(12));
1074  subelem[1]->set_node(5, elem->node_ptr(16));
1075  subelem[1]->set_node(6, elem->node_ptr(17));
1076  subelem[1]->set_node(7, elem->node_ptr(14));
1077  subelem[1]->set_node(8, elem->node_ptr(13));
1078  subelem[1]->set_node(9, elem->node_ptr(11));
1079 
1080  subelem[2]->set_node(0, elem->node_ptr(0));
1081  subelem[2]->set_node(1, elem->node_ptr(1));
1082  subelem[2]->set_node(2, elem->node_ptr(2));
1083  subelem[2]->set_node(3, elem->node_ptr(4));
1084 
1085  subelem[2]->set_node(4, elem->node_ptr(6));
1086  subelem[2]->set_node(5, elem->node_ptr(7));
1087  subelem[2]->set_node(6, elem->node_ptr(8));
1088  subelem[2]->set_node(7, elem->node_ptr(15));
1089  subelem[2]->set_node(8, elem->node_ptr(10));
1090  subelem[2]->set_node(9, elem->node_ptr(16));
1091  }
1092  }
1093  else // Split on 1-3 diagonal
1094  {
1095  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
1096 
1097  // Split on 0-5 diagonal
1098  if (split_first_diagonal(elem, 0,5, 2,3))
1099  {
1100  // 1-3 and 0-5 split implies 1-5 split
1101  libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
1102 
1103  subelem[0]->set_node(0, elem->node_ptr(1));
1104  subelem[0]->set_node(1, elem->node_ptr(3));
1105  subelem[0]->set_node(2, elem->node_ptr(4));
1106  subelem[0]->set_node(3, elem->node_ptr(5));
1107 
1108  subelem[0]->set_node(4, elem->node_ptr(15));
1109  subelem[0]->set_node(5, elem->node_ptr(12));
1110  subelem[0]->set_node(6, elem->node_ptr(10));
1111  subelem[0]->set_node(7, elem->node_ptr(16));
1112  subelem[0]->set_node(8, elem->node_ptr(14));
1113  subelem[0]->set_node(9, elem->node_ptr(13));
1114 
1115  subelem[1]->set_node(0, elem->node_ptr(1));
1116  subelem[1]->set_node(1, elem->node_ptr(0));
1117  subelem[1]->set_node(2, elem->node_ptr(3));
1118  subelem[1]->set_node(3, elem->node_ptr(5));
1119 
1120  subelem[1]->set_node(4, elem->node_ptr(6));
1121  subelem[1]->set_node(5, elem->node_ptr(9));
1122  subelem[1]->set_node(6, elem->node_ptr(15));
1123  subelem[1]->set_node(7, elem->node_ptr(16));
1124  subelem[1]->set_node(8, elem->node_ptr(17));
1125  subelem[1]->set_node(9, elem->node_ptr(14));
1126 
1127  subelem[2]->set_node(0, elem->node_ptr(0));
1128  subelem[2]->set_node(1, elem->node_ptr(1));
1129  subelem[2]->set_node(2, elem->node_ptr(2));
1130  subelem[2]->set_node(3, elem->node_ptr(5));
1131 
1132  subelem[2]->set_node(4, elem->node_ptr(6));
1133  subelem[2]->set_node(5, elem->node_ptr(7));
1134  subelem[2]->set_node(6, elem->node_ptr(8));
1135  subelem[2]->set_node(7, elem->node_ptr(17));
1136  subelem[2]->set_node(8, elem->node_ptr(16));
1137  subelem[2]->set_node(9, elem->node_ptr(11));
1138  }
1139  else // Split on 2-3 diagonal
1140  {
1141  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1142 
1143  // Split on 1-5 diagonal
1144  if (split_first_diagonal(elem, 1,5, 2,4))
1145  {
1146  subelem[0]->set_node(0, elem->node_ptr(0));
1147  subelem[0]->set_node(1, elem->node_ptr(1));
1148  subelem[0]->set_node(2, elem->node_ptr(2));
1149  subelem[0]->set_node(3, elem->node_ptr(3));
1150 
1151  subelem[0]->set_node(4, elem->node_ptr(6));
1152  subelem[0]->set_node(5, elem->node_ptr(7));
1153  subelem[0]->set_node(6, elem->node_ptr(8));
1154  subelem[0]->set_node(7, elem->node_ptr(9));
1155  subelem[0]->set_node(8, elem->node_ptr(15));
1156  subelem[0]->set_node(9, elem->node_ptr(17));
1157 
1158  subelem[1]->set_node(0, elem->node_ptr(3));
1159  subelem[1]->set_node(1, elem->node_ptr(1));
1160  subelem[1]->set_node(2, elem->node_ptr(2));
1161  subelem[1]->set_node(3, elem->node_ptr(5));
1162 
1163  subelem[1]->set_node(4, elem->node_ptr(15));
1164  subelem[1]->set_node(5, elem->node_ptr(7));
1165  subelem[1]->set_node(6, elem->node_ptr(17));
1166  subelem[1]->set_node(7, elem->node_ptr(14));
1167  subelem[1]->set_node(8, elem->node_ptr(16));
1168  subelem[1]->set_node(9, elem->node_ptr(11));
1169 
1170  subelem[2]->set_node(0, elem->node_ptr(1));
1171  subelem[2]->set_node(1, elem->node_ptr(3));
1172  subelem[2]->set_node(2, elem->node_ptr(4));
1173  subelem[2]->set_node(3, elem->node_ptr(5));
1174 
1175  subelem[2]->set_node(4, elem->node_ptr(15));
1176  subelem[2]->set_node(5, elem->node_ptr(12));
1177  subelem[2]->set_node(6, elem->node_ptr(10));
1178  subelem[2]->set_node(7, elem->node_ptr(16));
1179  subelem[2]->set_node(8, elem->node_ptr(14));
1180  subelem[2]->set_node(9, elem->node_ptr(13));
1181  }
1182  else // Split on 2-4 diagonal
1183  {
1184  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1185 
1186  subelem[0]->set_node(0, elem->node_ptr(0));
1187  subelem[0]->set_node(1, elem->node_ptr(1));
1188  subelem[0]->set_node(2, elem->node_ptr(2));
1189  subelem[0]->set_node(3, elem->node_ptr(3));
1190 
1191  subelem[0]->set_node(4, elem->node_ptr(6));
1192  subelem[0]->set_node(5, elem->node_ptr(7));
1193  subelem[0]->set_node(6, elem->node_ptr(8));
1194  subelem[0]->set_node(7, elem->node_ptr(9));
1195  subelem[0]->set_node(8, elem->node_ptr(15));
1196  subelem[0]->set_node(9, elem->node_ptr(17));
1197 
1198  subelem[1]->set_node(0, elem->node_ptr(2));
1199  subelem[1]->set_node(1, elem->node_ptr(3));
1200  subelem[1]->set_node(2, elem->node_ptr(4));
1201  subelem[1]->set_node(3, elem->node_ptr(5));
1202 
1203  subelem[1]->set_node(4, elem->node_ptr(17));
1204  subelem[1]->set_node(5, elem->node_ptr(12));
1205  subelem[1]->set_node(6, elem->node_ptr(16));
1206  subelem[1]->set_node(7, elem->node_ptr(11));
1207  subelem[1]->set_node(8, elem->node_ptr(14));
1208  subelem[1]->set_node(9, elem->node_ptr(13));
1209 
1210  subelem[2]->set_node(0, elem->node_ptr(3));
1211  subelem[2]->set_node(1, elem->node_ptr(1));
1212  subelem[2]->set_node(2, elem->node_ptr(2));
1213  subelem[2]->set_node(3, elem->node_ptr(4));
1214 
1215  subelem[2]->set_node(4, elem->node_ptr(15));
1216  subelem[2]->set_node(5, elem->node_ptr(7));
1217  subelem[2]->set_node(6, elem->node_ptr(17));
1218  subelem[2]->set_node(7, elem->node_ptr(12));
1219  subelem[2]->set_node(8, elem->node_ptr(10));
1220  subelem[2]->set_node(9, elem->node_ptr(16));
1221  }
1222  }
1223  }
1224 
1225  break;
1226  }
1227 
1228  // No need to split elements that are already simplicial:
1229  case EDGE2:
1230  case EDGE3:
1231  case EDGE4:
1232  case TRI3:
1233  case TRI6:
1234  case TRI7:
1235  case TET4:
1236  case TET10:
1237  case TET14:
1238  case INFEDGE2:
1239  // No way to split infinite quad/prism elements, so
1240  // hopefully no need to
1241  case INFQUAD4:
1242  case INFQUAD6:
1243  case INFPRISM6:
1244  case INFPRISM12:
1245  continue;
1246  // If we're left with an unimplemented hex we're probably
1247  // out of luck. TODO: implement hexes
1248  default:
1249  {
1250  libMesh::err << "Error, encountered unimplemented element "
1251  << Utility::enum_to_string<ElemType>(etype)
1252  << " in MeshTools::Modification::all_tri()..."
1253  << std::endl;
1254  libmesh_not_implemented();
1255  }
1256  } // end switch (etype)
1257 
1258  // Be sure the correct data is set for all subelems.
1259  const unsigned int nei = elem->n_extra_integers();
1260  for (unsigned int i=0; i != max_subelems; ++i)
1261  if (subelem[i]) {
1262  subelem[i]->processor_id() = elem->processor_id();
1263  subelem[i]->subdomain_id() = elem->subdomain_id();
1264 
1265  // Copy any extra element data. Since the subelements
1266  // haven't been added to the mesh yet any allocation has
1267  // to be done manually.
1268  subelem[i]->add_extra_integers(nei);
1269  for (unsigned int ei=0; ei != nei; ++ei)
1270  subelem[ei]->set_extra_integer(ei, elem->get_extra_integer(ei));
1271 
1272 
1273  // Copy any mapping data.
1274  subelem[i]->set_mapping_type(elem->mapping_type());
1275  subelem[i]->set_mapping_data(elem->mapping_data());
1276  }
1277 
1278  // On a mesh with boundary data, we need to move that data to
1279  // the new elements.
1280 
1281  // On a mesh which is distributed, we need to move
1282  // remote_elem links to the new elements.
1283  bool mesh_is_serial = mesh.is_serial();
1284 
1285  if (mesh_has_boundary_data || !mesh_is_serial)
1286  {
1287  // Container to key boundary IDs handed back by the BoundaryInfo object.
1288  std::vector<boundary_id_type> bc_ids;
1289 
1290  for (auto sn : elem->side_index_range())
1291  {
1292  mesh.get_boundary_info().boundary_ids(elem, sn, bc_ids);
1293 
1294  if (bc_ids.empty() && elem->neighbor_ptr(sn) != remote_elem)
1295  continue;
1296 
1297  // Make a sorted list of node ids for elem->side(sn)
1298  elem->build_side_ptr(elem_side, sn);
1299  std::vector<dof_id_type> elem_side_nodes(elem_side->n_nodes());
1300  for (unsigned int esn=0,
1301  n_esn = cast_int<unsigned int>(elem_side_nodes.size());
1302  esn != n_esn; ++esn)
1303  elem_side_nodes[esn] = elem_side->node_id(esn);
1304  std::sort(elem_side_nodes.begin(), elem_side_nodes.end());
1305 
1306  for (unsigned int i=0; i != max_subelems; ++i)
1307  if (subelem[i])
1308  {
1309  for (auto subside : subelem[i]->side_index_range())
1310  {
1311  subelem[i]->build_side_ptr(subside_elem, subside);
1312 
1313  // Make a list of *vertex* node ids for this subside, see if they are all present
1314  // in elem->side(sn). Note 1: we can't just compare elem->key(sn) to
1315  // subelem[i]->key(subside) in the Prism cases, since the new side is
1316  // a different type. Note 2: we only use vertex nodes since, in the future,
1317  // a Hex20 or Prism15's QUAD8 face may be split into two Tri6 faces, and the
1318  // original face will not contain the mid-edge node.
1319  std::vector<dof_id_type> subside_nodes(subside_elem->n_vertices());
1320  for (unsigned int ssn=0,
1321  n_ssn = cast_int<unsigned int>(subside_nodes.size());
1322  ssn != n_ssn; ++ssn)
1323  subside_nodes[ssn] = subside_elem->node_id(ssn);
1324  std::sort(subside_nodes.begin(), subside_nodes.end());
1325 
1326  // std::includes returns true if every element of the second sorted range is
1327  // contained in the first sorted range.
1328  if (std::includes(elem_side_nodes.begin(), elem_side_nodes.end(),
1329  subside_nodes.begin(), subside_nodes.end()))
1330  {
1331  for (const auto & b_id : bc_ids)
1332  if (b_id != BoundaryInfo::invalid_id)
1333  {
1334  new_bndry_ids.push_back(b_id);
1335  new_bndry_elements.push_back(subelem[i].get());
1336  new_bndry_sides.push_back(subside);
1337  }
1338 
1339  // If the original element had a RemoteElem neighbor on side 'sn',
1340  // then the subelem has one on side 'subside'.
1341  if (elem->neighbor_ptr(sn) == remote_elem)
1342  subelem[i]->set_neighbor(subside, const_cast<RemoteElem*>(remote_elem));
1343  }
1344  }
1345  } // end for loop over subelem
1346  } // end for loop over sides
1347 
1348  // Remove the original element from the BoundaryInfo structure.
1349  mesh.get_boundary_info().remove(elem);
1350 
1351  } // end if (mesh_has_boundary_data)
1352 
1353  // Determine new IDs for the split elements which will be
1354  // the same on all processors, therefore keeping the Mesh
1355  // in sync. Note: we offset the new IDs by max_orig_id to
1356  // avoid overwriting any of the original IDs.
1357  for (unsigned int i=0; i != max_subelems; ++i)
1358  if (subelem[i])
1359  {
1360  // Determine new IDs for the split elements which will be
1361  // the same on all processors, therefore keeping the Mesh
1362  // in sync. Note: we offset the new IDs by the max of the
1363  // pre-existing ids to avoid conflicting with originals.
1364  subelem[i]->set_id( max_orig_id + 6*elem->id() + i );
1365 
1366 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1367  subelem[i]->set_unique_id(max_unique_id + max_subelems*elem->unique_id() + i);
1368 #endif
1369 
1370  // Prepare to add the newly-created simplices
1371  new_elements.push_back(std::move(subelem[i]));
1372  }
1373 
1374  // Delete the original element
1375  mesh.delete_elem(elem);
1376  } // End for loop over elements
1377  } // end scope
1378 
1379 
1380  // Now, iterate over the new elements vector, and add them each to
1381  // the Mesh.
1382  for (auto & elem : new_elements)
1383  mesh.add_elem(std::move(elem));
1384 
1385  if (mesh_has_boundary_data)
1386  {
1387  // If the old mesh had boundary data, the new mesh better have
1388  // some. However, we can't assert that the size of
1389  // new_bndry_elements vector is > 0, since we may not have split
1390  // any elements actually on the boundary. We also can't assert
1391  // that the original number of boundary sides is equal to the
1392  // sum of the boundary sides currently in the mesh and the
1393  // newly-added boundary sides, since in 3D, we may have split a
1394  // boundary QUAD into two boundary TRIs. Therefore, we won't be
1395  // too picky about the actual number of BCs, and just assert that
1396  // there are some, somewhere.
1397 #ifndef NDEBUG
1398  bool nbe_nonempty = new_bndry_elements.size();
1399  mesh.comm().max(nbe_nonempty);
1400  libmesh_assert(nbe_nonempty ||
1402 #endif
1403 
1404  // We should also be sure that the lengths of the new boundary data vectors
1405  // are all the same.
1406  libmesh_assert_equal_to (new_bndry_elements.size(), new_bndry_sides.size());
1407  libmesh_assert_equal_to (new_bndry_sides.size(), new_bndry_ids.size());
1408 
1409  // Add the new boundary info to the mesh
1410  for (auto s : index_range(new_bndry_elements))
1411  mesh.get_boundary_info().add_side(new_bndry_elements[s],
1412  new_bndry_sides[s],
1413  new_bndry_ids[s]);
1414  }
1415 
1416  // In a DistributedMesh any newly added ghost node ids may be
1417  // inconsistent, and unique_ids of newly added ghost nodes remain
1418  // unset.
1419  // make_nodes_parallel_consistent() will fix all this.
1420  if (!mesh.is_serial())
1421  {
1422  mesh.comm().max(added_new_ghost_point);
1423 
1424  if (added_new_ghost_point)
1426  }
1427 
1428 
1429 
1430  // Prepare the newly created mesh for use.
1432 
1433  // Let the new_elements and new_bndry_elements vectors go out of scope.
1434 }
1435 
1436 
1438  const unsigned int n_iterations,
1439  const Real power)
1440 {
1444  libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1445 
1446  /*
1447  * Create a quickly-searchable list of boundary nodes.
1448  */
1449  std::unordered_set<dof_id_type> boundary_node_ids =
1451 
1452  // For avoiding extraneous element side allocation
1453  ElemSideBuilder side_builder;
1454 
1455  for (unsigned int iter=0; iter<n_iterations; iter++)
1456  {
1457  /*
1458  * loop over the mesh refinement level
1459  */
1460  unsigned int n_levels = MeshTools::n_levels(mesh);
1461  for (unsigned int refinement_level=0; refinement_level != n_levels;
1462  refinement_level++)
1463  {
1464  // initialize the storage (have to do it on every level to get empty vectors
1465  std::vector<Point> new_positions;
1466  std::vector<Real> weight;
1467  new_positions.resize(mesh.n_nodes());
1468  weight.resize(mesh.n_nodes());
1469 
1470  {
1471  // Loop over the elements to calculate new node positions
1472  for (const auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1473  mesh.level_elements_end(refinement_level)))
1474  {
1475  /*
1476  * We relax all nodes on level 0 first
1477  * If the element is refined (level > 0), we interpolate the
1478  * parents nodes with help of the embedding matrix
1479  */
1480  if (refinement_level == 0)
1481  {
1482  for (auto s : elem->side_index_range())
1483  {
1484  /*
1485  * Only operate on sides which are on the
1486  * boundary or for which the current element's
1487  * id is greater than its neighbor's.
1488  * Sides get only built once.
1489  */
1490  if ((elem->neighbor_ptr(s) != nullptr) &&
1491  (elem->id() > elem->neighbor_ptr(s)->id()))
1492  {
1493  const Elem & side = side_builder(*elem, s);
1494  const Node & node0 = side.node_ref(0);
1495  const Node & node1 = side.node_ref(1);
1496 
1497  Real node_weight = 1.;
1498  // calculate the weight of the nodes
1499  if (power > 0)
1500  {
1501  Point diff = node0-node1;
1502  node_weight = std::pow(diff.norm(), power);
1503  }
1504 
1505  const dof_id_type id0 = node0.id(), id1 = node1.id();
1506  new_positions[id0].add_scaled( node1, node_weight );
1507  new_positions[id1].add_scaled( node0, node_weight );
1508  weight[id0] += node_weight;
1509  weight[id1] += node_weight;
1510  }
1511  } // element neighbor loop
1512  }
1513 #ifdef LIBMESH_ENABLE_AMR
1514  else // refinement_level > 0
1515  {
1516  /*
1517  * Find the positions of the hanging nodes of refined elements.
1518  * We do this by calculating their position based on the parent
1519  * (one level less refined) element, and the embedding matrix
1520  */
1521 
1522  const Elem * parent = elem->parent();
1523 
1524  /*
1525  * find out which child I am
1526  */
1527  unsigned int c = parent->which_child_am_i(elem);
1528  /*
1529  *loop over the childs (that is, the current elements) nodes
1530  */
1531  for (auto nc : elem->node_index_range())
1532  {
1533  /*
1534  * the new position of the node
1535  */
1536  Point point;
1537  for (auto n : parent->node_index_range())
1538  {
1539  /*
1540  * The value from the embedding matrix
1541  */
1542  const Real em_val = parent->embedding_matrix(c,nc,n);
1543 
1544  if (em_val != 0.)
1545  point.add_scaled (parent->point(n), em_val);
1546  }
1547 
1548  const dof_id_type id = elem->node_ptr(nc)->id();
1549  new_positions[id] = point;
1550  weight[id] = 1.;
1551  }
1552  } // if element refinement_level
1553 #endif // #ifdef LIBMESH_ENABLE_AMR
1554 
1555  } // element loop
1556 
1557  /*
1558  * finally reposition the vertex nodes
1559  */
1560  for (auto nid : make_range(mesh.n_nodes()))
1561  if (!boundary_node_ids.count(nid) && weight[nid] > 0.)
1562  mesh.node_ref(nid) = new_positions[nid]/weight[nid];
1563  }
1564 
1565  // Now handle the additional second_order nodes by calculating
1566  // their position based on the vertex positions
1567  // we do a second loop over the level elements
1568  for (auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1569  mesh.level_elements_end(refinement_level)))
1570  {
1571  const unsigned int son_begin = elem->n_vertices();
1572  const unsigned int son_end = elem->n_nodes();
1573  for (unsigned int n=son_begin; n<son_end; n++)
1574  {
1575  const unsigned int n_adjacent_vertices =
1577 
1578  Point point;
1579  for (unsigned int v=0; v<n_adjacent_vertices; v++)
1580  point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
1581 
1582  const dof_id_type id = elem->node_ptr(n)->id();
1583  mesh.node_ref(id) = point/n_adjacent_vertices;
1584  }
1585  }
1586  } // refinement_level loop
1587  } // end iteration
1588 }
1589 
1590 
1591 
1592 #ifdef LIBMESH_ENABLE_AMR
1594 {
1596 
1597  // Algorithm:
1598  // .) For each active element in the mesh: construct a
1599  // copy which is the same in every way *except* it is
1600  // a level 0 element. Store the pointers to these in
1601  // a separate vector. Save any boundary information as well.
1602  // Delete the active element from the mesh.
1603  // .) Loop over all (remaining) elements in the mesh, delete them.
1604  // .) Add the level-0 copies back to the mesh
1605 
1606  // Temporary storage for new element pointers
1607  std::vector<std::unique_ptr<Elem>> new_elements;
1608 
1609  // BoundaryInfo Storage for element ids, sides, and BC ids
1610  std::vector<Elem *> saved_boundary_elements;
1611  std::vector<boundary_id_type> saved_bc_ids;
1612  std::vector<unsigned short int> saved_bc_sides;
1613 
1614  // Container to catch boundary ids passed back by BoundaryInfo
1615  std::vector<boundary_id_type> bc_ids;
1616 
1617  // Reserve a reasonable amt. of space for each
1618  new_elements.reserve(mesh.n_active_elem());
1619  saved_boundary_elements.reserve(mesh.get_boundary_info().n_boundary_conds());
1620  saved_bc_ids.reserve(mesh.get_boundary_info().n_boundary_conds());
1621  saved_bc_sides.reserve(mesh.get_boundary_info().n_boundary_conds());
1622 
1623  for (auto & elem : mesh.active_element_ptr_range())
1624  {
1625  // Make a new element of the same type
1626  auto copy = Elem::build(elem->type());
1627 
1628  // Set node pointers (they still point to nodes in the original mesh)
1629  for (auto n : elem->node_index_range())
1630  copy->set_node(n, elem->node_ptr(n));
1631 
1632  // Copy over ids
1633  copy->processor_id() = elem->processor_id();
1634  copy->subdomain_id() = elem->subdomain_id();
1635 
1636  // Retain the original element's ID(s) as well, otherwise
1637  // the Mesh may try to create them for you...
1638  copy->set_id( elem->id() );
1639 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1640  copy->set_unique_id(elem->unique_id());
1641 #endif
1642 
1643  // This element could have boundary info or DistributedMesh
1644  // remote_elem links as well. We need to save the (elem,
1645  // side, bc_id) triples and those links
1646  for (auto s : elem->side_index_range())
1647  {
1648  if (elem->neighbor_ptr(s) == remote_elem)
1649  copy->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
1650 
1651  mesh.get_boundary_info().boundary_ids(elem, s, bc_ids);
1652  for (const auto & bc_id : bc_ids)
1653  if (bc_id != BoundaryInfo::invalid_id)
1654  {
1655  saved_boundary_elements.push_back(copy.get());
1656  saved_bc_ids.push_back(bc_id);
1657  saved_bc_sides.push_back(s);
1658  }
1659  }
1660 
1661  // Copy any extra element data. Since the copy hasn't been
1662  // added to the mesh yet any allocation has to be done manually.
1663  const unsigned int nei = elem->n_extra_integers();
1664  copy->add_extra_integers(nei);
1665  for (unsigned int i=0; i != nei; ++i)
1666  copy->set_extra_integer(i, elem->get_extra_integer(i));
1667 
1668  // Copy any mapping data.
1669  copy->set_mapping_type(elem->mapping_type());
1670  copy->set_mapping_data(elem->mapping_data());
1671 
1672  // We're done with this element
1673  mesh.delete_elem(elem);
1674 
1675  // But save the copy
1676  new_elements.push_back(std::move(copy));
1677  }
1678 
1679  // Make sure we saved the same number of boundary conditions
1680  // in each vector.
1681  libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
1682  libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
1683 
1684  // Loop again, delete any remaining elements
1685  for (auto & elem : mesh.element_ptr_range())
1686  mesh.delete_elem(elem);
1687 
1688  // Add the copied (now level-0) elements back to the mesh
1689  for (auto & new_elem : new_elements)
1690  {
1691  // Save the original ID, because the act of adding the Elem can
1692  // change new_elem's id!
1693  dof_id_type orig_id = new_elem->id();
1694 
1695  Elem * added_elem = mesh.add_elem(std::move(new_elem));
1696 
1697  // If the Elem, as it was re-added to the mesh, now has a
1698  // different ID (this is unlikely, so it's just an assert)
1699  // the boundary information will no longer be correct.
1700  libmesh_assert_equal_to (orig_id, added_elem->id());
1701 
1702  // Avoid compiler warnings in opt mode.
1703  libmesh_ignore(added_elem, orig_id);
1704  }
1705 
1706  // Finally, also add back the saved boundary information
1707  for (auto e : index_range(saved_boundary_elements))
1708  mesh.get_boundary_info().add_side(saved_boundary_elements[e],
1709  saved_bc_sides[e],
1710  saved_bc_ids[e]);
1711 
1712  // Trim unused and renumber nodes and elements
1714 }
1715 #endif // #ifdef LIBMESH_ENABLE_AMR
1716 
1717 
1718 
1720  const boundary_id_type old_id,
1721  const boundary_id_type new_id)
1722 {
1723  // This is just a shim around the member implementation, now
1724  mesh.get_boundary_info().renumber_id(old_id, new_id);
1725 }
1726 
1727 
1728 
1730  const subdomain_id_type old_id,
1731  const subdomain_id_type new_id)
1732 {
1733  if (old_id == new_id)
1734  {
1735  // If the IDs are the same, this is a no-op.
1736  return;
1737  }
1738 
1739  for (auto & elem : mesh.element_ptr_range())
1740  {
1741  if (elem->subdomain_id() == old_id)
1742  elem->subdomain_id() = new_id;
1743  }
1744 }
1745 
1746 
1747 } // namespace libMesh
OStreamProxy err
unsigned char mapping_data() const
Definition: elem.h:3136
std::size_t n_boundary_conds() const
ElemType
Defines an enum for geometric element types.
const Elem * parent() const
Definition: elem.h:3030
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:629
bool is_prepared() const
Definition: mesh_base.h:198
virtual dof_id_type n_active_elem() const =0
A Node is like a Point, but with more information.
Definition: node.h:52
virtual unique_id_type parallel_max_unique_id() const =0
auto norm() const -> decltype(std::norm(T()))
Definition: type_vector.h:907
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
void orient(BoundaryInfo *boundary_info)
Flips the element (by swapping node and neighbor pointers) to have a mapping Jacobian of opposite sig...
Definition: elem.h:3286
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2710
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i)=0
void scale(MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
Scales the mesh.
void distort(MeshBase &mesh, const Real factor, const bool perturb_boundary=false)
Randomly perturb the nodal locations.
void set_spatial_dimension(unsigned char d)
Sets the "spatial dimension" of the Mesh.
Definition: mesh_base.C:550
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly ecreated (or read) mesh for use.
Definition: mesh_base.C:759
void remove(const Node *node)
Removes the boundary conditions associated with node node, if any exist.
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
unique_id_type unique_id() const
Definition: dof_object.h:844
const Parallel::Communicator & comm() const
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
Fills a user-provided std::vector with the boundary ids associated with Node node.
The libMesh namespace provides an interface to certain functionality in the library.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
void renumber_id(boundary_id_type old_id, boundary_id_type new_id)
Changes all entities (nodes, sides, edges, shellfaces) with boundary id old_id to instead be labeled ...
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.
This is the MeshBase class.
Definition: mesh_base.h:75
void add(const TypeVector< T2 > &)
Add to this vector without creating a temporary.
Definition: type_vector.h:605
std::unordered_set< dof_id_type > find_boundary_nodes(const MeshBase &mesh)
Calling this function on a 2D mesh will convert all the elements to triangles.
Definition: mesh_tools.C:517
void permute_elements(MeshBase &mesh)
Randomly permute the nodal ordering of each element (without twisting the element mapping)...
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 Real embedding_matrix(const unsigned int child_num, const unsigned int child_node_num, const unsigned int parent_node_num) const =0
TensorValue< Real > RealTensorValue
Useful typedefs to allow transparent switching between Real and Complex data types.
void all_tri(MeshBase &mesh)
Subdivides any non-simplex elements in a Mesh to produce simplex (triangular in 2D, tetrahedral in 3D) elements.
void smooth(MeshBase &, unsigned int, Real)
Smooth the mesh with a simple Laplace smoothing algorithm.
virtual bool is_serial() const
Definition: mesh_base.h:211
TypeVector< T > unit() const
Definition: type_vector.h:1104
void libmesh_ignore(const Args &...)
T pow(const T &x)
Definition: utility.h:328
int8_t boundary_id_type
Definition: id_types.h:51
static const boundary_id_type invalid_id
Number used for internal use.
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
ElemMappingType mapping_type() const
Definition: elem.h:3120
This is the MeshCommunication class.
const Node & node_ref(const unsigned int i) const
Definition: elem.h:2529
dof_id_type id() const
Definition: dof_object.h:828
virtual Real hmin() const
Definition: elem.C:702
virtual unsigned int n_nodes() const =0
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:437
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:444
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
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:802
virtual const Node * query_node_ptr(const dof_id_type i) const =0
virtual dof_id_type max_elem_id() const =0
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libmesh_assert(ctx)
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:482
static TensorValue< Real > intrinsic_rotation_matrix(Real phi, Real theta, Real psi)
Generate the intrinsic rotation matrix associated with the provided Euler angles. ...
Definition: tensor_value.h:330
void orient_elements(MeshBase &mesh)
Redo the nodal ordering of each element as necessary to give the element Jacobian a positive orientat...
Helper for building element sides that minimizes the construction of new elements.
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int n) const
Definition: elem.C:3040
auto norm(const T &a) -> decltype(std::abs(a))
Definition: tensor_tools.h:74
void translate(MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.)
Translates the mesh.
RealTensorValue rotate(MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
Rotates the mesh in the xy plane.
unsigned int which_child_am_i(const Elem *e) const
Definition: elem.h:3192
SimpleRange< NodeRefIter > node_ref_range()
Returns a range with all nodes of an element, usable in range-based for loops.
Definition: elem.h:2665
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2598
virtual unsigned int n_vertices() const =0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual const Elem * query_elem_ptr(const dof_id_type i) const =0
subdomain_id_type subdomain_id() const
Definition: elem.h:2582
void max(const T &r, T &o, Request &req) const
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2507
void make_nodes_parallel_consistent(MeshBase &)
Copy processor_ids and ids on ghost nodes from their local processors.
virtual bool is_replicated() const
Definition: mesh_base.h:233
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
Add side side of element number elem with boundary id id to the boundary information data structure...
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const
Definition: elem.C:3048
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
virtual void permute(unsigned int perm_num)=0
Permutes the element (by swapping node and neighbor pointers) according to the specified index...
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
void redistribute(MeshBase &mesh, const FunctionBase< Real > &mapfunc)
Deterministically perturb the nodal locations.
virtual std::unique_ptr< FunctionBase< Output > > clone() const =0
IntRange< unsigned short > node_index_range() const
Definition: elem.h:2683
unsigned int n_extra_integers() const
Returns how many extra integers are associated to the DofObject.
Definition: dof_object.h:1170
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:596
virtual const Point & point(const dof_id_type i) const =0
Defines a dense vector for use in Finite Element-type computations.
Definition: dof_map.h:74
virtual dof_id_type max_node_id() const =0
virtual dof_id_type n_elem() const =0
processor_id_type processor_id() const
processor_id_type processor_id() const
Definition: dof_object.h:905
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:2475
const Point & point(const unsigned int i) const
Definition: elem.h:2453
uint8_t unique_id_type
Definition: id_types.h:86
void flatten(MeshBase &mesh)
Removes all the refinement tree structure of Mesh, leaving only the highest-level (most-refined) elem...
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
virtual unsigned int n_permutations() const =0
Returns the number of independent permutations of element nodes - e.g.
virtual dof_id_type n_nodes() const =0
virtual std::vector< unsigned int > nodes_on_side(const unsigned int) const =0
dof_id_type get_extra_integer(const unsigned int index) const
Gets the value on this object of the extra integer associated with index, which should have been obta...
Definition: dof_object.h:1102
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
uint8_t dof_id_type
Definition: id_types.h:67
const RemoteElem * remote_elem
Definition: remote_elem.C:57