libMesh
Functions
libMesh::MeshTools::Modification Namespace Reference

Tools for Mesh modification. More...

Functions

void distort (MeshBase &mesh, const Real factor, const bool perturb_boundary=false)
 Randomly perturb the nodal locations. More...
 
void permute_elements (MeshBase &mesh)
 Randomly permute the nodal ordering of each element (without twisting the element mapping). More...
 
void orient_elements (MeshBase &mesh)
 Redo the nodal ordering of each element as necessary to give the element Jacobian a positive orientation. More...
 
void redistribute (MeshBase &mesh, const FunctionBase< Real > &mapfunc)
 Deterministically perturb the nodal locations. More...
 
void translate (MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.)
 Translates the mesh. More...
 
RealTensorValue rotate (MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
 Rotates the mesh in the xy plane. More...
 
void scale (MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
 Scales the mesh. More...
 
void all_tri (MeshBase &mesh)
 Subdivides any non-simplex elements in a Mesh to produce simplex (triangular in 2D, tetrahedral in 3D) elements. More...
 
void all_rbb (MeshBase &mesh)
 Converts all element geometric mappings from the default Lagrange to the more flexible Rational-Bezier-Bernstein. More...
 
void smooth (MeshBase &, unsigned int, Real)
 Smooth the mesh with a simple Laplace smoothing algorithm. More...
 
void flatten (MeshBase &mesh)
 Removes all the refinement tree structure of Mesh, leaving only the highest-level (most-refined) elements. More...
 
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. More...
 
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. More...
 

Detailed Description

Tools for Mesh modification.

Author
Benjamin S. Kirk
Date
2004

Function Documentation

◆ all_rbb()

void libMesh::MeshTools::Modification::all_rbb ( MeshBase &  mesh)

Converts all element geometric mappings from the default Lagrange to the more flexible Rational-Bezier-Bernstein.

When elements have curved edges and/or faces, node weights are chosen so that the new edges interpolate the old edge node locations with a circular arc.

Definition at line 1533 of file mesh_modification.C.

References libMesh::MeshBase::add_node_datum(), libMesh::Elem::build_edge_ptr(), libMesh::Elem::build_side_ptr(), libMesh::Elem::default_order(), libMesh::Elem::dim(), libMesh::EDGE3, libMesh::Elem::edge_index_range(), libMesh::FIRST, libMesh::DofObject::get_extra_datum(), libMesh::HEX27, libMesh::Elem::infinite(), libMesh::Elem::level(), libMesh::make_range(), mesh, libMesh::Elem::n_edges(), libMesh::Elem::n_faces(), libMesh::Elem::n_nodes(), libMesh::Elem::n_vertices(), libMesh::Elem::node_ref(), libMesh::Elem::point(), libMesh::QUAD9, libMesh::RATIONAL_BERNSTEIN_MAP, libMesh::Real, libMesh::MeshBase::set_default_mapping_data(), libMesh::MeshBase::set_default_mapping_type(), libMesh::DofObject::set_extra_datum(), libMesh::Elem::set_mapping_data(), libMesh::Elem::set_mapping_type(), libMesh::Elem::side_index_range(), libMesh::Elem::side_type(), libMesh::TRI6, and libMesh::Elem::type().

Referenced by AllRBBTest::test_box(), AllRBBTest::test_circle(), AllRBBTest::test_cylinder(), and AllRBBTest::test_disk().

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 }
libMesh::Node
A Node is like a Point, but with more information.
Definition: node.h:52
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::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::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::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::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
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::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::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::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39

◆ all_tri()

void libMesh::MeshTools::Modification::all_tri ( MeshBase &  mesh)

Subdivides any non-simplex elements in a Mesh to produce simplex (triangular in 2D, tetrahedral in 3D) elements.

Note
Only supports coarse / unrefined meshes. A uniformly refined mesh can be used only after a flatten() removes its coarser layers.

Definition at line 434 of file mesh_modification.C.

References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::build(), libMesh::Elem::build_side_ptr(), libMesh::C0POLYGON, libMesh::C0POLYHEDRON, libMesh::ParallelObject::comm(), libMesh::MeshBase::delete_elem(), libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, libMesh::err, libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::get_extra_integer(), libMesh::HEX8, libMesh::DofObject::id(), libMesh::index_range(), libMesh::INFEDGE2, libMesh::INFPRISM12, libMesh::INFPRISM6, libMesh::INFQUAD4, libMesh::INFQUAD6, libMesh::DofObject::invalid_id, libMesh::BoundaryInfo::invalid_id, libMesh::MeshBase::is_prepared(), libMesh::MeshBase::is_replicated(), libMesh::MeshBase::is_serial(), libMesh::libmesh_assert(), libMesh::MeshCommunication::make_nodes_parallel_consistent(), libMesh::Elem::mapping_data(), libMesh::Elem::mapping_type(), TIMPI::Communicator::max(), libMesh::MeshBase::max_elem_id(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::MeshBase::n_elem(), libMesh::DofObject::n_extra_integers(), libMesh::Polyhedron::n_subelements(), libMesh::Polygon::n_subtriangles(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_id(), libMesh::Elem::node_ptr(), libMesh::Elem::nodes_on_side(), libMesh::TensorTools::norm(), libMesh::MeshBase::parallel_max_unique_id(), libMesh::Elem::parent(), libMesh::Elem::point(), libMesh::MeshBase::point(), libMesh::MeshBase::prepare_for_use(), libMesh::PRISM18, libMesh::PRISM20, libMesh::PRISM21, libMesh::PRISM6, libMesh::ParallelObject::processor_id(), libMesh::DofObject::processor_id(), libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::remote_elem, libMesh::BoundaryInfo::remove(), libMesh::Elem::side_index_range(), libMesh::Elem::subdomain_id(), libMesh::Polyhedron::subelement(), libMesh::Polygon::subtriangle(), libMesh::TET10, libMesh::TET14, libMesh::TET4, libMesh::TRI3, libMesh::TRI6, libMesh::TRI7, libMesh::Elem::type(), and libMesh::DofObject::unique_id().

Referenced by OverlappingFunctorTest::checkCouplingFunctorTri(), OverlappingFunctorTest::checkCouplingFunctorTriUnifRef(), main(), AllTriTest::test_helper_2D(), AllTriTest::test_helper_3D(), AllTriTest::test_helper_c0polyhedron(), AllTriTest::testAllTriC0Polygon(), AllTriTest::testAllTriC0PolygonOctagon(), and libMesh::MeshTetInterface::volume_to_surface_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 }
libMesh::err
OStreamProxy err
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::MeshBase::is_prepared
bool is_prepared() const
Definition: mesh_base.C:1057
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::Polygon::n_subtriangles
unsigned int n_subtriangles() const
Definition: face_polygon.h:212
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::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::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:170
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::Polygon::subtriangle
virtual std::array< int, 3 > subtriangle(unsigned int i) const
Definition: face_polygon.h:217
libMesh::MeshBase::is_serial
virtual bool is_serial() const
Definition: mesh_base.h:347
libMesh::MeshBase::delete_elem
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
libMesh::MeshCommunication
This is the MeshCommunication class.
Definition: mesh_communication.h:53
libMesh::MeshBase::add_elem
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
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::MeshBase::max_elem_id
virtual dof_id_type max_elem_id() const =0
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::C0Polygon
The C0Polygon is an element in 2D with an arbitrary (but fixed) number of first-order (EDGE2) sides...
Definition: face_c0polygon.h:51
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::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::MeshBase::mesh_dimension
unsigned int mesh_dimension() const
Definition: mesh_base.C:430
libMesh::MeshBase::point
virtual const Point & point(const dof_id_type i) 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::n_elem
virtual dof_id_type n_elem() const =0
libMesh::ParallelObject::processor_id
processor_id_type processor_id() const
Definition: parallel_object.h:114
libMesh::unique_id_type
uint8_t unique_id_type
Definition: id_types.h:86
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::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

◆ change_boundary_id()

void libMesh::MeshTools::Modification::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 at line 2093 of file mesh_modification.C.

References libMesh::MeshBase::get_boundary_info(), mesh, and libMesh::BoundaryInfo::renumber_id().

Referenced by MeshStitchTest::renameAndShift().

2096 {
2097  // This is just a shim around the member implementation, now
2098  mesh.get_boundary_info().renumber_id(old_id, new_id);
2099 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
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

◆ change_subdomain_id()

void libMesh::MeshTools::Modification::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 at line 2103 of file mesh_modification.C.

References libMesh::MeshBase::element_stored_range(), mesh, libMesh::Threads::parallel_for(), libMesh::Elem::subdomain_id(), and libMesh::MeshBase::unset_has_cached_elem_data().

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 }
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::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::StoredRange
The StoredRange class defines a contiguous, divisible set of objects.
Definition: stored_range.h:54
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::element_stored_range
const ElemRange & element_stored_range()
Definition: mesh_base.C:1913

◆ distort()

void libMesh::MeshTools::Modification::distort ( MeshBase &  mesh,
const Real  factor,
const bool  perturb_boundary = false 
)

Randomly perturb the nodal locations.

This function will move each node factor fraction of its minimum neighboring node separation distance. Nodes on the boundary are not moved by default, however they may be by setting the flag perturb_boundary true.

Definition at line 146 of file mesh_modification.C.

References libMesh::MeshBase::clear_point_locator(), libMesh::MeshTools::find_boundary_nodes(), libMesh::Elem::hmin(), libMesh::libmesh_assert(), libMesh::make_range(), libMesh::MeshBase::max_node_id(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::Elem::node_ref_range(), libMesh::MeshBase::query_node_ptr(), and libMesh::TypeVector< T >::unit().

Referenced by main(), DistortTest::perturb_and_check(), VolumeTest::test_true_centroid_and_volume(), and VolumeTest::testQuad4TrueCentroid().

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 }
libMesh::Node
A Node is like a Point, but with more information.
Definition: node.h:52
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
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::TypeVector::unit
TypeVector< T > unit() const
Definition: type_vector.h:1141
libMesh::MeshBase::query_node_ptr
virtual const Node * query_node_ptr(const dof_id_type i) const =0
libMesh::MeshBase::clear_point_locator
void clear_point_locator()
Releases the current PointLocator object.
Definition: mesh_base.C:1859
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
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::MeshBase::mesh_dimension
unsigned int mesh_dimension() const
Definition: mesh_base.C:430
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::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0

◆ flatten()

void libMesh::MeshTools::Modification::flatten ( MeshBase &  mesh)

Removes all the refinement tree structure of Mesh, leaving only the highest-level (most-refined) elements.

This is useful when you want to write out a uniformly-refined grid to be treated later as an initial mesh.

Note
Many functions in LibMesh assume a conforming (with no hanging nodes) grid exists at some level, so you probably only want to do this on meshes which have been uniformly refined.

Definition at line 1967 of file mesh_modification.C.

References libMesh::MeshBase::add_elem(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::build(), libMesh::MeshBase::delete_elem(), libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::get_extra_integer(), libMesh::DofObject::id(), libMesh::index_range(), libMesh::BoundaryInfo::invalid_id, libMesh::MeshBase::is_prepared(), libMesh::MeshBase::is_replicated(), libMesh::libmesh_assert(), libMesh::libmesh_ignore(), libMesh::Elem::mapping_data(), libMesh::Elem::mapping_type(), mesh, libMesh::MeshBase::n_active_elem(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::DofObject::n_extra_integers(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_index_range(), libMesh::Elem::node_ptr(), libMesh::MeshBase::prepare_for_use(), libMesh::DofObject::processor_id(), libMesh::remote_elem, libMesh::Elem::side_index_range(), libMesh::Elem::subdomain_id(), libMesh::Elem::type(), and libMesh::DofObject::unique_id().

Referenced by main().

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 }
libMesh::BoundaryInfo::n_boundary_conds
std::size_t n_boundary_conds() const
Definition: boundary_info.C:2594
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::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::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::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::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:170
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:535
libMesh::MeshBase::delete_elem
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
libMesh::DofObject::id
dof_id_type id() const
Definition: dof_object.h:819
libMesh::MeshBase::add_elem
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
libMesh::libmesh_assert
libmesh_assert(ctx)
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::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::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

◆ orient_elements()

void libMesh::MeshTools::Modification::orient_elements ( MeshBase &  mesh)

Redo the nodal ordering of each element as necessary to give the element Jacobian a positive orientation.

This function does not currently handle meshes with any element refinement.

Definition at line 271 of file mesh_modification.C.

References libMesh::MeshBase::get_boundary_info(), mesh, libMesh::MeshTools::n_levels(), and libMesh::Elem::orient().

Referenced by ElemTest< elem_type >::test_orient_elements().

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 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
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::MeshTools::n_levels
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:826
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57

◆ permute_elements()

void libMesh::MeshTools::Modification::permute_elements ( MeshBase &  mesh)

Randomly permute the nodal ordering of each element (without twisting the element mapping).

This is useful for regression testing with a variety of element orientations.

This function does not currently handle meshes with any element refinement.

This function does not currently permute BoundaryInfo data associated with element sides, which will likely be scrambled by the permutation.

Definition at line 231 of file mesh_modification.C.

References libMesh::make_range(), libMesh::MeshBase::max_elem_id(), mesh, libMesh::MeshTools::n_levels(), libMesh::Elem::n_permutations(), libMesh::Elem::permute(), and libMesh::MeshBase::query_elem_ptr().

Referenced by main(), and FETestBase< order, family, elem_type, N_x >::setUp().

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 }
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::MeshTools::n_levels
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:826
libMesh::MeshBase::max_elem_id
virtual dof_id_type max_elem_id() const =0
libMesh::MeshBase::query_elem_ptr
virtual const Elem * query_elem_ptr(const dof_id_type i) const =0
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::Elem::n_permutations
virtual unsigned int n_permutations() const =0
Returns the number of independent permutations of element nodes - e.g.

◆ redistribute()

void libMesh::MeshTools::Modification::redistribute ( MeshBase &  mesh,
const FunctionBase< Real > &  mapfunc 
)

Deterministically perturb the nodal locations.

This function will move each node from it's current x/y/z coordinates to a new x/y/z coordinate given by the first LIBMESH_DIM components of the specified function mapfunc

Nodes on the boundary are also moved.

Currently, non-vertex nodes are moved in the same way as vertex nodes, according to (newx,newy,newz) = mapfunc(x,y,z). This behavior is often suboptimal for higher order geometries and may be subject to change in future libMesh versions.

Definition at line 289 of file mesh_modification.C.

References libMesh::MeshBase::clear_point_locator(), libMesh::FunctionBase< Output >::clone(), libMesh::libmesh_assert(), mesh, libMesh::MeshBase::n_elem(), and libMesh::MeshBase::n_nodes().

Referenced by libMesh::MeshTools::Generation::build_cube(), FETestBase< order, family, elem_type, N_x >::setUp(), MeshSmootherTest::testLaplaceSmoother(), MeshSmootherTest::testVariationalSmoother(), and MeshSmootherTest::testVariationalSmootherRegression().

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 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::MeshBase::clear_point_locator
void clear_point_locator()
Releases the current PointLocator object.
Definition: mesh_base.C:1859
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::FunctionBase::clone
virtual std::unique_ptr< FunctionBase< Output > > clone() const =0
libMesh::DenseVector
Defines a dense vector for use in Finite Element-type computations.
Definition: dof_map.h:74
libMesh::MeshBase::n_elem
virtual dof_id_type n_elem() const =0
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0

◆ rotate()

RealTensorValue libMesh::MeshTools::Modification::rotate ( MeshBase &  mesh,
const Real  phi,
const Real  theta = 0.,
const Real  psi = 0. 
)

Rotates the mesh in the xy plane.

The rotation is counter-clock-wise (mathematical definition). The angle is in degrees (360 make a full circle) Rotates the mesh in 3D space. Here the standard Euler angles are adopted (http://mathworld.wolfram.com/EulerAngles.html) The angles are in degrees (360 make a full circle)

Returns
the 3x3 rotation matrix implied by (phi, theta, psi)

Definition at line 361 of file mesh_modification.C.

References libMesh::MeshBase::clear_point_locator(), libMesh::TensorValue< T >::intrinsic_rotation_matrix(), libMesh::libmesh_ignore(), mesh, and libMesh::MeshBase::set_spatial_dimension().

Referenced by main(), libMesh::Prism6::permute(), libMesh::Hex8::permute(), libMesh::Hex20::permute(), libMesh::Tet10::permute(), libMesh::Tet14::permute(), libMesh::Prism15::permute(), libMesh::Tet4::permute(), libMesh::Hex27::permute(), libMesh::Prism18::permute(), libMesh::Prism20::permute(), libMesh::Prism21::permute(), FETestBase< order, family, elem_type, N_x >::setUp(), BBoxTest::test_no_degenerate(), and BBoxTest::test_two_degenerate().

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 }
libMesh::MeshBase::set_spatial_dimension
void set_spatial_dimension(unsigned char d)
Sets the "spatial dimension" of the Mesh.
Definition: mesh_base.C:613
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::RealTensorValue
TensorValue< Real > RealTensorValue
Useful typedefs to allow transparent switching between Real and Complex data types.
Definition: hp_coarsentest.h:48
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:535
libMesh::MeshBase::clear_point_locator
void clear_point_locator()
Releases the current PointLocator object.
Definition: mesh_base.C:1859
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39

◆ scale()

void libMesh::MeshTools::Modification::scale ( MeshBase &  mesh,
const Real  xs,
const Real  ys = 0.,
const Real  zs = 0. 
)

Scales the mesh.

The grid points are scaled in the x direction by xs, in the y direction by ys, etc... If only xs is specified then the scaling is assumed uniform in all directions.

Definition at line 393 of file mesh_modification.C.

References libMesh::MeshBase::clear_point_locator(), mesh, and libMesh::Real.

Referenced by libMesh::PetscVector< libMesh::Number >::add(), main(), libMesh::DenseVector< Output >::operator*=(), libMesh::DenseMatrix< Real >::operator*=(), libMesh::PetscMatrix< T >::scale(), and libMesh::SparseMatrix< ValOut >::scale().

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 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::MeshBase::clear_point_locator
void clear_point_locator()
Releases the current PointLocator object.
Definition: mesh_base.C:1859
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144

◆ smooth()

void libMesh::MeshTools::Modification::smooth ( MeshBase &  mesh,
unsigned int  n_iterations,
Real  power 
)

Smooth the mesh with a simple Laplace smoothing algorithm.

The mesh is smoothed n_iterations times. If the parameter power is 0, each node is moved to the average position of the neighboring connected nodes. If power > 0, the node positions are weighted by their distance. The positions of higher order nodes, and nodes living in refined elements, are calculated from the vertex positions of their parent nodes. Only works in 2D.

Author
Martin Luthi (luthi.nosp@m.@gi..nosp@m.alask.nosp@m.a.ed.nosp@m.u)
Date
2005

This implementation assumes every element "side" has only 2 nodes.

Definition at line 1807 of file mesh_modification.C.

References libMesh::TypeVector< T >::add(), libMesh::TypeVector< T >::add_scaled(), libMesh::as_range(), libMesh::MeshBase::clear_point_locator(), libMesh::Elem::embedding_matrix(), libMesh::MeshTools::find_boundary_nodes(), libMesh::DofObject::id(), libMesh::make_range(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::MeshTools::n_levels(), libMesh::MeshBase::n_nodes(), libMesh::Elem::n_nodes(), libMesh::Elem::n_second_order_adjacent_vertices(), libMesh::Elem::n_vertices(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_index_range(), libMesh::Elem::node_ptr(), libMesh::Elem::node_ref(), libMesh::MeshBase::node_ref(), libMesh::TypeVector< T >::norm(), libMesh::Elem::parent(), libMesh::Elem::point(), libMesh::Utility::pow(), libMesh::Real, libMesh::Elem::second_order_adjacent_vertex(), libMesh::Elem::side_index_range(), libMesh::MeshTools::weight(), and libMesh::Elem::which_child_am_i().

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 }
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::Node
A Node is like a Point, but with more information.
Definition: node.h:52
libMesh::TypeVector::norm
auto norm() const
Definition: type_vector.h:908
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::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::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::Utility::pow
T pow(const T &x)
Definition: utility.h:296
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::MeshTools::weight
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:444
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::n_levels
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:826
libMesh::MeshBase::clear_point_locator
void clear_point_locator()
Releases the current PointLocator object.
Definition: mesh_base.C:1859
libMesh::ElemSideBuilder
Helper for building element sides that minimizes the construction of new elements.
Definition: elem_side_builder.h:50
libMesh::Elem::which_child_am_i
unsigned int which_child_am_i(const Elem *e) const
Definition: elem.h:3206
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
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::MeshBase::mesh_dimension
unsigned int mesh_dimension() const
Definition: mesh_base.C:430
libMesh::Elem::node_index_range
IntRange< unsigned short > node_index_range() const
Definition: elem.h:2697
libMesh::MeshBase::node_ref
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:735
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:2459
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67

◆ translate()

void libMesh::MeshTools::Modification::translate ( MeshBase &  mesh,
const Real  xt = 0.,
const Real  yt = 0.,
const Real  zt = 0. 
)

Translates the mesh.

The grid points are translated in the x direction by xt, in the y direction by yt, etc...

Definition at line 322 of file mesh_modification.C.

References libMesh::MeshBase::clear_point_locator(), and mesh.

Referenced by MeshStitchTest::testMeshStitchElemsets(), and MeshTriangulationTest::testTriangulatorRoundHole().

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 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::MeshBase::clear_point_locator
void clear_point_locator()
Releases the current PointLocator object.
Definition: mesh_base.C:1859
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
Generated on Wed Jun 24 2026 21:22:00 for libMesh by   doxygen 1.8.14