libMesh
mesh_modification.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2026 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/cell_c0polyhedron.h"
34 #include "libmesh/cell_polyhedron.h"
35 #include "libmesh/elem_range.h"
36 #include "libmesh/face_c0polygon.h"
37 #include "libmesh/face_polygon.h"
38 #include "libmesh/face_tri3.h"
39 #include "libmesh/face_tri6.h"
40 #include "libmesh/libmesh_logging.h"
41 #include "libmesh/mesh_communication.h"
42 #include "libmesh/mesh_modification.h"
43 #include "libmesh/mesh_tools.h"
44 #include "libmesh/parallel.h"
45 #include "libmesh/remote_elem.h"
46 #include "libmesh/enum_to_string.h"
47 #include "libmesh/unstructured_mesh.h"
48 #include "libmesh/elem_side_builder.h"
49 #include "libmesh/tensor_value.h"
50 
51 namespace
52 {
53 using namespace libMesh;
54 
55 bool split_first_diagonal(const Elem * elem,
56  unsigned int diag_1_node_1,
57  unsigned int diag_1_node_2,
58  unsigned int diag_2_node_1,
59  unsigned int diag_2_node_2)
60 {
61  return ((elem->node_id(diag_1_node_1) > elem->node_id(diag_2_node_1) &&
62  elem->node_id(diag_1_node_1) > elem->node_id(diag_2_node_2)) ||
63  (elem->node_id(diag_1_node_2) > elem->node_id(diag_2_node_1) &&
64  elem->node_id(diag_1_node_2) > elem->node_id(diag_2_node_2)));
65 }
66 
67 
68 // Return the local index of the vertex on \p elem with the highest
69 // node id.
70 unsigned int highest_vertex_on(const Elem * elem)
71 {
72  unsigned int highest_n = 0;
73  dof_id_type highest_n_id = elem->node_id(0);
74  for (auto n : make_range(1u, elem->n_vertices()))
75  {
76  const dof_id_type n_id = elem->node_id(n);
77  if (n_id > highest_n_id)
78  {
79  highest_n = n;
80  highest_n_id = n_id;
81  }
82  }
83 
84  return highest_n;
85 }
86 
87 
88 static const std::array<std::array<unsigned int, 3>, 8> opposing_nodes =
89 {{ {2,5,7},{3,4,6},{0,5,7},{1,4,6},{1,3,6},{0,2,7},{1,3,4},{0,2,5} }};
90 
91 
92 // Find the highest id on these side nodes of this element
93 std::pair<unsigned int, unsigned int>
94 split_diagonal(const Elem * elem,
95  const std::vector<unsigned int> & nodes_on_side)
96 {
97  libmesh_assert_equal_to(elem->type(), HEX8);
98 
99  unsigned int highest_n = nodes_on_side.front();
100  dof_id_type highest_n_id = elem->node_id(nodes_on_side.front());
101  for (auto n : nodes_on_side)
102  {
103  const dof_id_type n_id = elem->node_id(n);
104  if (n_id > highest_n_id)
105  {
106  highest_n = n;
107  highest_n_id = n_id;
108  }
109  }
110 
111  for (auto n : nodes_on_side)
112  {
113  for (auto n2 : opposing_nodes[highest_n])
114  if (n2 == n)
115  return std::make_pair(highest_n, n2);
116  }
117 
118  libmesh_error();
119 
120  return std::make_pair(libMesh::invalid_uint, libMesh::invalid_uint);
121 }
122 
123 
124 // Reconstruct a C++20 feature in C++14
125 template <typename T>
126 struct reversion_wrapper { T& iterable; };
127 
128 template <typename T>
129 auto begin (reversion_wrapper<T> w) {return std::rbegin(w.iterable);}
130 
131 template <typename T>
132 auto end (reversion_wrapper<T> w) {return std::rend(w.iterable);}
133 
134 template <typename T>
135 reversion_wrapper<T> reverse(T&& iterable) {return {iterable};}
136 
137 }
138 
139 
140 namespace libMesh
141 {
142 
143 
144 // ------------------------------------------------------------
145 // MeshTools::Modification functions for mesh modification
146 void MeshTools::Modification::distort (MeshBase & mesh,
147  const Real factor,
148  const bool perturb_boundary)
149 {
150  libmesh_assert (mesh.n_nodes());
151  libmesh_assert (mesh.n_elem());
152  libmesh_assert ((factor >= 0.) && (factor <= 1.));
153 
154  LOG_SCOPE("distort()", "MeshTools::Modification");
155 
156  // If we are not perturbing boundary nodes, make a
157  // quickly-searchable list of node ids we can check against.
158  std::unordered_set<dof_id_type> boundary_node_ids;
159  if (!perturb_boundary)
160  boundary_node_ids = MeshTools::find_boundary_nodes (mesh);
161 
162  // Now calculate the minimum distance to
163  // neighboring nodes for each node.
164  // hmin holds these distances.
165  std::vector<float> hmin (mesh.max_node_id(),
166  std::numeric_limits<float>::max());
167 
168  for (const auto & elem : mesh.active_element_ptr_range())
169  for (auto & n : elem->node_ref_range())
170  hmin[n.id()] = std::min(hmin[n.id()],
171  static_cast<float>(elem->hmin()));
172 
173  // Now actually move the nodes
174  {
175  const unsigned int seed = 123456;
176 
177  // seed the random number generator.
178  // We'll loop from 1 to n_nodes on every processor, even those
179  // that don't have a particular node, so that the pseudorandom
180  // numbers will be the same everywhere.
181  std::srand(seed);
182 
183  // If the node is on the boundary or
184  // the node is not used by any element (hmin[n]<1.e20)
185  // then we should not move it.
186  // [Note: Testing for (in)equality might be wrong
187  // (different types, namely float and double)]
188  for (auto n : make_range(mesh.max_node_id()))
189  if ((perturb_boundary || !boundary_node_ids.count(n)) && hmin[n] < 1.e20)
190  {
191  // the direction, random but unit normalized
192  Point dir (static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
193  (mesh.mesh_dimension() > 1) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.,
194  ((mesh.mesh_dimension() == 3) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.));
195 
196  dir(0) = (dir(0)-.5)*2.;
197 #if LIBMESH_DIM > 1
198  if (mesh.mesh_dimension() > 1)
199  dir(1) = (dir(1)-.5)*2.;
200 #endif
201 #if LIBMESH_DIM > 2
202  if (mesh.mesh_dimension() == 3)
203  dir(2) = (dir(2)-.5)*2.;
204 #endif
205 
206  dir = dir.unit();
207 
208  Node * node = mesh.query_node_ptr(n);
209  if (!node)
210  continue;
211 
212  (*node)(0) += dir(0)*factor*hmin[n];
213 #if LIBMESH_DIM > 1
214  if (mesh.mesh_dimension() > 1)
215  (*node)(1) += dir(1)*factor*hmin[n];
216 #endif
217 #if LIBMESH_DIM > 2
218  if (mesh.mesh_dimension() == 3)
219  (*node)(2) += dir(2)*factor*hmin[n];
220 #endif
221  }
222  }
223 
224  // We haven't changed any topology, but just changing geometry could
225  // have invalidated a point locator.
226  mesh.clear_point_locator();
227 }
228 
229 
230 
231 void MeshTools::Modification::permute_elements(MeshBase & mesh)
232 {
233  LOG_SCOPE("permute_elements()", "MeshTools::Modification");
234 
235  // We don't yet support doing permute() on a parent element, which
236  // would require us to consistently permute all its children and
237  // give them different local child numbers.
238  unsigned int n_levels = MeshTools::n_levels(mesh);
239  if (n_levels > 1)
240  libmesh_error();
241 
242  const unsigned int seed = 123456;
243 
244  // seed the random number generator.
245  // We'll loop from 1 to max_elem_id on every processor, even those
246  // that don't have a particular element, so that the pseudorandom
247  // numbers will be the same everywhere.
248  std::srand(seed);
249 
250 
251  for (auto e_id : make_range(mesh.max_elem_id()))
252  {
253  int my_rand = std::rand();
254 
255  Elem * elem = mesh.query_elem_ptr(e_id);
256 
257  if (!elem)
258  continue;
259 
260  const unsigned int max_permutation = elem->n_permutations();
261  if (!max_permutation)
262  continue;
263 
264  const unsigned int perm = my_rand % max_permutation;
265 
266  elem->permute(perm);
267  }
268 }
269 
270 
271 void MeshTools::Modification::orient_elements(MeshBase & mesh)
272 {
273  LOG_SCOPE("orient_elements()", "MeshTools::Modification");
274 
275  // We don't yet support doing orient() on a parent element, which
276  // would require us to consistently orient all its children and
277  // give them different local child numbers.
278  unsigned int n_levels = MeshTools::n_levels(mesh);
279  if (n_levels > 1)
280  libmesh_not_implemented_msg("orient_elements() does not support refined meshes");
281 
282  BoundaryInfo & boundary_info = mesh.get_boundary_info();
283  for (auto elem : mesh.element_ptr_range())
284  elem->orient(&boundary_info);
285 }
286 
287 
288 
289 void MeshTools::Modification::redistribute (MeshBase & mesh,
290  const FunctionBase<Real> & mapfunc)
291 {
292  libmesh_assert (mesh.n_nodes());
293  libmesh_assert (mesh.n_elem());
294 
295  LOG_SCOPE("redistribute()", "MeshTools::Modification");
296 
297  DenseVector<Real> output_vec(LIBMESH_DIM);
298 
299  // FIXME - we should thread this later.
300  std::unique_ptr<FunctionBase<Real>> myfunc = mapfunc.clone();
301 
302  for (auto & node : mesh.node_ptr_range())
303  {
304  (*myfunc)(*node, output_vec);
305 
306  (*node)(0) = output_vec(0);
307 #if LIBMESH_DIM > 1
308  (*node)(1) = output_vec(1);
309 #endif
310 #if LIBMESH_DIM > 2
311  (*node)(2) = output_vec(2);
312 #endif
313  }
314 
315  // We haven't changed any topology, but just changing geometry could
316  // have invalidated a point locator.
317  mesh.clear_point_locator();
318 }
319 
320 
321 
322 void MeshTools::Modification::translate (MeshBase & mesh,
323  const Real xt,
324  const Real yt,
325  const Real zt)
326 {
327  const Point p(xt, yt, zt);
328 
329  for (auto & node : mesh.node_ptr_range())
330  *node += p;
331 
332  // We haven't changed any topology, but just changing geometry could
333  // have invalidated a point locator.
334  mesh.clear_point_locator();
335 }
336 
337 
338 // void MeshTools::Modification::rotate2D (MeshBase & mesh,
339 // const Real alpha)
340 // {
341 // libmesh_assert_not_equal_to (mesh.mesh_dimension(), 1);
342 
343 // const Real pi = std::acos(-1);
344 // const Real a = alpha/180.*pi;
345 // for (unsigned int n=0; n<mesh.n_nodes(); n++)
346 // {
347 // const Point p = mesh.node_ref(n);
348 // const Real x = p(0);
349 // const Real y = p(1);
350 // const Real z = p(2);
351 // mesh.node_ref(n) = Point(std::cos(a)*x - std::sin(a)*y,
352 // std::sin(a)*x + std::cos(a)*y,
353 // z);
354 // }
355 
356 // }
357 
358 
359 
360 RealTensorValue
361 MeshTools::Modification::rotate (MeshBase & mesh,
362  const Real phi,
363  const Real theta,
364  const Real psi)
365 {
366  // We won't change any topology, but just changing geometry could
367  // invalidate a point locator.
368  mesh.clear_point_locator();
369 
370 #if LIBMESH_DIM == 3
371  const auto R = RealTensorValue::intrinsic_rotation_matrix(phi, theta, psi);
372 
373  if (theta)
374  mesh.set_spatial_dimension(3);
375 
376  for (auto & node : mesh.node_ptr_range())
377  {
378  Point & pt = *node;
379  pt = R * pt;
380  }
381 
382  return R;
383 
384 #else
385  libmesh_ignore(mesh, phi, theta, psi);
386  libmesh_error_msg("MeshTools::Modification::rotate() requires libMesh to be compiled with LIBMESH_DIM==3");
387  // We'll never get here
388  return RealTensorValue();
389 #endif
390 }
391 
392 
393 void MeshTools::Modification::scale (MeshBase & mesh,
394  const Real xs,
395  const Real ys,
396  const Real zs)
397 {
398  const Real x_scale = xs;
399  Real y_scale = ys;
400  Real z_scale = zs;
401 
402  if (ys == 0.)
403  {
404  libmesh_assert_equal_to (zs, 0.);
405 
406  y_scale = z_scale = x_scale;
407  }
408 
409  // Scale the x coordinate in all dimensions
410  for (auto & node : mesh.node_ptr_range())
411  (*node)(0) *= x_scale;
412 
413  // Only scale the y coordinate in 2 and 3D
414  if (LIBMESH_DIM < 2)
415  return;
416 
417  for (auto & node : mesh.node_ptr_range())
418  (*node)(1) *= y_scale;
419 
420  // Only scale the z coordinate in 3D
421  if (LIBMESH_DIM < 3)
422  return;
423 
424  for (auto & node : mesh.node_ptr_range())
425  (*node)(2) *= z_scale;
426 
427  // We haven't changed any topology, but just changing geometry could
428  // have invalidated a point locator.
429  mesh.clear_point_locator();
430 }
431 
432 
433 
434 void MeshTools::Modification::all_tri (MeshBase & mesh)
435 {
436  if (!mesh.is_replicated() && !mesh.is_prepared())
437  mesh.prepare_for_use();
438 
439  // The number of elements in the original mesh before any additions
440  // or deletions.
441  const dof_id_type n_orig_elem = mesh.n_elem();
442  const dof_id_type max_orig_id = mesh.max_elem_id();
443 
444  // We store pointers to the newly created elements in a vector
445  // until they are ready to be added to the mesh. This is because
446  // adding new elements on the fly can cause reallocation and invalidation
447  // of existing mesh element_iterators.
448  std::vector<std::unique_ptr<Elem>> new_elements;
449 
450  unsigned int max_subelems = 1; // in 1D nothing needs to change
451  if (mesh.mesh_dimension() == 2) // in 2D quads can split into 2 tris
452  max_subelems = 2;
453  if (mesh.mesh_dimension() == 3) // in 3D hexes can split into 6 tets
454  max_subelems = 6;
455 
456  // 2D polygons and 3D polyhedra can be split into an arbitrary
457  // number of triangles/tetrahedra depending on their topology, so we
458  // have to scan the mesh to find the largest split we will need.
459  for (const Elem * elem : mesh.element_ptr_range())
460  {
461  if (const Polygon * poly = dynamic_cast<const Polygon *>(elem))
462  max_subelems = std::max(max_subelems, poly->n_subtriangles());
463  else if (const Polyhedron * polyhedron = dynamic_cast<const Polyhedron *>(elem))
464  max_subelems = std::max(max_subelems, polyhedron->n_subelements());
465  }
466  mesh.comm().max(max_subelems);
467 
468  new_elements.reserve (max_subelems*n_orig_elem);
469 
470  // If the original mesh has *side* boundary data, we carry that over
471  // to the new mesh with triangular elements. We currently only
472  // support bringing over side-based BCs to the all-tri mesh, but
473  // that could probably be extended to node and edge-based BCs as
474  // well.
475  const bool mesh_has_boundary_data = (mesh.get_boundary_info().n_boundary_conds() > 0);
476 
477  // Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
478  std::vector<Elem *> new_bndry_elements;
479  std::vector<unsigned short int> new_bndry_sides;
480  std::vector<boundary_id_type> new_bndry_ids;
481 
482  // We may need to add new points if we run into a 1.5th order
483  // element; if we do that on a DistributedMesh in a ghost element then
484  // we will need to fix their ids / unique_ids
485  bool added_new_ghost_point = false;
486 
487  // Iterate over the elements, splitting:
488  // QUADs into pairs of conforming triangles
489  // PYRAMIDs into pairs of conforming tets,
490  // PRISMs into triplets of conforming tets, and
491  // HEXs into quintets or sextets of conforming tets.
492  // We split on the shortest diagonal to give us better
493  // triangle quality in 2D, and we split based on node ids
494  // to guarantee consistency in 3D.
495  // C0POLYGONs into their sub-triangles
496  // C0POLYHEDRA into their sub-elements (currently only tets)
497 
498  // FIXME: This algorithm does not work on refined grids!
499  {
500 #ifdef LIBMESH_ENABLE_UNIQUE_ID
501  unique_id_type max_unique_id = mesh.parallel_max_unique_id();
502 #endif
503 
504  // For avoiding extraneous allocation when building side elements
505  std::unique_ptr<const Elem> elem_side, subside_elem;
506 
507  for (auto & elem : mesh.element_ptr_range())
508  {
509  const ElemType etype = elem->type();
510 
511  // all_tri currently only works on coarse meshes
512  if (elem->parent())
513  libmesh_not_implemented_msg("Cannot convert a refined element into simplices\n");
514 
515  // The new elements we will split the quad into. Reserving for the maximum
516  // number of sub-elements created for each element
517  std::vector<std::unique_ptr<Elem>> subelem(max_subelems);
518 
519  switch (etype)
520  {
521  case QUAD4:
522  {
523  subelem[0] = Elem::build(TRI3);
524  subelem[1] = Elem::build(TRI3);
525 
526  // Check for possible edge swap
527  if ((elem->point(0) - elem->point(2)).norm() <
528  (elem->point(1) - elem->point(3)).norm())
529  {
530  subelem[0]->set_node(0, elem->node_ptr(0));
531  subelem[0]->set_node(1, elem->node_ptr(1));
532  subelem[0]->set_node(2, elem->node_ptr(2));
533 
534  subelem[1]->set_node(0, elem->node_ptr(0));
535  subelem[1]->set_node(1, elem->node_ptr(2));
536  subelem[1]->set_node(2, elem->node_ptr(3));
537  }
538 
539  else
540  {
541  subelem[0]->set_node(0, elem->node_ptr(0));
542  subelem[0]->set_node(1, elem->node_ptr(1));
543  subelem[0]->set_node(2, elem->node_ptr(3));
544 
545  subelem[1]->set_node(0, elem->node_ptr(1));
546  subelem[1]->set_node(1, elem->node_ptr(2));
547  subelem[1]->set_node(2, elem->node_ptr(3));
548  }
549 
550 
551  break;
552  }
553 
554  case QUAD8:
555  {
556  if (elem->processor_id() != mesh.processor_id())
557  added_new_ghost_point = true;
558 
559  subelem[0] = Elem::build(TRI6);
560  subelem[1] = Elem::build(TRI6);
561 
562  // Add a new node at the center (vertex average) of the element.
563  Node * new_node = mesh.add_point((mesh.point(elem->node_id(0)) +
564  mesh.point(elem->node_id(1)) +
565  mesh.point(elem->node_id(2)) +
566  mesh.point(elem->node_id(3)))/4,
567  DofObject::invalid_id,
568  elem->processor_id());
569 
570  // Check for possible edge swap
571  if ((elem->point(0) - elem->point(2)).norm() <
572  (elem->point(1) - elem->point(3)).norm())
573  {
574  subelem[0]->set_node(0, elem->node_ptr(0));
575  subelem[0]->set_node(1, elem->node_ptr(1));
576  subelem[0]->set_node(2, elem->node_ptr(2));
577  subelem[0]->set_node(3, elem->node_ptr(4));
578  subelem[0]->set_node(4, elem->node_ptr(5));
579  subelem[0]->set_node(5, new_node);
580 
581  subelem[1]->set_node(0, elem->node_ptr(0));
582  subelem[1]->set_node(1, elem->node_ptr(2));
583  subelem[1]->set_node(2, elem->node_ptr(3));
584  subelem[1]->set_node(3, new_node);
585  subelem[1]->set_node(4, elem->node_ptr(6));
586  subelem[1]->set_node(5, elem->node_ptr(7));
587 
588  }
589 
590  else
591  {
592  subelem[0]->set_node(0, elem->node_ptr(3));
593  subelem[0]->set_node(1, elem->node_ptr(0));
594  subelem[0]->set_node(2, elem->node_ptr(1));
595  subelem[0]->set_node(3, elem->node_ptr(7));
596  subelem[0]->set_node(4, elem->node_ptr(4));
597  subelem[0]->set_node(5, new_node);
598 
599  subelem[1]->set_node(0, elem->node_ptr(1));
600  subelem[1]->set_node(1, elem->node_ptr(2));
601  subelem[1]->set_node(2, elem->node_ptr(3));
602  subelem[1]->set_node(3, elem->node_ptr(5));
603  subelem[1]->set_node(4, elem->node_ptr(6));
604  subelem[1]->set_node(5, new_node);
605  }
606 
607  break;
608  }
609 
610  case QUAD9:
611  {
612  subelem[0] = Elem::build(TRI6);
613  subelem[1] = Elem::build(TRI6);
614 
615  // Check for possible edge swap
616  if ((elem->point(0) - elem->point(2)).norm() <
617  (elem->point(1) - elem->point(3)).norm())
618  {
619  subelem[0]->set_node(0, elem->node_ptr(0));
620  subelem[0]->set_node(1, elem->node_ptr(1));
621  subelem[0]->set_node(2, elem->node_ptr(2));
622  subelem[0]->set_node(3, elem->node_ptr(4));
623  subelem[0]->set_node(4, elem->node_ptr(5));
624  subelem[0]->set_node(5, elem->node_ptr(8));
625 
626  subelem[1]->set_node(0, elem->node_ptr(0));
627  subelem[1]->set_node(1, elem->node_ptr(2));
628  subelem[1]->set_node(2, elem->node_ptr(3));
629  subelem[1]->set_node(3, elem->node_ptr(8));
630  subelem[1]->set_node(4, elem->node_ptr(6));
631  subelem[1]->set_node(5, elem->node_ptr(7));
632  }
633 
634  else
635  {
636  subelem[0]->set_node(0, elem->node_ptr(0));
637  subelem[0]->set_node(1, elem->node_ptr(1));
638  subelem[0]->set_node(2, elem->node_ptr(3));
639  subelem[0]->set_node(3, elem->node_ptr(4));
640  subelem[0]->set_node(4, elem->node_ptr(8));
641  subelem[0]->set_node(5, elem->node_ptr(7));
642 
643  subelem[1]->set_node(0, elem->node_ptr(1));
644  subelem[1]->set_node(1, elem->node_ptr(2));
645  subelem[1]->set_node(2, elem->node_ptr(3));
646  subelem[1]->set_node(3, elem->node_ptr(5));
647  subelem[1]->set_node(4, elem->node_ptr(6));
648  subelem[1]->set_node(5, elem->node_ptr(8));
649  }
650 
651  break;
652  }
653 
654  case HEX8:
655  {
656  BoundaryInfo & boundary_info = mesh.get_boundary_info();
657 
658  // Hexes all split into six tetrahedra
659  subelem[0] = Elem::build(TET4);
660  subelem[1] = Elem::build(TET4);
661  subelem[2] = Elem::build(TET4);
662  subelem[3] = Elem::build(TET4);
663  subelem[4] = Elem::build(TET4);
664  subelem[5] = Elem::build(TET4);
665 
666  // On faces, we choose the node with the highest
667  // global id, and we split on the diagonal which
668  // includes that node. This ensures that (even in
669  // parallel, even on distributed meshes) the same
670  // diagonal split will be chosen for elements on either
671  // side of the same quad face.
672  const unsigned int highest_n = highest_vertex_on(elem);
673 
674  // opposing_node[n] is the local node number of the node
675  // on the farthest corner of a hex8 from local node n
676  static const std::array<unsigned int, 8> opposing_node =
677  {6, 7, 4, 5, 2, 3, 0, 1};
678 
679  static const std::vector<std::vector<unsigned int>> sides_opposing_highest =
680  {{2,3,5},{3,4,5},{1,4,5},{1,2,5},{0,2,3},{0,3,4},{0,1,4},{0,1,2}};
681  static const std::vector<std::vector<unsigned int>> nodes_neighboring_highest =
682  {{1,3,4},{0,2,5},{1,3,6},{0,2,7},{0,5,7},{1,4,6},{2,5,7},{3,4,6}};
683 
684  // Start by looking in three directions away from the
685  // highest-id node. In each direction there will be two
686  // different possibilities for the split depending on
687  // how the opposing face nodes are numbered.
688  //
689  // This is tricky enough that I'm not going to worry
690  // about manually keeping tets oriented; we'll just call
691  // orient() on each as we go.
692 
693  unsigned int next_subelem = 0;
694  for (auto side : sides_opposing_highest[highest_n])
695  {
696  const std::vector<unsigned int> nodes_on_side =
697  elem->nodes_on_side(side);
698 
699  auto [dn, dn2] = split_diagonal(elem, nodes_on_side);
700 
701  unsigned int split_on_neighbor = false;
702  for (auto n : nodes_neighboring_highest[highest_n])
703  if (dn == n || dn2 == n)
704  {
705  split_on_neighbor = true;
706  break;
707  }
708 
709  // Add one or two elements for each opposing side,
710  // depending on whether the diagonal split there
711  // connects to the neighboring diagonal split or
712  // not.
713  if (split_on_neighbor)
714  {
715  subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
716  subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
717  subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
718  for (auto n : nodes_on_side)
719  if (n != dn && n != dn2)
720  {
721  subelem[next_subelem]->set_node(3, elem->node_ptr(n));
722  break;
723  }
724  subelem[next_subelem]->orient(&boundary_info);
725  ++next_subelem;
726 
727  subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
728  subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
729  subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
730  for (auto n : reverse(nodes_on_side))
731  if (n != dn && n != dn2)
732  {
733  subelem[next_subelem]->set_node(3, elem->node_ptr(n));
734  break;
735  }
736  subelem[next_subelem]->orient(&boundary_info);
737  ++next_subelem;
738  }
739  else
740  {
741  subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
742  subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
743  subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
744  for (auto n : nodes_on_side)
745  for (auto n2 : nodes_neighboring_highest[highest_n])
746  if (n == n2)
747  {
748  subelem[next_subelem]->set_node(3, elem->node_ptr(n));
749  goto break_both_loops;
750  }
751 
752  break_both_loops:
753  subelem[next_subelem]->orient(&boundary_info);
754  ++next_subelem;
755  }
756  }
757 
758  // At this point we've created between 3 and 6 tets.
759  // What's left to do depends on how many.
760 
761  // If we just chopped off three vertices into three
762  // tets, then the best way to split this hex would be
763  // the symmetric five-split. Chop off the opposing
764  // vertex too, and then the remaining interior is our
765  // final tet.
766  if (next_subelem == 3)
767  {
768  subelem[next_subelem]->set_node(0, elem->node_ptr(opposing_nodes[highest_n][0]));
769  subelem[next_subelem]->set_node(1, elem->node_ptr(opposing_nodes[highest_n][1]));
770  subelem[next_subelem]->set_node(2, elem->node_ptr(opposing_nodes[highest_n][2]));
771  subelem[next_subelem]->set_node(3, elem->node_ptr(opposing_node[highest_n]));
772  subelem[next_subelem]->orient(&boundary_info);
773  ++next_subelem;
774 
775  subelem[next_subelem]->set_node(0, elem->node_ptr(opposing_nodes[highest_n][0]));
776  subelem[next_subelem]->set_node(1, elem->node_ptr(opposing_nodes[highest_n][1]));
777  subelem[next_subelem]->set_node(2, elem->node_ptr(opposing_nodes[highest_n][2]));
778  subelem[next_subelem]->set_node(3, elem->node_ptr(highest_n));
779  subelem[next_subelem]->orient(&boundary_info);
780  ++next_subelem;
781 
782  // We don't need the 6th tet after all
783  subelem[next_subelem].reset();
784  ++next_subelem;
785  }
786 
787  // If we just chopped off one (or two) vertices into
788  // tets, then the remaining gap is best (or only) filled
789  // by pairing another tet with each.
790  if (next_subelem == 4 ||
791  next_subelem == 5)
792  {
793  for (auto side : sides_opposing_highest[highest_n])
794  {
795  const std::vector<unsigned int> nodes_on_side =
796  elem->nodes_on_side(side);
797 
798  auto [dn, dn2] = split_diagonal(elem, nodes_on_side);
799 
800  unsigned int split_on_neighbor = false;
801  for (auto n : nodes_neighboring_highest[highest_n])
802  if (dn == n || dn2 == n)
803  {
804  split_on_neighbor = true;
805  break;
806  }
807 
808  // The two !split_on_neighbor sides are where we
809  // need the two remaining tets
810  if (!split_on_neighbor)
811  {
812  subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
813  subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
814  subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
815  subelem[next_subelem]->set_node(3, elem->node_ptr(opposing_node[highest_n]));
816  subelem[next_subelem]->orient(&boundary_info);
817  ++next_subelem;
818  }
819  }
820  }
821 
822  // Whether we got there by creating six tets from the
823  // first for loop or by patching up the split afterward,
824  // we should have considered six tets (possibly
825  // including one deleted one...) at this point.
826  libmesh_assert(next_subelem == 6);
827 
828  break;
829  }
830 
831  case PRISM6:
832  {
833  // Prisms all split into three tetrahedra
834  subelem[0] = Elem::build(TET4);
835  subelem[1] = Elem::build(TET4);
836  subelem[2] = Elem::build(TET4);
837 
838  // Triangular faces are not split.
839 
840  // On quad faces, we choose the node with the highest
841  // global id, and we split on the diagonal which
842  // includes that node. This ensures that (even in
843  // parallel, even on distributed meshes) the same
844  // diagonal split will be chosen for elements on either
845  // side of the same quad face. It also ensures that we
846  // always have a mix of "clockwise" and
847  // "counterclockwise" split faces (two of one and one
848  // of the other on each prism; this is useful since the
849  // alternative all-clockwise or all-counterclockwise
850  // face splittings can't be turned into tets without
851  // adding more nodes
852 
853  // Split on 0-4 diagonal
854  if (split_first_diagonal(elem, 0,4, 1,3))
855  {
856  // Split on 0-5 diagonal
857  if (split_first_diagonal(elem, 0,5, 2,3))
858  {
859  // Split on 1-5 diagonal
860  if (split_first_diagonal(elem, 1,5, 2,4))
861  {
862  subelem[0]->set_node(0, elem->node_ptr(0));
863  subelem[0]->set_node(1, elem->node_ptr(4));
864  subelem[0]->set_node(2, elem->node_ptr(5));
865  subelem[0]->set_node(3, elem->node_ptr(3));
866 
867  subelem[1]->set_node(0, elem->node_ptr(0));
868  subelem[1]->set_node(1, elem->node_ptr(4));
869  subelem[1]->set_node(2, elem->node_ptr(1));
870  subelem[1]->set_node(3, elem->node_ptr(5));
871 
872  subelem[2]->set_node(0, elem->node_ptr(0));
873  subelem[2]->set_node(1, elem->node_ptr(1));
874  subelem[2]->set_node(2, elem->node_ptr(2));
875  subelem[2]->set_node(3, elem->node_ptr(5));
876  }
877  else // Split on 2-4 diagonal
878  {
879  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
880 
881  subelem[0]->set_node(0, elem->node_ptr(0));
882  subelem[0]->set_node(1, elem->node_ptr(4));
883  subelem[0]->set_node(2, elem->node_ptr(5));
884  subelem[0]->set_node(3, elem->node_ptr(3));
885 
886  subelem[1]->set_node(0, elem->node_ptr(0));
887  subelem[1]->set_node(1, elem->node_ptr(4));
888  subelem[1]->set_node(2, elem->node_ptr(2));
889  subelem[1]->set_node(3, elem->node_ptr(5));
890 
891  subelem[2]->set_node(0, elem->node_ptr(0));
892  subelem[2]->set_node(1, elem->node_ptr(1));
893  subelem[2]->set_node(2, elem->node_ptr(2));
894  subelem[2]->set_node(3, elem->node_ptr(4));
895  }
896  }
897  else // Split on 2-3 diagonal
898  {
899  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
900 
901  // 0-4 and 2-3 split implies 2-4 split
902  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
903 
904  subelem[0]->set_node(0, elem->node_ptr(0));
905  subelem[0]->set_node(1, elem->node_ptr(4));
906  subelem[0]->set_node(2, elem->node_ptr(2));
907  subelem[0]->set_node(3, elem->node_ptr(3));
908 
909  subelem[1]->set_node(0, elem->node_ptr(3));
910  subelem[1]->set_node(1, elem->node_ptr(4));
911  subelem[1]->set_node(2, elem->node_ptr(2));
912  subelem[1]->set_node(3, elem->node_ptr(5));
913 
914  subelem[2]->set_node(0, elem->node_ptr(0));
915  subelem[2]->set_node(1, elem->node_ptr(1));
916  subelem[2]->set_node(2, elem->node_ptr(2));
917  subelem[2]->set_node(3, elem->node_ptr(4));
918  }
919  }
920  else // Split on 1-3 diagonal
921  {
922  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
923 
924  // Split on 0-5 diagonal
925  if (split_first_diagonal(elem, 0,5, 2,3))
926  {
927  // 1-3 and 0-5 split implies 1-5 split
928  libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
929 
930  subelem[0]->set_node(0, elem->node_ptr(1));
931  subelem[0]->set_node(1, elem->node_ptr(3));
932  subelem[0]->set_node(2, elem->node_ptr(4));
933  subelem[0]->set_node(3, elem->node_ptr(5));
934 
935  subelem[1]->set_node(0, elem->node_ptr(1));
936  subelem[1]->set_node(1, elem->node_ptr(0));
937  subelem[1]->set_node(2, elem->node_ptr(3));
938  subelem[1]->set_node(3, elem->node_ptr(5));
939 
940  subelem[2]->set_node(0, elem->node_ptr(0));
941  subelem[2]->set_node(1, elem->node_ptr(1));
942  subelem[2]->set_node(2, elem->node_ptr(2));
943  subelem[2]->set_node(3, elem->node_ptr(5));
944  }
945  else // Split on 2-3 diagonal
946  {
947  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
948 
949  // Split on 1-5 diagonal
950  if (split_first_diagonal(elem, 1,5, 2,4))
951  {
952  subelem[0]->set_node(0, elem->node_ptr(0));
953  subelem[0]->set_node(1, elem->node_ptr(1));
954  subelem[0]->set_node(2, elem->node_ptr(2));
955  subelem[0]->set_node(3, elem->node_ptr(3));
956 
957  subelem[1]->set_node(0, elem->node_ptr(3));
958  subelem[1]->set_node(1, elem->node_ptr(1));
959  subelem[1]->set_node(2, elem->node_ptr(2));
960  subelem[1]->set_node(3, elem->node_ptr(5));
961 
962  subelem[2]->set_node(0, elem->node_ptr(1));
963  subelem[2]->set_node(1, elem->node_ptr(3));
964  subelem[2]->set_node(2, elem->node_ptr(4));
965  subelem[2]->set_node(3, elem->node_ptr(5));
966  }
967  else // Split on 2-4 diagonal
968  {
969  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
970 
971  subelem[0]->set_node(0, elem->node_ptr(0));
972  subelem[0]->set_node(1, elem->node_ptr(1));
973  subelem[0]->set_node(2, elem->node_ptr(2));
974  subelem[0]->set_node(3, elem->node_ptr(3));
975 
976  subelem[1]->set_node(0, elem->node_ptr(2));
977  subelem[1]->set_node(1, elem->node_ptr(3));
978  subelem[1]->set_node(2, elem->node_ptr(4));
979  subelem[1]->set_node(3, elem->node_ptr(5));
980 
981  subelem[2]->set_node(0, elem->node_ptr(3));
982  subelem[2]->set_node(1, elem->node_ptr(1));
983  subelem[2]->set_node(2, elem->node_ptr(2));
984  subelem[2]->set_node(3, elem->node_ptr(4));
985  }
986  }
987  }
988 
989  break;
990  }
991 
992  case PRISM20:
993  case PRISM21:
994  libmesh_experimental(); // We should upgrade this to TET14...
995  libmesh_fallthrough();
996  case PRISM18:
997  {
998  subelem[0] = Elem::build(TET10);
999  subelem[1] = Elem::build(TET10);
1000  subelem[2] = Elem::build(TET10);
1001 
1002  // Split on 0-4 diagonal
1003  if (split_first_diagonal(elem, 0,4, 1,3))
1004  {
1005  // Split on 0-5 diagonal
1006  if (split_first_diagonal(elem, 0,5, 2,3))
1007  {
1008  // Split on 1-5 diagonal
1009  if (split_first_diagonal(elem, 1,5, 2,4))
1010  {
1011  subelem[0]->set_node(0, elem->node_ptr(0));
1012  subelem[0]->set_node(1, elem->node_ptr(4));
1013  subelem[0]->set_node(2, elem->node_ptr(5));
1014  subelem[0]->set_node(3, elem->node_ptr(3));
1015 
1016  subelem[0]->set_node(4, elem->node_ptr(15));
1017  subelem[0]->set_node(5, elem->node_ptr(13));
1018  subelem[0]->set_node(6, elem->node_ptr(17));
1019  subelem[0]->set_node(7, elem->node_ptr(9));
1020  subelem[0]->set_node(8, elem->node_ptr(12));
1021  subelem[0]->set_node(9, elem->node_ptr(14));
1022 
1023  subelem[1]->set_node(0, elem->node_ptr(0));
1024  subelem[1]->set_node(1, elem->node_ptr(4));
1025  subelem[1]->set_node(2, elem->node_ptr(1));
1026  subelem[1]->set_node(3, elem->node_ptr(5));
1027 
1028  subelem[1]->set_node(4, elem->node_ptr(15));
1029  subelem[1]->set_node(5, elem->node_ptr(10));
1030  subelem[1]->set_node(6, elem->node_ptr(6));
1031  subelem[1]->set_node(7, elem->node_ptr(17));
1032  subelem[1]->set_node(8, elem->node_ptr(13));
1033  subelem[1]->set_node(9, elem->node_ptr(16));
1034 
1035  subelem[2]->set_node(0, elem->node_ptr(0));
1036  subelem[2]->set_node(1, elem->node_ptr(1));
1037  subelem[2]->set_node(2, elem->node_ptr(2));
1038  subelem[2]->set_node(3, elem->node_ptr(5));
1039 
1040  subelem[2]->set_node(4, elem->node_ptr(6));
1041  subelem[2]->set_node(5, elem->node_ptr(7));
1042  subelem[2]->set_node(6, elem->node_ptr(8));
1043  subelem[2]->set_node(7, elem->node_ptr(17));
1044  subelem[2]->set_node(8, elem->node_ptr(16));
1045  subelem[2]->set_node(9, elem->node_ptr(11));
1046  }
1047  else // Split on 2-4 diagonal
1048  {
1049  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1050 
1051  subelem[0]->set_node(0, elem->node_ptr(0));
1052  subelem[0]->set_node(1, elem->node_ptr(4));
1053  subelem[0]->set_node(2, elem->node_ptr(5));
1054  subelem[0]->set_node(3, elem->node_ptr(3));
1055 
1056  subelem[0]->set_node(4, elem->node_ptr(15));
1057  subelem[0]->set_node(5, elem->node_ptr(13));
1058  subelem[0]->set_node(6, elem->node_ptr(17));
1059  subelem[0]->set_node(7, elem->node_ptr(9));
1060  subelem[0]->set_node(8, elem->node_ptr(12));
1061  subelem[0]->set_node(9, elem->node_ptr(14));
1062 
1063  subelem[1]->set_node(0, elem->node_ptr(0));
1064  subelem[1]->set_node(1, elem->node_ptr(4));
1065  subelem[1]->set_node(2, elem->node_ptr(2));
1066  subelem[1]->set_node(3, elem->node_ptr(5));
1067 
1068  subelem[1]->set_node(4, elem->node_ptr(15));
1069  subelem[1]->set_node(5, elem->node_ptr(16));
1070  subelem[1]->set_node(6, elem->node_ptr(8));
1071  subelem[1]->set_node(7, elem->node_ptr(17));
1072  subelem[1]->set_node(8, elem->node_ptr(13));
1073  subelem[1]->set_node(9, elem->node_ptr(11));
1074 
1075  subelem[2]->set_node(0, elem->node_ptr(0));
1076  subelem[2]->set_node(1, elem->node_ptr(1));
1077  subelem[2]->set_node(2, elem->node_ptr(2));
1078  subelem[2]->set_node(3, elem->node_ptr(4));
1079 
1080  subelem[2]->set_node(4, elem->node_ptr(6));
1081  subelem[2]->set_node(5, elem->node_ptr(7));
1082  subelem[2]->set_node(6, elem->node_ptr(8));
1083  subelem[2]->set_node(7, elem->node_ptr(15));
1084  subelem[2]->set_node(8, elem->node_ptr(10));
1085  subelem[2]->set_node(9, elem->node_ptr(16));
1086  }
1087  }
1088  else // Split on 2-3 diagonal
1089  {
1090  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1091 
1092  // 0-4 and 2-3 split implies 2-4 split
1093  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1094 
1095  subelem[0]->set_node(0, elem->node_ptr(0));
1096  subelem[0]->set_node(1, elem->node_ptr(4));
1097  subelem[0]->set_node(2, elem->node_ptr(2));
1098  subelem[0]->set_node(3, elem->node_ptr(3));
1099 
1100  subelem[0]->set_node(4, elem->node_ptr(15));
1101  subelem[0]->set_node(5, elem->node_ptr(16));
1102  subelem[0]->set_node(6, elem->node_ptr(8));
1103  subelem[0]->set_node(7, elem->node_ptr(9));
1104  subelem[0]->set_node(8, elem->node_ptr(12));
1105  subelem[0]->set_node(9, elem->node_ptr(17));
1106 
1107  subelem[1]->set_node(0, elem->node_ptr(3));
1108  subelem[1]->set_node(1, elem->node_ptr(4));
1109  subelem[1]->set_node(2, elem->node_ptr(2));
1110  subelem[1]->set_node(3, elem->node_ptr(5));
1111 
1112  subelem[1]->set_node(4, elem->node_ptr(12));
1113  subelem[1]->set_node(5, elem->node_ptr(16));
1114  subelem[1]->set_node(6, elem->node_ptr(17));
1115  subelem[1]->set_node(7, elem->node_ptr(14));
1116  subelem[1]->set_node(8, elem->node_ptr(13));
1117  subelem[1]->set_node(9, elem->node_ptr(11));
1118 
1119  subelem[2]->set_node(0, elem->node_ptr(0));
1120  subelem[2]->set_node(1, elem->node_ptr(1));
1121  subelem[2]->set_node(2, elem->node_ptr(2));
1122  subelem[2]->set_node(3, elem->node_ptr(4));
1123 
1124  subelem[2]->set_node(4, elem->node_ptr(6));
1125  subelem[2]->set_node(5, elem->node_ptr(7));
1126  subelem[2]->set_node(6, elem->node_ptr(8));
1127  subelem[2]->set_node(7, elem->node_ptr(15));
1128  subelem[2]->set_node(8, elem->node_ptr(10));
1129  subelem[2]->set_node(9, elem->node_ptr(16));
1130  }
1131  }
1132  else // Split on 1-3 diagonal
1133  {
1134  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
1135 
1136  // Split on 0-5 diagonal
1137  if (split_first_diagonal(elem, 0,5, 2,3))
1138  {
1139  // 1-3 and 0-5 split implies 1-5 split
1140  libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
1141 
1142  subelem[0]->set_node(0, elem->node_ptr(1));
1143  subelem[0]->set_node(1, elem->node_ptr(3));
1144  subelem[0]->set_node(2, elem->node_ptr(4));
1145  subelem[0]->set_node(3, elem->node_ptr(5));
1146 
1147  subelem[0]->set_node(4, elem->node_ptr(15));
1148  subelem[0]->set_node(5, elem->node_ptr(12));
1149  subelem[0]->set_node(6, elem->node_ptr(10));
1150  subelem[0]->set_node(7, elem->node_ptr(16));
1151  subelem[0]->set_node(8, elem->node_ptr(14));
1152  subelem[0]->set_node(9, elem->node_ptr(13));
1153 
1154  subelem[1]->set_node(0, elem->node_ptr(1));
1155  subelem[1]->set_node(1, elem->node_ptr(0));
1156  subelem[1]->set_node(2, elem->node_ptr(3));
1157  subelem[1]->set_node(3, elem->node_ptr(5));
1158 
1159  subelem[1]->set_node(4, elem->node_ptr(6));
1160  subelem[1]->set_node(5, elem->node_ptr(9));
1161  subelem[1]->set_node(6, elem->node_ptr(15));
1162  subelem[1]->set_node(7, elem->node_ptr(16));
1163  subelem[1]->set_node(8, elem->node_ptr(17));
1164  subelem[1]->set_node(9, elem->node_ptr(14));
1165 
1166  subelem[2]->set_node(0, elem->node_ptr(0));
1167  subelem[2]->set_node(1, elem->node_ptr(1));
1168  subelem[2]->set_node(2, elem->node_ptr(2));
1169  subelem[2]->set_node(3, elem->node_ptr(5));
1170 
1171  subelem[2]->set_node(4, elem->node_ptr(6));
1172  subelem[2]->set_node(5, elem->node_ptr(7));
1173  subelem[2]->set_node(6, elem->node_ptr(8));
1174  subelem[2]->set_node(7, elem->node_ptr(17));
1175  subelem[2]->set_node(8, elem->node_ptr(16));
1176  subelem[2]->set_node(9, elem->node_ptr(11));
1177  }
1178  else // Split on 2-3 diagonal
1179  {
1180  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1181 
1182  // Split on 1-5 diagonal
1183  if (split_first_diagonal(elem, 1,5, 2,4))
1184  {
1185  subelem[0]->set_node(0, elem->node_ptr(0));
1186  subelem[0]->set_node(1, elem->node_ptr(1));
1187  subelem[0]->set_node(2, elem->node_ptr(2));
1188  subelem[0]->set_node(3, elem->node_ptr(3));
1189 
1190  subelem[0]->set_node(4, elem->node_ptr(6));
1191  subelem[0]->set_node(5, elem->node_ptr(7));
1192  subelem[0]->set_node(6, elem->node_ptr(8));
1193  subelem[0]->set_node(7, elem->node_ptr(9));
1194  subelem[0]->set_node(8, elem->node_ptr(15));
1195  subelem[0]->set_node(9, elem->node_ptr(17));
1196 
1197  subelem[1]->set_node(0, elem->node_ptr(3));
1198  subelem[1]->set_node(1, elem->node_ptr(1));
1199  subelem[1]->set_node(2, elem->node_ptr(2));
1200  subelem[1]->set_node(3, elem->node_ptr(5));
1201 
1202  subelem[1]->set_node(4, elem->node_ptr(15));
1203  subelem[1]->set_node(5, elem->node_ptr(7));
1204  subelem[1]->set_node(6, elem->node_ptr(17));
1205  subelem[1]->set_node(7, elem->node_ptr(14));
1206  subelem[1]->set_node(8, elem->node_ptr(16));
1207  subelem[1]->set_node(9, elem->node_ptr(11));
1208 
1209  subelem[2]->set_node(0, elem->node_ptr(1));
1210  subelem[2]->set_node(1, elem->node_ptr(3));
1211  subelem[2]->set_node(2, elem->node_ptr(4));
1212  subelem[2]->set_node(3, elem->node_ptr(5));
1213 
1214  subelem[2]->set_node(4, elem->node_ptr(15));
1215  subelem[2]->set_node(5, elem->node_ptr(12));
1216  subelem[2]->set_node(6, elem->node_ptr(10));
1217  subelem[2]->set_node(7, elem->node_ptr(16));
1218  subelem[2]->set_node(8, elem->node_ptr(14));
1219  subelem[2]->set_node(9, elem->node_ptr(13));
1220  }
1221  else // Split on 2-4 diagonal
1222  {
1223  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1224 
1225  subelem[0]->set_node(0, elem->node_ptr(0));
1226  subelem[0]->set_node(1, elem->node_ptr(1));
1227  subelem[0]->set_node(2, elem->node_ptr(2));
1228  subelem[0]->set_node(3, elem->node_ptr(3));
1229 
1230  subelem[0]->set_node(4, elem->node_ptr(6));
1231  subelem[0]->set_node(5, elem->node_ptr(7));
1232  subelem[0]->set_node(6, elem->node_ptr(8));
1233  subelem[0]->set_node(7, elem->node_ptr(9));
1234  subelem[0]->set_node(8, elem->node_ptr(15));
1235  subelem[0]->set_node(9, elem->node_ptr(17));
1236 
1237  subelem[1]->set_node(0, elem->node_ptr(2));
1238  subelem[1]->set_node(1, elem->node_ptr(3));
1239  subelem[1]->set_node(2, elem->node_ptr(4));
1240  subelem[1]->set_node(3, elem->node_ptr(5));
1241 
1242  subelem[1]->set_node(4, elem->node_ptr(17));
1243  subelem[1]->set_node(5, elem->node_ptr(12));
1244  subelem[1]->set_node(6, elem->node_ptr(16));
1245  subelem[1]->set_node(7, elem->node_ptr(11));
1246  subelem[1]->set_node(8, elem->node_ptr(14));
1247  subelem[1]->set_node(9, elem->node_ptr(13));
1248 
1249  subelem[2]->set_node(0, elem->node_ptr(3));
1250  subelem[2]->set_node(1, elem->node_ptr(1));
1251  subelem[2]->set_node(2, elem->node_ptr(2));
1252  subelem[2]->set_node(3, elem->node_ptr(4));
1253 
1254  subelem[2]->set_node(4, elem->node_ptr(15));
1255  subelem[2]->set_node(5, elem->node_ptr(7));
1256  subelem[2]->set_node(6, elem->node_ptr(17));
1257  subelem[2]->set_node(7, elem->node_ptr(12));
1258  subelem[2]->set_node(8, elem->node_ptr(10));
1259  subelem[2]->set_node(9, elem->node_ptr(16));
1260  }
1261  }
1262  }
1263 
1264  break;
1265  }
1266 
1267  case C0POLYGON:
1268  {
1269  // Split a C0Polygon into the triangles defined by its
1270  // current triangulation. This relies on the user having
1271  // a valid triangulation (the constructor sets a default
1272  // one, and the user can refresh it via retriangulate()
1273  // after moving nodes).
1274  const C0Polygon * polygon = cast_ptr<const C0Polygon *>(elem);
1275  const unsigned int n_subtri = polygon->n_subtriangles();
1276  for (unsigned int t = 0; t != n_subtri; ++t)
1277  {
1278  const std::array<int, 3> tri = polygon->subtriangle(t);
1279  if (tri[0] < 0 || tri[1] < 0 || tri[2] < 0)
1280  libmesh_not_implemented_msg
1281  ("Cannot convert a C0Polygon whose triangulation\n"
1282  "introduces special (non-vertex) points");
1283  subelem[t] = Elem::build(TRI3);
1284  subelem[t]->set_node(0, elem->node_ptr(tri[0]));
1285  subelem[t]->set_node(1, elem->node_ptr(tri[1]));
1286  subelem[t]->set_node(2, elem->node_ptr(tri[2]));
1287  }
1288 
1289  break;
1290  }
1291 
1292  case C0POLYHEDRON:
1293  {
1294  // Split a C0Polyhedron into the tetrahedra defined by its
1295  // current tetrahedralization. If the polyhedron required
1296  // a mid-element node, the user is expected to have added
1297  // that node to the mesh during construction; we just
1298  // reference it via the polyhedron's node pointers.
1299  const C0Polyhedron * polyhedron =
1300  cast_ptr<const C0Polyhedron *>(elem);
1301  const unsigned int n_sub = polyhedron->n_subelements();
1302  for (unsigned int t = 0; t != n_sub; ++t)
1303  {
1304  const std::array<int, 4> tet = polyhedron->subelement(t);
1305  if (tet[0] < 0 || tet[1] < 0 || tet[2] < 0 || tet[3] < 0)
1306  libmesh_not_implemented_msg
1307  ("Cannot convert a C0Polyhedron whose triangulation\n"
1308  "introduces special (non-vertex) points");
1309  subelem[t] = Elem::build(TET4);
1310  subelem[t]->set_node(0, elem->node_ptr(tet[0]));
1311  subelem[t]->set_node(1, elem->node_ptr(tet[1]));
1312  subelem[t]->set_node(2, elem->node_ptr(tet[2]));
1313  subelem[t]->set_node(3, elem->node_ptr(tet[3]));
1314  }
1315  // There is a concern that two neighbor polyhedra might have
1316  // a triangulation of a side that does not match. But the
1317  // default triangulation is based on the side's triangulation
1318  // and the side element is supposed to be shared (that's why
1319  // shared pointers to polygons are used to build the polyhedra).
1320  // So the default one should work.
1321 
1322  break;
1323  }
1324 
1325  // No need to split elements that are already simplicial:
1326  case EDGE2:
1327  case EDGE3:
1328  case EDGE4:
1329  case TRI3:
1330  case TRI6:
1331  case TRI7:
1332  case TET4:
1333  case TET10:
1334  case TET14:
1335  case INFEDGE2:
1336  // No way to split infinite quad/prism elements, so
1337  // hopefully no need to
1338  case INFQUAD4:
1339  case INFQUAD6:
1340  case INFPRISM6:
1341  case INFPRISM12:
1342  continue;
1343  // If we're left with an unimplemented hex we're probably
1344  // out of luck. TODO: implement hexes
1345  default:
1346  {
1347  libMesh::err << "Error, encountered unimplemented element "
1348  << Utility::enum_to_string<ElemType>(etype)
1349  << " in MeshTools::Modification::all_tri()..."
1350  << std::endl;
1351  libmesh_not_implemented();
1352  }
1353  } // end switch (etype)
1354 
1355  // Be sure the correct data is set for all subelems.
1356  const unsigned int nei = elem->n_extra_integers();
1357  for (unsigned int i=0; i != max_subelems; ++i)
1358  if (subelem[i]) {
1359  subelem[i]->processor_id() = elem->processor_id();
1360  subelem[i]->subdomain_id() = elem->subdomain_id();
1361 
1362  // Copy any extra element data. Since the subelements
1363  // haven't been added to the mesh yet any allocation has
1364  // to be done manually.
1365  subelem[i]->add_extra_integers(nei);
1366  for (unsigned int ei=0; ei != nei; ++ei)
1367  subelem[ei]->set_extra_integer(ei, elem->get_extra_integer(ei));
1368 
1369 
1370  // Copy any mapping data.
1371  subelem[i]->set_mapping_type(elem->mapping_type());
1372  subelem[i]->set_mapping_data(elem->mapping_data());
1373  }
1374 
1375  // On a mesh with boundary data, we need to move that data to
1376  // the new elements.
1377 
1378  // On a mesh which is distributed, we need to move
1379  // remote_elem links to the new elements.
1380  bool mesh_is_serial = mesh.is_serial();
1381 
1382  if (mesh_has_boundary_data || !mesh_is_serial)
1383  {
1384  // Container to key boundary IDs handed back by the BoundaryInfo object.
1385  std::vector<boundary_id_type> bc_ids;
1386 
1387  for (auto sn : elem->side_index_range())
1388  {
1389  mesh.get_boundary_info().boundary_ids(elem, sn, bc_ids);
1390 
1391  if (bc_ids.empty() && elem->neighbor_ptr(sn) != remote_elem)
1392  continue;
1393 
1394  // Make a sorted list of node ids for elem->side(sn)
1395  elem->build_side_ptr(elem_side, sn);
1396  std::vector<dof_id_type> elem_side_nodes(elem_side->n_nodes());
1397  for (unsigned int esn=0,
1398  n_esn = cast_int<unsigned int>(elem_side_nodes.size());
1399  esn != n_esn; ++esn)
1400  elem_side_nodes[esn] = elem_side->node_id(esn);
1401  std::sort(elem_side_nodes.begin(), elem_side_nodes.end());
1402 
1403  for (unsigned int i=0; i != max_subelems; ++i)
1404  if (subelem[i])
1405  {
1406  for (auto subside : subelem[i]->side_index_range())
1407  {
1408  subelem[i]->build_side_ptr(subside_elem, subside);
1409 
1410  // Make a list of *vertex* node ids for this subside, see if they are all present
1411  // in elem->side(sn). Note 1: we can't just compare elem->key(sn) to
1412  // subelem[i]->key(subside) in the Prism cases, since the new side is
1413  // a different type. Note 2: we only use vertex nodes since, in the future,
1414  // a Hex20 or Prism15's QUAD8 face may be split into two Tri6 faces, and the
1415  // original face will not contain the mid-edge node.
1416  std::vector<dof_id_type> subside_nodes(subside_elem->n_vertices());
1417  for (unsigned int ssn=0,
1418  n_ssn = cast_int<unsigned int>(subside_nodes.size());
1419  ssn != n_ssn; ++ssn)
1420  subside_nodes[ssn] = subside_elem->node_id(ssn);
1421  std::sort(subside_nodes.begin(), subside_nodes.end());
1422 
1423  // std::includes returns true if every element of the second sorted range is
1424  // contained in the first sorted range.
1425  if (std::includes(elem_side_nodes.begin(), elem_side_nodes.end(),
1426  subside_nodes.begin(), subside_nodes.end()))
1427  {
1428  for (const auto & b_id : bc_ids)
1429  if (b_id != BoundaryInfo::invalid_id)
1430  {
1431  new_bndry_ids.push_back(b_id);
1432  new_bndry_elements.push_back(subelem[i].get());
1433  new_bndry_sides.push_back(subside);
1434  }
1435 
1436  // If the original element had a RemoteElem neighbor on side 'sn',
1437  // then the subelem has one on side 'subside'.
1438  if (elem->neighbor_ptr(sn) == remote_elem)
1439  subelem[i]->set_neighbor(subside, const_cast<RemoteElem*>(remote_elem));
1440  }
1441  }
1442  } // end for loop over subelem
1443  } // end for loop over sides
1444 
1445  // Remove the original element from the BoundaryInfo structure.
1446  mesh.get_boundary_info().remove(elem);
1447 
1448  } // end if (mesh_has_boundary_data)
1449 
1450  // Determine new IDs for the split elements which will be
1451  // the same on all processors, therefore keeping the Mesh
1452  // in sync. Note: we offset the new IDs by max_orig_id to
1453  // avoid overwriting any of the original IDs.
1454  for (unsigned int i=0; i != max_subelems; ++i)
1455  if (subelem[i])
1456  {
1457  // Determine new IDs for the split elements which will be
1458  // the same on all processors, therefore keeping the Mesh
1459  // in sync. Note: we offset the new IDs by the max of the
1460  // pre-existing ids to avoid conflicting with originals.
1461  subelem[i]->set_id( max_orig_id + max_subelems*elem->id() + i );
1462 
1463 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1464  subelem[i]->set_unique_id(max_unique_id + max_subelems*elem->unique_id() + i);
1465 #endif
1466 
1467  // Prepare to add the newly-created simplices
1468  new_elements.push_back(std::move(subelem[i]));
1469  }
1470 
1471  // Delete the original element
1472  mesh.delete_elem(elem);
1473  } // End for loop over elements
1474  } // end scope
1475 
1476 
1477  // Now, iterate over the new elements vector, and add them each to
1478  // the Mesh.
1479  for (auto & elem : new_elements)
1480  mesh.add_elem(std::move(elem));
1481 
1482  if (mesh_has_boundary_data)
1483  {
1484  // If the old mesh had boundary data, the new mesh better have
1485  // some. However, we can't assert that the size of
1486  // new_bndry_elements vector is > 0, since we may not have split
1487  // any elements actually on the boundary. We also can't assert
1488  // that the original number of boundary sides is equal to the
1489  // sum of the boundary sides currently in the mesh and the
1490  // newly-added boundary sides, since in 3D, we may have split a
1491  // boundary QUAD into two boundary TRIs. Therefore, we won't be
1492  // too picky about the actual number of BCs, and just assert that
1493  // there are some, somewhere.
1494 #ifndef NDEBUG
1495  bool nbe_nonempty = new_bndry_elements.size();
1496  mesh.comm().max(nbe_nonempty);
1497  libmesh_assert(nbe_nonempty ||
1498  mesh.get_boundary_info().n_boundary_conds()>0);
1499 #endif
1500 
1501  // We should also be sure that the lengths of the new boundary data vectors
1502  // are all the same.
1503  libmesh_assert_equal_to (new_bndry_elements.size(), new_bndry_sides.size());
1504  libmesh_assert_equal_to (new_bndry_sides.size(), new_bndry_ids.size());
1505 
1506  // Add the new boundary info to the mesh
1507  for (auto s : index_range(new_bndry_elements))
1508  mesh.get_boundary_info().add_side(new_bndry_elements[s],
1509  new_bndry_sides[s],
1510  new_bndry_ids[s]);
1511  }
1512 
1513  // In a DistributedMesh any newly added ghost node ids may be
1514  // inconsistent, and unique_ids of newly added ghost nodes remain
1515  // unset.
1516  // make_nodes_parallel_consistent() will fix all this.
1517  if (!mesh.is_serial())
1518  {
1519  mesh.comm().max(added_new_ghost_point);
1520 
1521  if (added_new_ghost_point)
1522  MeshCommunication().make_nodes_parallel_consistent (mesh);
1523  }
1524 
1525  // Prepare the newly created mesh for use.
1526  mesh.prepare_for_use();
1527 
1528  // Let the new_elements and new_bndry_elements vectors go out of scope.
1529 }
1530 
1531 
1532 
1533 void MeshTools::Modification::all_rbb (MeshBase & mesh)
1534 {
1535  // By default, use 1.0 as the weight on every RATIONAL_BERNSTEIN
1536  // mapped node
1537  const Real default_weight = 1.0;
1538 
1539  const auto weight_index =
1540  (mesh.add_node_datum<Real>("rational_weight", true,
1541  &default_weight));
1542 
1543  mesh.set_default_mapping_type(RATIONAL_BERNSTEIN_MAP);
1544  mesh.set_default_mapping_data(weight_index);
1545 
1546  // Out of loop to reduce heap allocations
1547  std::unique_ptr<Elem> edge_ptr, face_ptr;
1548 
1549  for (auto & elem : mesh.element_ptr_range())
1550  {
1551  if (elem->level())
1552  libmesh_not_implemented_msg
1553  ("all_rbb() currently only supports flat meshes with no refinement levels");
1554 
1555 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1556  if (elem->infinite())
1557  libmesh_not_implemented_msg
1558  ("all_rbb() currently only supports finite geometric elements");
1559 #endif
1560 
1561  elem->set_mapping_type(RATIONAL_BERNSTEIN_MAP);
1562  elem->set_mapping_data(weight_index);
1563 
1564  // Nothing to do unless we have curves to correct
1565  if (elem->default_order() == FIRST)
1566  continue;
1567 
1568  // Modify the center node of an "edge" - possibly an actual edge
1569  // element's node, possibly a center node between points on a
1570  // face's or cell's edge - for RBB interpolation. This should fit
1571  // a circular curve exactly in cases where the original nodes are
1572  // equispaced and the outer nodes' weights are equal, and should
1573  // be a good fit otherwise.
1574  //
1575  // We want to use this to interpolate "internal" conceptual
1576  // edges of a Hex27 too, so we'll handle the cases where w0 and
1577  // w1 aren't 1, as well as the cases where the Nodes n0 and n1
1578  // are already control points which don't match their
1579  // corresponding physical points.
1580  auto make_edge_rbb = [default_weight, weight_index]
1581  (const Node & n0, const Node & n1, Node & n_center,
1582  const Point & p0, const Point & p1)
1583  {
1584  // Skip edges we've already modified; the center node for
1585  // these is no longer at the curve point we wish to
1586  // interpolate, it should already be at the control point that
1587  // accomplishes the interpolation.
1588  const Real old_weight = n_center.get_extra_datum<Real>(weight_index);
1589  if (old_weight != default_weight)
1590  return;
1591 
1592  Point & p2 = n_center;
1593 
1594  const Real w0 = n0.get_extra_datum<Real>(weight_index);
1595  const Real w1 = n1.get_extra_datum<Real>(weight_index);
1596 
1597  const Point e02 = p2-p0,
1598  e21 = p1-p2;
1599  const Real chord_02_len_sq = e02.norm_sq(),
1600  chord_21_len_sq = e21.norm_sq();
1601 
1602  // First find the cosine of phi, the angle between our two
1603  // subchords (turning from the direction of one to the
1604  // direction of the other; this is the supplementary angle to
1605  // the angle at the midpoint). This is the same as half of
1606  // the angle of our circular arc, which nicely enough is also
1607  // the angle we take cos and sec of in NURBS formulae
1608  const Real cos_phi = (e02*e21)/std::sqrt(chord_02_len_sq*chord_21_len_sq);
1609 
1610  // There's a way to do really large arcs using negative
1611  // weights, but we're going to get lousy approximation quality
1612  // from isoparametric elements if we go too low, as well as
1613  // bad numerics here, so let's just disallow it.
1614  if (cos_phi < 0.5)
1615  libmesh_not_implemented_msg
1616  ("all_rbb() is not recommended for extremely sharp curves on one edge");
1617 
1618  const Real w_center = cos_phi*std::sqrt(w0*w1);
1619 
1620  n_center.set_extra_datum<Real>(weight_index, w_center);
1621 
1622  // Now let's get the control point location. This comes from
1623  // a lot of back-and-forth with Gemini, but fortunately I'm
1624  // rewriting it after I've already added unit tests that
1625  // should scream if it's badly wrong.
1626  const Real w_mid = w0/4 + w1/4 + w_center/2;
1627  n_center *= 2*w_mid;
1628  n_center -= (w0 * p0 + w1 * p1)/2;
1629  n_center /= w_center;
1630  };
1631 
1632  auto make_face_rbb = [weight_index] (Elem & face)
1633  {
1634  // Prisms and pyramids may need to skip some faces while
1635  // adjusting others
1636  if (face.type() == TRI6)
1637  return;
1638 
1639  if (face.type() != QUAD9)
1640  libmesh_not_implemented_msg
1641  ("all_rbb() currently only supports mid-face nodes on Quad9 faces");
1642 
1643  // We only use [4,8) but matching indices is nice and stack is
1644  // cheap.
1645  Real w[9];
1646 
1647  for (unsigned int i : make_range(4u, 8u))
1648  w[i] = face.node_ref(i).get_extra_datum<Real>(weight_index);
1649 
1650  // We can't currently handle arbitrary vertex weights
1651 #ifndef NDEBUG
1652  for (unsigned int i : make_range(4u))
1653  libmesh_assert_equal_to
1654  (face.node_ref(i).get_extra_datum<Real>(weight_index), 1);
1655 #endif
1656 
1657  // For the mid-face point, if we want to exactly match
1658  // any cylinders and cones and spheres, we're actually already
1659  // entirely constrained by the other points.
1660  //
1661  // This formula gives the minimum-energy Steiner surface based
1662  // on the outer 8 points.
1663  //
1664  // That's an isogeometric representation of a cylinder aligned
1665  // to either axis, or of a sphere where the quad edges are on
1666  // latitude/longitude lines, or of a cone where two edges are
1667  // segments of cone generating lines and the other two are
1668  // arcs perpendicular to the axis.
1669  //
1670  // It's not perfectly isogeometric for the spheres we generate
1671  // (where the quad edges are all great circles), but it should
1672  // still converge asymptotically faster than non-rational
1673  // quadratic Lagrange.
1674  const Point xi_avg = (face.point(7) + face.point(5))/2;
1675  const Point eta_avg = (face.point(4) + face.point(6))/2;
1676  const Point vertex_avg = (face.point(0) + face.point(1) +
1677  face.point(2) + face.point(3))/4;
1678 
1679  const Real w_xi = (w[7] + w[5])/2;
1680  const Real w_eta = (w[4] + w[6])/2;
1681  const Real w_mid = w_xi * w_eta;
1682 
1683  Node & midnode = face.node_ref(8);
1684  midnode.set_extra_datum<Real>(weight_index, w_mid);
1685  midnode = ((1+w_mid)/(w_xi+w_eta) * (w_xi*xi_avg + w_eta*eta_avg) - vertex_avg)/w_mid;
1686  };
1687 
1688  // If we're on a Hex27, our formula for the mid-volume node
1689  // relies on the locations of the mid-face points. We could
1690  // re-calculate those later but let's just save them now.
1691  Point midfacepts[6];
1692  if (elem->type() == HEX27)
1693  for (auto i : make_range(6))
1694  midfacepts[i] = elem->point(20+i);
1695 
1696  // Check each edge for a curve, and adjust it if needed.
1697  for (auto e : elem->edge_index_range())
1698  {
1699  elem->build_edge_ptr(edge_ptr, e);
1700 
1701  // We should add EDGE4 once we have QUAD16/TRI10/HEX64 to
1702  // use it
1703  if (edge_ptr->type() != EDGE3)
1704  libmesh_not_implemented_msg
1705  ("all_rbb() currently only supports meshes with 2- and/or 3-node edges");
1706 
1707  make_edge_rbb(edge_ptr->node_ref(0), edge_ptr->node_ref(1),
1708  edge_ptr->node_ref(2),
1709  edge_ptr->node_ref(0), edge_ptr->node_ref(1));
1710 
1711  }
1712 
1713  // If we're in 3D, we may have face nodes that also need to be
1714  // adjusted to replace an interpolated curve with a spline
1715  // curve. We know what to do with a quad face, but we'll have
1716  // to scream and die if we see a Tri7 face node.
1717  bool check_face_points = (elem->dim() > 2) &&
1718  (elem->n_nodes() > elem->n_edges() + elem->n_vertices());
1719 
1720  if (check_face_points)
1721  for (auto f : elem->side_index_range())
1722  {
1723  // Prisms and pyramids may need to skip some faces while
1724  // adjusting others
1725  if (elem->side_type(f) == TRI6)
1726  continue;
1727 
1728  elem->build_side_ptr(face_ptr, f);
1729 
1730  make_face_rbb(*face_ptr);
1731  }
1732 
1733  bool check_interior_points =
1734  elem->n_nodes() > elem->n_edges() + elem->n_vertices() + elem->n_faces();
1735 
1736  if (check_interior_points)
1737  {
1738  if (elem->type() == EDGE3)
1739  {
1740  make_edge_rbb(elem->node_ref(0), elem->node_ref(1),
1741  elem->node_ref(2),
1742  elem->node_ref(0), elem->node_ref(1));
1743  }
1744  else if (elem->dim() == 2)
1745  {
1746  make_face_rbb(*elem);
1747  }
1748  else if (elem->type() == HEX27)
1749  {
1750  // We still have the midnode left to go. We want
1751  // something here that will preserve the tensor product
1752  // structure for 2.5D extrusions of IGA faces, but also
1753  // be at least near to the minimum-energy control point
1754  // and weight for general cases. We'll treat opposing
1755  // mid-face nodes as the endpoints of a (more general
1756  // than our edges, since they might have non-1 weights)
1757  // Edge3, and see what we'd need on the midnode to
1758  // interpolate the center point with them. If we've got
1759  // something isogeometric like an extrusion then our
1760  // results should agree; for a quick-but-good output in
1761  // general we'll take an average.
1762  const int opposite_sides[3][2] = {{0,5}, {1,3}, {2,4}};
1763 
1764  Node & midnode = elem->node_ref(26);
1765  const Point original_midpoint = midnode;
1766 
1767  // Averaging in projective space
1768  Point sum_weighted_point = 0;
1769  Real sum_weight = 0;
1770 
1771  for (int i : make_range(3))
1772  {
1773  Node & n0 = elem->node_ref(20+opposite_sides[i][0]);
1774  Node & n1 = elem->node_ref(20+opposite_sides[i][1]);
1775 
1776  make_edge_rbb(n0, n1, midnode,
1777  midfacepts[opposite_sides[i][0]],
1778  midfacepts[opposite_sides[i][1]]);
1779 
1780  const Real midweight =
1781  midnode.get_extra_datum<Real>(weight_index);
1782  sum_weight += midweight;
1783  sum_weighted_point += midweight * midnode;
1784 
1785  // Reset for next run
1786  midnode = original_midpoint;
1787  midnode.set_extra_datum<Real>(weight_index,
1788  default_weight);
1789 
1790  }
1791 
1792  const Real midweight = sum_weight/3;
1793  midnode.set_extra_datum<Real>(weight_index,
1794  midweight);
1795 
1796  midnode = sum_weighted_point / 3 / midweight;
1797  }
1798  else
1799  libmesh_not_implemented_msg
1800  ("all_rbb() doesn't yet support " << elem->type());
1801  }
1802  }
1803 }
1804 
1805 
1806 
1807 void MeshTools::Modification::smooth (MeshBase & mesh,
1808  const unsigned int n_iterations,
1809  const Real power)
1810 {
1814  libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1815 
1816  /*
1817  * Create a quickly-searchable list of boundary nodes.
1818  */
1819  std::unordered_set<dof_id_type> boundary_node_ids =
1820  MeshTools::find_boundary_nodes (mesh);
1821 
1822  // For avoiding extraneous element side allocation
1823  ElemSideBuilder side_builder;
1824 
1825  for (unsigned int iter=0; iter<n_iterations; iter++)
1826  {
1827  /*
1828  * loop over the mesh refinement level
1829  */
1830  unsigned int n_levels = MeshTools::n_levels(mesh);
1831  for (unsigned int refinement_level=0; refinement_level != n_levels;
1832  refinement_level++)
1833  {
1834  // initialize the storage (have to do it on every level to get empty vectors
1835  std::vector<Point> new_positions;
1836  std::vector<Real> weight;
1837  new_positions.resize(mesh.n_nodes());
1838  weight.resize(mesh.n_nodes());
1839 
1840  {
1841  // Loop over the elements to calculate new node positions
1842  for (const auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1843  mesh.level_elements_end(refinement_level)))
1844  {
1845  /*
1846  * We relax all nodes on level 0 first
1847  * If the element is refined (level > 0), we interpolate the
1848  * parents nodes with help of the embedding matrix
1849  */
1850  if (refinement_level == 0)
1851  {
1852  for (auto s : elem->side_index_range())
1853  {
1854  /*
1855  * Only operate on sides which are on the
1856  * boundary or for which the current element's
1857  * id is greater than its neighbor's.
1858  * Sides get only built once.
1859  */
1860  if ((elem->neighbor_ptr(s) != nullptr) &&
1861  (elem->id() > elem->neighbor_ptr(s)->id()))
1862  {
1863  const Elem & side = side_builder(*elem, s);
1864  const Node & node0 = side.node_ref(0);
1865  const Node & node1 = side.node_ref(1);
1866 
1867  Real node_weight = 1.;
1868  // calculate the weight of the nodes
1869  if (power > 0)
1870  {
1871  Point diff = node0-node1;
1872  node_weight = std::pow(diff.norm(), power);
1873  }
1874 
1875  const dof_id_type id0 = node0.id(), id1 = node1.id();
1876  new_positions[id0].add_scaled( node1, node_weight );
1877  new_positions[id1].add_scaled( node0, node_weight );
1878  weight[id0] += node_weight;
1879  weight[id1] += node_weight;
1880  }
1881  } // element neighbor loop
1882  }
1883 #ifdef LIBMESH_ENABLE_AMR
1884  else // refinement_level > 0
1885  {
1886  /*
1887  * Find the positions of the hanging nodes of refined elements.
1888  * We do this by calculating their position based on the parent
1889  * (one level less refined) element, and the embedding matrix
1890  */
1891 
1892  const Elem * parent = elem->parent();
1893 
1894  /*
1895  * find out which child I am
1896  */
1897  unsigned int c = parent->which_child_am_i(elem);
1898  /*
1899  *loop over the childs (that is, the current elements) nodes
1900  */
1901  for (auto nc : elem->node_index_range())
1902  {
1903  /*
1904  * the new position of the node
1905  */
1906  Point point;
1907  for (auto n : parent->node_index_range())
1908  {
1909  /*
1910  * The value from the embedding matrix
1911  */
1912  const Real em_val = parent->embedding_matrix(c,nc,n);
1913 
1914  if (em_val != 0.)
1915  point.add_scaled (parent->point(n), em_val);
1916  }
1917 
1918  const dof_id_type id = elem->node_ptr(nc)->id();
1919  new_positions[id] = point;
1920  weight[id] = 1.;
1921  }
1922  } // if element refinement_level
1923 #endif // #ifdef LIBMESH_ENABLE_AMR
1924 
1925  } // element loop
1926 
1927  /*
1928  * finally reposition the vertex nodes
1929  */
1930  for (auto nid : make_range(mesh.n_nodes()))
1931  if (!boundary_node_ids.count(nid) && weight[nid] > 0.)
1932  mesh.node_ref(nid) = new_positions[nid]/weight[nid];
1933  }
1934 
1935  // Now handle the additional second_order nodes by calculating
1936  // their position based on the vertex positions
1937  // we do a second loop over the level elements
1938  for (auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1939  mesh.level_elements_end(refinement_level)))
1940  {
1941  const unsigned int son_begin = elem->n_vertices();
1942  const unsigned int son_end = elem->n_nodes();
1943  for (unsigned int n=son_begin; n<son_end; n++)
1944  {
1945  const unsigned int n_adjacent_vertices =
1946  elem->n_second_order_adjacent_vertices(n);
1947 
1948  Point point;
1949  for (unsigned int v=0; v<n_adjacent_vertices; v++)
1950  point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
1951 
1952  const dof_id_type id = elem->node_ptr(n)->id();
1953  mesh.node_ref(id) = point/n_adjacent_vertices;
1954  }
1955  }
1956  } // refinement_level loop
1957  } // end iteration
1958 
1959  // We haven't changed any topology, but just changing geometry could
1960  // have invalidated a point locator.
1961  mesh.clear_point_locator();
1962 }
1963 
1964 
1965 
1966 #ifdef LIBMESH_ENABLE_AMR
1967 void MeshTools::Modification::flatten(MeshBase & mesh)
1968 {
1969  libmesh_assert(mesh.is_prepared() || mesh.is_replicated());
1970 
1971  // Algorithm:
1972  // .) For each active element in the mesh: construct a
1973  // copy which is the same in every way *except* it is
1974  // a level 0 element. Store the pointers to these in
1975  // a separate vector. Save any boundary information as well.
1976  // Delete the active element from the mesh.
1977  // .) Loop over all (remaining) elements in the mesh, delete them.
1978  // .) Add the level-0 copies back to the mesh
1979 
1980  // Temporary storage for new element pointers
1981  std::vector<std::unique_ptr<Elem>> new_elements;
1982 
1983  // BoundaryInfo Storage for element ids, sides, and BC ids
1984  std::vector<Elem *> saved_boundary_elements;
1985  std::vector<boundary_id_type> saved_bc_ids;
1986  std::vector<unsigned short int> saved_bc_sides;
1987 
1988  // Container to catch boundary ids passed back by BoundaryInfo
1989  std::vector<boundary_id_type> bc_ids;
1990 
1991  // Reserve a reasonable amt. of space for each
1992  new_elements.reserve(mesh.n_active_elem());
1993  saved_boundary_elements.reserve(mesh.get_boundary_info().n_boundary_conds());
1994  saved_bc_ids.reserve(mesh.get_boundary_info().n_boundary_conds());
1995  saved_bc_sides.reserve(mesh.get_boundary_info().n_boundary_conds());
1996 
1997  for (auto & elem : mesh.active_element_ptr_range())
1998  {
1999  // Make a new element of the same type
2000  auto copy = Elem::build(elem->type());
2001 
2002  // Set node pointers (they still point to nodes in the original mesh)
2003  for (auto n : elem->node_index_range())
2004  copy->set_node(n, elem->node_ptr(n));
2005 
2006  // Copy over ids
2007  copy->processor_id() = elem->processor_id();
2008  copy->subdomain_id() = elem->subdomain_id();
2009 
2010  // Retain the original element's ID(s) as well, otherwise
2011  // the Mesh may try to create them for you...
2012  copy->set_id( elem->id() );
2013 #ifdef LIBMESH_ENABLE_UNIQUE_ID
2014  copy->set_unique_id(elem->unique_id());
2015 #endif
2016 
2017  // This element could have boundary info or DistributedMesh
2018  // remote_elem links as well. We need to save the (elem,
2019  // side, bc_id) triples and those links
2020  for (auto s : elem->side_index_range())
2021  {
2022  if (elem->neighbor_ptr(s) == remote_elem)
2023  copy->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
2024 
2025  mesh.get_boundary_info().boundary_ids(elem, s, bc_ids);
2026  for (const auto & bc_id : bc_ids)
2027  if (bc_id != BoundaryInfo::invalid_id)
2028  {
2029  saved_boundary_elements.push_back(copy.get());
2030  saved_bc_ids.push_back(bc_id);
2031  saved_bc_sides.push_back(s);
2032  }
2033  }
2034 
2035  // Copy any extra element data. Since the copy hasn't been
2036  // added to the mesh yet any allocation has to be done manually.
2037  const unsigned int nei = elem->n_extra_integers();
2038  copy->add_extra_integers(nei);
2039  for (unsigned int i=0; i != nei; ++i)
2040  copy->set_extra_integer(i, elem->get_extra_integer(i));
2041 
2042  // Copy any mapping data.
2043  copy->set_mapping_type(elem->mapping_type());
2044  copy->set_mapping_data(elem->mapping_data());
2045 
2046  // We're done with this element
2047  mesh.delete_elem(elem);
2048 
2049  // But save the copy
2050  new_elements.push_back(std::move(copy));
2051  }
2052 
2053  // Make sure we saved the same number of boundary conditions
2054  // in each vector.
2055  libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
2056  libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
2057 
2058  // Loop again, delete any remaining elements
2059  for (auto & elem : mesh.element_ptr_range())
2060  mesh.delete_elem(elem);
2061 
2062  // Add the copied (now level-0) elements back to the mesh
2063  for (auto & new_elem : new_elements)
2064  {
2065  // Save the original ID, because the act of adding the Elem can
2066  // change new_elem's id!
2067  dof_id_type orig_id = new_elem->id();
2068 
2069  Elem * added_elem = mesh.add_elem(std::move(new_elem));
2070 
2071  // If the Elem, as it was re-added to the mesh, now has a
2072  // different ID (this is unlikely, so it's just an assert)
2073  // the boundary information will no longer be correct.
2074  libmesh_assert_equal_to (orig_id, added_elem->id());
2075 
2076  // Avoid compiler warnings in opt mode.
2077  libmesh_ignore(added_elem, orig_id);
2078  }
2079 
2080  // Finally, also add back the saved boundary information
2081  for (auto e : index_range(saved_boundary_elements))
2082  mesh.get_boundary_info().add_side(saved_boundary_elements[e],
2083  saved_bc_sides[e],
2084  saved_bc_ids[e]);
2085 
2086  // Trim unused and renumber nodes and elements
2087  mesh.prepare_for_use();
2088 }
2089 #endif // #ifdef LIBMESH_ENABLE_AMR
2090 
2091 
2092 
2093 void MeshTools::Modification::change_boundary_id (MeshBase & mesh,
2094  const boundary_id_type old_id,
2095  const boundary_id_type new_id)
2096 {
2097  // This is just a shim around the member implementation, now
2098  mesh.get_boundary_info().renumber_id(old_id, new_id);
2099 }
2100 
2101 
2102 
2103 void MeshTools::Modification::change_subdomain_id (MeshBase & mesh,
2104  const subdomain_id_type old_id,
2105  const subdomain_id_type new_id)
2106 {
2107  if (old_id == new_id)
2108  {
2109  // If the IDs are the same, this is a no-op.
2110  return;
2111  }
2112 
2113  Threads::parallel_for
2114  (mesh.element_stored_range(),
2115  [old_id, new_id](const ElemRange & range)
2116  {
2117  for (Elem * elem : range)
2118  if (elem->subdomain_id() == old_id)
2119  elem->subdomain_id() = new_id;
2120  });
2121 
2122  // We just invalidated mesh.get_subdomain_ids(), but it might not be
2123  // efficient to fix that here.
2124  mesh.unset_has_cached_elem_data();
2125 }
2126 
2127 
2128 } // namespace libMesh
libMesh::err
OStreamProxy err
libMesh::Elem::mapping_data
unsigned char mapping_data() const
Definition: elem.h:3150
libMesh::Threads::parallel_for
void parallel_for(const Range &range, const Body &body, unsigned int n_threads=libMesh::n_threads())
Execute the provided function object in parallel on the specified range.
Definition: threads_none.h:73
libMesh::BoundaryInfo::n_boundary_conds
std::size_t n_boundary_conds() const
Definition: boundary_info.C:2594
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::Elem::parent
const Elem * parent() const
Definition: elem.h:3044
libMesh::TypeVector::add_scaled
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:629
libMesh::MeshBase::is_prepared
bool is_prepared() const
Definition: mesh_base.C:1057
libMesh::MeshBase::n_active_elem
virtual dof_id_type n_active_elem() const =0
libMesh::Node
A Node is like a Point, but with more information.
Definition: node.h:52
libMesh::Polyhedron
The Polyhedron is an element in 3D with an arbitrary number of polygonal faces.
Definition: cell_polyhedron.h:42
libMesh::MeshBase::parallel_max_unique_id
virtual unique_id_type parallel_max_unique_id() const =0
libMesh::TypeVector::norm
auto norm() const
Definition: type_vector.h:908
libMesh::invalid_uint
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:303
libMesh::Elem::orient
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:3300
libMesh::Elem::set_mapping_type
void set_mapping_type(const ElemMappingType type)
Sets the value of the mapping type for the element.
Definition: elem.h:3142
libMesh::Elem::side_type
virtual ElemType side_type(const unsigned int s) const =0
libMesh::Polygon::n_subtriangles
unsigned int n_subtriangles() const
Definition: face_polygon.h:212
libMesh::Elem::side_index_range
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2724
libMesh::Elem::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i)=0
libMesh::MeshTools::Modification::scale
void scale(MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
Scales the mesh.
Definition: mesh_modification.C:393
libMesh::MeshTools::Modification::all_rbb
void all_rbb(MeshBase &mesh)
Converts all element geometric mappings from the default Lagrange to the more flexible Rational-Bezie...
Definition: mesh_modification.C:1533
libMesh::MeshTools::Modification::distort
void distort(MeshBase &mesh, const Real factor, const bool perturb_boundary=false)
Randomly perturb the nodal locations.
Definition: mesh_modification.C:146
libMesh::MeshBase::set_spatial_dimension
void set_spatial_dimension(unsigned char d)
Sets the "spatial dimension" of the Mesh.
Definition: mesh_base.C:613
libMesh::MeshBase::prepare_for_use
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly created (or read) mesh for use.
Definition: mesh_base.C:824
libMesh::BoundaryInfo::remove
void remove(const Node *node)
Removes the boundary conditions associated with node node, if any exist.
Definition: boundary_info.C:1786
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::DofObject::unique_id
unique_id_type unique_id() const
Definition: dof_object.h:835
libMesh::ParallelObject::comm
const Parallel::Communicator & comm() const
Definition: parallel_object.h:97
libMesh::BoundaryInfo::boundary_ids
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.
Definition: boundary_info.C:1348
libMesh::StoredRange
The StoredRange class defines a contiguous, divisible set of objects.
Definition: stored_range.h:54
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:170
libMesh::BoundaryInfo::renumber_id
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 ...
Definition: boundary_info.C:2070
libMesh::MeshBase::add_point
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.
libMesh::MeshBase::unset_has_cached_elem_data
void unset_has_cached_elem_data()
Tells this we have done some operation (e.g.
Definition: mesh_base.h:281
libMesh::MeshBase
This is the MeshBase class.
Definition: mesh_base.h:80
libMesh::TypeVector::add
void add(const TypeVector< T2 > &)
Add to this vector without creating a temporary.
Definition: type_vector.h:605
libMesh::MeshTools::find_boundary_nodes
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:524
libMesh::MeshTools::Modification::permute_elements
void permute_elements(MeshBase &mesh)
Randomly permute the nodal ordering of each element (without twisting the element mapping)...
Definition: mesh_modification.C:231
libMesh::Polygon::subtriangle
virtual std::array< int, 3 > subtriangle(unsigned int i) const
Definition: face_polygon.h:217
libMesh::Elem::edge_index_range
IntRange< unsigned short > edge_index_range() const
Definition: elem.h:2706
libMesh::MeshTools::Modification::change_subdomain_id
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.
Definition: mesh_modification.C:2103
libMesh::Elem::embedding_matrix
virtual Real embedding_matrix(const unsigned int child_num, const unsigned int child_node_num, const unsigned int parent_node_num) const =0
libMesh::RealTensorValue
TensorValue< Real > RealTensorValue
Useful typedefs to allow transparent switching between Real and Complex data types.
Definition: hp_coarsentest.h:48
libMesh::MeshTools::Modification::all_tri
void all_tri(MeshBase &mesh)
Subdivides any non-simplex elements in a Mesh to produce simplex (triangular in 2D, tetrahedral in 3D) elements.
Definition: mesh_modification.C:434
libMesh::MeshTools::Modification::smooth
void smooth(MeshBase &, unsigned int, Real)
Smooth the mesh with a simple Laplace smoothing algorithm.
Definition: mesh_modification.C:1807
libMesh::MeshBase::is_serial
virtual bool is_serial() const
Definition: mesh_base.h:347
libMesh::TypeVector::unit
TypeVector< T > unit() const
Definition: type_vector.h:1141
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:535
libMesh::Utility::pow
T pow(const T &x)
Definition: utility.h:296
libMesh::boundary_id_type
int8_t boundary_id_type
Definition: id_types.h:51
libMesh::BoundaryInfo::invalid_id
static const boundary_id_type invalid_id
Number used for internal use.
Definition: boundary_info.h:912
libMesh::MeshBase::delete_elem
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
libMesh::Elem::mapping_type
ElemMappingType mapping_type() const
Definition: elem.h:3134
libMesh::MeshCommunication
This is the MeshCommunication class.
Definition: mesh_communication.h:53
libMesh::Elem::node_ref
const Node & node_ref(const unsigned int i) const
Definition: elem.h:2535
libMesh::DofObject::id
dof_id_type id() const
Definition: dof_object.h:819
libMesh::Elem::hmin
virtual Real hmin() const
Definition: elem.C:683
libMesh::DofObject::invalid_id
static constexpr dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:473
libMesh::Elem::n_nodes
virtual unsigned int n_nodes() const =0
libMesh::MeshBase::element_stored_range
const ElemRange & element_stored_range()
Definition: mesh_base.C:1913
libMesh::MeshTools::weight
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:444
libMesh::MeshBase::add_elem
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
libMesh::Elem::build
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:442
libMesh::Polyhedron::n_subelements
unsigned int n_subelements() const
Definition: cell_polyhedron.h:281
libMesh::Polyhedron::subelement
virtual std::array< int, 4 > subelement(unsigned int i) const
Definition: cell_polyhedron.h:289
libMesh::Polygon
The Polygon is an element in 2D with an arbitrary (but fixed) number of sides.
Definition: face_polygon.h:39
libMesh::as_range
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
libMesh::MeshTools::Modification::change_boundary_id
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.
Definition: mesh_modification.C:2093
libMesh::MeshTools::n_levels
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:826
libMesh::MeshBase::query_node_ptr
virtual const Node * query_node_ptr(const dof_id_type i) const =0
libMesh::Elem::set_mapping_data
void set_mapping_data(const unsigned char data)
Sets the value of the mapping data for the element.
Definition: elem.h:3158
libMesh::MeshBase::set_default_mapping_type
void set_default_mapping_type(const ElemMappingType type)
Set the default master space to physical space mapping basis functions to be used on newly added elem...
Definition: mesh_base.h:940
libMesh::MeshBase::max_elem_id
virtual dof_id_type max_elem_id() const =0
libMesh::MeshBase::clear_point_locator
void clear_point_locator()
Releases the current PointLocator object.
Definition: mesh_base.C:1859
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::TensorValue::intrinsic_rotation_matrix
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
libMesh::DofObject::set_extra_datum
void set_extra_datum(const unsigned int index, const T value)
Sets the value on this object of the extra datum associated with index, which should have been obtain...
Definition: dof_object.h:1102
libMesh::MeshTools::Modification::orient_elements
void orient_elements(MeshBase &mesh)
Redo the nodal ordering of each element as necessary to give the element Jacobian a positive orientat...
Definition: mesh_modification.C:271
libMesh::ElemSideBuilder
Helper for building element sides that minimizes the construction of new elements.
Definition: elem_side_builder.h:50
libMesh::Elem::n_edges
virtual unsigned int n_edges() const =0
libMesh::Elem::n_second_order_adjacent_vertices
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int n) const
Definition: elem.C:3073
libMesh::MeshTools::Modification::translate
void translate(MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.)
Translates the mesh.
Definition: mesh_modification.C:322
libMesh::MeshBase::add_node_datum
unsigned int add_node_datum(const std::string &name, bool allocate_data=true, const T *default_value=nullptr)
Register a datum (of type T) to be added to each node in the mesh.
Definition: mesh_base.h:2646
libMesh::MeshTools::Modification::rotate
RealTensorValue rotate(MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
Rotates the mesh in the xy plane.
Definition: mesh_modification.C:361
libMesh::Elem::which_child_am_i
unsigned int which_child_am_i(const Elem *e) const
Definition: elem.h:3206
libMesh::Elem::node_ref_range
SimpleRange< NodeRefIter > node_ref_range()
Returns a range with all nodes of an element, usable in range-based for loops.
Definition: elem.h:2679
libMesh::C0Polygon
The C0Polygon is an element in 2D with an arbitrary (but fixed) number of first-order (EDGE2) sides...
Definition: face_c0polygon.h:51
libMesh::Elem::neighbor_ptr
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2612
libMesh::Elem::level
unsigned int level() const
Definition: elem.h:3088
libMesh::Elem::n_vertices
virtual unsigned int n_vertices() const =0
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::MeshBase::query_elem_ptr
virtual const Elem * query_elem_ptr(const dof_id_type i) const =0
libMesh::Elem::subdomain_id
subdomain_id_type subdomain_id() const
Definition: elem.h:2588
TIMPI::Communicator::max
void max(const T &r, T &o, Request &req) const
libMesh::TensorTools::norm
auto norm(const T &a)
Definition: tensor_tools.h:74
libMesh::Elem::dim
virtual unsigned short dim() const =0
libMesh::Elem::node_ptr
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2513
libMesh::MeshCommunication::make_nodes_parallel_consistent
void make_nodes_parallel_consistent(MeshBase &)
Copy processor_ids and ids on ghost nodes from their local processors.
Definition: mesh_communication.C:2042
libMesh::MeshBase::is_replicated
virtual bool is_replicated() const
Definition: mesh_base.h:369
libMesh::BoundaryInfo::add_side
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...
Definition: boundary_info.C:1210
libMesh::Elem::second_order_adjacent_vertex
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const
Definition: elem.C:3081
libMesh::make_range
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:176
libMesh::Elem::permute
virtual void permute(unsigned int perm_num)=0
Permutes the element (by swapping node and neighbor pointers) according to the specified index...
libMesh::MeshBase::mesh_dimension
unsigned int mesh_dimension() const
Definition: mesh_base.C:430
libMesh::MeshTools::Modification::redistribute
void redistribute(MeshBase &mesh, const FunctionBase< Real > &mapfunc)
Deterministically perturb the nodal locations.
Definition: mesh_modification.C:289
libMesh::FunctionBase::clone
virtual std::unique_ptr< FunctionBase< Output > > clone() const =0
libMesh::MeshBase::set_default_mapping_data
void set_default_mapping_data(const unsigned char data)
Set the default master space to physical space mapping basis functions to be used on newly added elem...
Definition: mesh_base.h:958
libMesh::Elem::node_index_range
IntRange< unsigned short > node_index_range() const
Definition: elem.h:2697
libMesh::DofObject::n_extra_integers
unsigned int n_extra_integers() const
Returns how many extra integers are associated to the DofObject.
Definition: dof_object.h:1146
libMesh::MeshBase::node_ref
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:735
libMesh::MeshBase::point
virtual const Point & point(const dof_id_type i) const =0
libMesh::DenseVector
Defines a dense vector for use in Finite Element-type computations.
Definition: dof_map.h:74
libMesh::Elem::n_faces
virtual unsigned int n_faces() const =0
libMesh::DofObject::get_extra_datum
T get_extra_datum(const unsigned int index) const
Gets the value on this object of the extra datum associated with index, which should have been obtain...
Definition: dof_object.h:1122
libMesh::Elem::infinite
virtual bool infinite() const =0
libMesh::C0Polyhedron
The C0Polyhedron is an element in 3D with an arbitrary (but fixed) number of polygonal first-order (C...
Definition: cell_c0polyhedron.h:47
libMesh::MeshBase::max_node_id
virtual dof_id_type max_node_id() const =0
libMesh::MeshBase::n_elem
virtual dof_id_type n_elem() const =0
libMesh::Elem::build_edge_ptr
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i)=0
libMesh::ParallelObject::processor_id
processor_id_type processor_id() const
Definition: parallel_object.h:114
libMesh::Elem::default_order
virtual Order default_order() const =0
libMesh::DofObject::processor_id
processor_id_type processor_id() const
Definition: dof_object.h:881
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::Elem::node_id
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:2481
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:2459
libMesh::unique_id_type
uint8_t unique_id_type
Definition: id_types.h:86
libMesh::MeshTools::Modification::flatten
void flatten(MeshBase &mesh)
Removes all the refinement tree structure of Mesh, leaving only the highest-level (most-refined) elem...
Definition: mesh_modification.C:1967
libMesh::index_range
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:153
libMesh::Elem::n_permutations
virtual unsigned int n_permutations() const =0
Returns the number of independent permutations of element nodes - e.g.
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0
libMesh::Elem::nodes_on_side
virtual std::vector< unsigned int > nodes_on_side(const unsigned int) const =0
libMesh::DofObject::get_extra_integer
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:1078
libMesh::TensorValue
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
Definition: exact_solution.h:46
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
libMesh::remote_elem
const RemoteElem * remote_elem
Definition: remote_elem.C:57
Generated on Wed Jun 24 2026 08:04:14 for libMesh by   doxygen 1.8.14