Tools for Mesh modification. More...

Functions
void	distort (MeshBase &mesh, const Real factor, const bool perturb_boundary=false)
	Randomly perturb the nodal locations.

void	permute_elements (MeshBase &mesh)
	Randomly permute the nodal ordering of each element (without twisting the element mapping).

void	orient_elements (MeshBase &mesh)
	Redo the nodal ordering of each element as necessary to give the element Jacobian a positive orientation.

void	redistribute (MeshBase &mesh, const FunctionBase< Real > &mapfunc)
	Deterministically perturb the nodal locations.

void	translate (MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.)
	Translates the mesh.

RealTensorValue	rotate (MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
	Rotates the mesh in 3D space.

void	scale (MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
	Scales the mesh.

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

void	all_rbb (MeshBase &mesh)
	Converts all element geometric mappings from the default Lagrange to the more flexible Rational-Bezier-Bernstein.

void	smooth (MeshBase &, unsigned int, Real)
	Smooth the mesh with a simple Laplace smoothing algorithm.

void	interpolate_surface (MeshBase &mesh, const Surface &surface, std::set< std::size_t > ids={}, bool use_boundary_nodes=true)
	Move nodes in `mesh` to their closest points on the specified `surface`.

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

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.

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.

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 1260 of file mesh_modification.C.

{
  LOG_SCOPE("all_rbb()", "MeshTools::Modification");
 
  // By default, use 1.0 as the weight on every RATIONAL_BERNSTEIN
  // mapped node
  const Real default_weight = 1.0;
 
  const auto weight_index =
    (mesh.add_node_datum<Real>("rational_weight", true,
                               &default_weight));
 
  mesh.set_default_mapping_type(RATIONAL_BERNSTEIN_MAP);
  mesh.set_default_mapping_data(weight_index);
 
  // Out of loop to reduce heap allocations
  std::unique_ptr<Elem> edge_ptr, face_ptr;
 
  for (auto & elem : mesh.element_ptr_range())
    {
      if (elem->level())
        libmesh_not_implemented_msg
          ("all_rbb() currently only supports flat meshes with no refinement levels");
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
      if (elem->infinite())
        libmesh_not_implemented_msg
          ("all_rbb() currently only supports finite geometric elements");
#endif
 
      elem->set_mapping_type(RATIONAL_BERNSTEIN_MAP);
      elem->set_mapping_data(weight_index);
 
      // Nothing to do unless we have curves to correct
      if (elem->default_order() == FIRST)
        continue;
 
      // Modify the center node of an "edge" - possibly an actual edge
      // element's node, possibly a center node between points on a
      // face's or cell's edge - for RBB interpolation.  This should fit
      // a circular curve exactly in cases where the original nodes are
      // equispaced and the outer nodes' weights are equal, and should
      // be a good fit otherwise.
      //
      // We want to use this to interpolate "internal" conceptual
      // edges of a Hex27 too, so we'll handle the cases where w0 and
      // w1 aren't 1, as well as the cases where the Nodes n0 and n1
      // are already control points which don't match their
      // corresponding physical points.
      auto make_edge_rbb = [default_weight, weight_index]
        (const Node & n0, const Node & n1, Node & n_center,
         const Point & p0, const Point & p1)
      {
        // Skip edges we've already modified; the center node for
        // these is no longer at the curve point we wish to
        // interpolate, it should already be at the control point that
        // accomplishes the interpolation.
        const Real old_weight = n_center.get_extra_datum<Real>(weight_index);
        if (old_weight != default_weight)
          return;
 
        Point & p2 = n_center;
 
        const Real w0 = n0.get_extra_datum<Real>(weight_index);
        const Real w1 = n1.get_extra_datum<Real>(weight_index);
 
        const Point e02 = p2-p0,
                    e21 = p1-p2;
        const Real chord_02_len_sq = e02.norm_sq(),
                   chord_21_len_sq = e21.norm_sq();
 
        // First find the cosine of phi, the angle between our two
        // subchords (turning from the direction of one to the
        // direction of the other; this is the supplementary angle to
        // the angle at the midpoint).  This is the same as half of
        // the angle of our circular arc, which nicely enough is also
        // the angle we take cos and sec of in NURBS formulae
        const Real cos_phi = (e02*e21)/std::sqrt(chord_02_len_sq*chord_21_len_sq);
 
        // There's a way to do really large arcs using negative
        // weights, but we're going to get lousy approximation quality
        // from isoparametric elements if we go too low, as well as
        // bad numerics here, so let's just disallow it.
        if (cos_phi < 0.5)
          libmesh_not_implemented_msg
            ("all_rbb() is not recommended for extremely sharp curves on one edge");
 
        const Real w_center = cos_phi*std::sqrt(w0*w1);
 
        n_center.set_extra_datum<Real>(weight_index, w_center);
 
        // Now let's get the control point location.  This comes from
        // a lot of back-and-forth with Gemini, but fortunately I'm
        // rewriting it after I've already added unit tests that
        // should scream if it's badly wrong.
        const Real w_mid = w0/4 + w1/4 + w_center/2;
        n_center *= 2*w_mid;
        n_center -= (w0 * p0 + w1 * p1)/2;
        n_center /= w_center;
      };
 
      auto make_face_rbb = [weight_index] (Elem & face)
      {
        // Prisms and pyramids may need to skip some faces while
        // adjusting others
        if (face.type() == TRI6)
          return;
 
        if (face.type() != QUAD9)
          libmesh_not_implemented_msg
            ("all_rbb() currently only supports mid-face nodes on Quad9 faces");
 
        // We only use [4,8) but matching indices is nice and stack is
        // cheap.
        Real w[9];
 
        for (unsigned int i : make_range(4u, 8u))
          w[i] = face.node_ref(i).get_extra_datum<Real>(weight_index);
 
        // We can't currently handle arbitrary vertex weights
#ifndef NDEBUG
        for (unsigned int i : make_range(4u))
          libmesh_assert_equal_to
            (face.node_ref(i).get_extra_datum<Real>(weight_index), 1);
#endif
 
        // For the mid-face point, if we want to exactly match
        // any cylinders and cones and spheres, we're actually already
        // entirely constrained by the other points.
        //
        // This formula gives the minimum-energy Steiner surface based
        // on the outer 8 points.
        //
        // That's an isogeometric representation of a cylinder aligned
        // to either axis, or of a sphere where the quad edges are on
        // latitude/longitude lines, or of a cone where two edges are
        // segments of cone generating lines and the other two are
        // arcs perpendicular to the axis.
        //
        // It's not perfectly isogeometric for the spheres we generate
        // (where the quad edges are all great circles), but it should
        // still converge asymptotically faster than non-rational
        // quadratic Lagrange.
        const Point xi_avg = (face.point(7) + face.point(5))/2;
        const Point eta_avg = (face.point(4) + face.point(6))/2;
        const Point vertex_avg = (face.point(0) + face.point(1) +
                                  face.point(2) + face.point(3))/4;
 
        const Real w_xi  = (w[7] + w[5])/2;
        const Real w_eta = (w[4] + w[6])/2;
        const Real w_mid = w_xi * w_eta;
 
        Node & midnode = face.node_ref(8);
        midnode.set_extra_datum<Real>(weight_index, w_mid);
        midnode = ((1+w_mid)/(w_xi+w_eta) * (w_xi*xi_avg + w_eta*eta_avg) - vertex_avg)/w_mid;
      };
 
      // If we're on a Hex27, our formula for the mid-volume node
      // relies on the locations of the mid-face points.  We could
      // re-calculate those later but let's just save them now.
      Point midfacepts[6];
      if (elem->type() == HEX27)
        for (auto i : make_range(6))
          midfacepts[i] = elem->point(20+i);
 
      // Check each edge for a curve, and adjust it if needed.
      for (auto e : elem->edge_index_range())
        {
          elem->build_edge_ptr(edge_ptr, e);
 
          // We should add EDGE4 once we have QUAD16/TRI10/HEX64 to
          // use it
          if (edge_ptr->type() != EDGE3)
            libmesh_not_implemented_msg
              ("all_rbb() currently only supports meshes with 2- and/or 3-node edges");
 
          make_edge_rbb(edge_ptr->node_ref(0), edge_ptr->node_ref(1),
                        edge_ptr->node_ref(2),
                        edge_ptr->node_ref(0), edge_ptr->node_ref(1));
 
        }
 
      // If we're in 3D, we may have face nodes that also need to be
      // adjusted to replace an interpolated curve with a spline
      // curve.  We know what to do with a quad face, but we'll have
      // to scream and die if we see a Tri7 face node.
      bool check_face_points = (elem->dim() > 2) &&
        (elem->n_nodes() > elem->n_edges() + elem->n_vertices());
 
      if (check_face_points)
        for (auto f : elem->side_index_range())
          {
            // Prisms and pyramids may need to skip some faces while
            // adjusting others
            if (elem->side_type(f) == TRI6)
              continue;
 
            elem->build_side_ptr(face_ptr, f);
 
            make_face_rbb(*face_ptr);
          }
 
      bool check_interior_points =
        elem->n_nodes() > elem->n_edges() + elem->n_vertices() + elem->n_faces();
 
      if (check_interior_points)
        {
          if (elem->type() == EDGE3)
            {
              make_edge_rbb(elem->node_ref(0), elem->node_ref(1),
                            elem->node_ref(2),
                            elem->node_ref(0), elem->node_ref(1));
            }
          else if (elem->dim() == 2)
            {
              make_face_rbb(*elem);
            }
          else if (elem->type() == HEX27)
            {
              // We still have the midnode left to go.  We want
              // something here that will preserve the tensor product
              // structure for 2.5D extrusions of IGA faces, but also
              // be at least near to the minimum-energy control point
              // and weight for general cases.  We'll treat opposing
              // mid-face nodes as the endpoints of a (more general
              // than our edges, since they might have non-1 weights)
              // Edge3, and see what we'd need on the midnode to
              // interpolate the center point with them.  If we've got
              // something isogeometric like an extrusion then our
              // results should agree; for a quick-but-good output in
              // general we'll take an average.
              const int opposite_sides[3][2] = {{0,5}, {1,3}, {2,4}};
 
              Node & midnode = elem->node_ref(26);
              const Point original_midpoint = midnode;
 
              // Averaging in projective space
              Point sum_weighted_point = 0;
              Real sum_weight = 0;
 
              for (int i : make_range(3))
                {
                  Node & n0 = elem->node_ref(20+opposite_sides[i][0]);
                  Node & n1 = elem->node_ref(20+opposite_sides[i][1]);
 
                  make_edge_rbb(n0, n1, midnode,
                                midfacepts[opposite_sides[i][0]],
                                midfacepts[opposite_sides[i][1]]);
 
                  const Real midweight =
                    midnode.get_extra_datum<Real>(weight_index);
                  sum_weight += midweight;
                  sum_weighted_point += midweight * midnode;
 
                  // Reset for next run
                  midnode = original_midpoint;
                  midnode.set_extra_datum<Real>(weight_index,
                                                default_weight);
 
                }
 
              const Real midweight = sum_weight/3;
              midnode.set_extra_datum<Real>(weight_index,
                                            midweight);
 
              midnode = sum_weighted_point / 3 / midweight;
            }
          else
            libmesh_not_implemented_msg
              ("all_rbb() doesn't yet support " << elem->type());
        }
    }
}

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::TypeVector< T >::norm_sq(), 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(), AllRBBTest::test_disk(), and AllRBBTest::test_sphere().

◆ 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 449 of file mesh_modification.C.

{
  LOG_SCOPE("all_tri()", "MeshTools::Modification");
 
  if (!mesh.is_replicated() && !mesh.is_prepared())
    mesh.prepare_for_use();
 
  // The number of elements in the original mesh before any additions
  // or deletions.
  const dof_id_type n_orig_elem = mesh.n_elem();
  const dof_id_type max_orig_id = mesh.max_elem_id();
 
  // We store pointers to the newly created elements in a vector
  // until they are ready to be added to the mesh.  This is because
  // adding new elements on the fly can cause reallocation and invalidation
  // of existing mesh element_iterators.
  std::vector<std::unique_ptr<Elem>> new_elements;
 
  unsigned int max_subelems = 1;  // in 1D nothing needs to change
  if (mesh.mesh_dimension() == 2) // in 2D quads can split into 2 tris
    max_subelems = 2;
  if (mesh.mesh_dimension() == 3) // in 3D hexes can split into 6 tets
    max_subelems = 6;
 
  // 2D polygons and 3D polyhedra can be split into an arbitrary
  // number of triangles/tetrahedra depending on their topology, so we
  // have to scan the mesh to find the largest split we will need.
  for (const Elem * elem : mesh.element_ptr_range())
    {
      if (const Polygon * poly = dynamic_cast<const Polygon *>(elem))
        max_subelems = std::max(max_subelems, poly->n_subtriangles());
      else if (const Polyhedron * polyhedron = dynamic_cast<const Polyhedron *>(elem))
        max_subelems = std::max(max_subelems, polyhedron->n_subelements());
    }
  mesh.comm().max(max_subelems);
 
  new_elements.reserve (max_subelems*n_orig_elem);
 
  // If the original mesh has *side* boundary data, we carry that over
  // to the new mesh with triangular elements.  We currently only
  // support bringing over side-based BCs to the all-tri mesh, but
  // that could probably be extended to node and edge-based BCs as
  // well.
  const bool mesh_has_boundary_data = (mesh.get_boundary_info().n_boundary_conds() > 0);
 
  // Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
  std::vector<Elem *> new_bndry_elements;
  std::vector<unsigned short int> new_bndry_sides;
  std::vector<boundary_id_type> new_bndry_ids;
 
  // We may need to add new points if we run into a 1.5th order
  // element; if we do that on a DistributedMesh in a ghost element then
  // we will need to fix their ids / unique_ids
  bool added_new_ghost_point = false;
 
  // Iterate over the elements, splitting:
  // QUADs into pairs of conforming triangles
  // PYRAMIDs into pairs of conforming tets,
  // PRISMs into triplets of conforming tets, and
  // HEXs into quintets or sextets of conforming tets.
  // We split on the shortest diagonal to give us better
  // triangle quality in 2D, and we split based on node ids
  // to guarantee consistency in 3D.
  // C0POLYGONs into their sub-triangles
  // C0POLYHEDRA into their sub-elements (currently only tets)
 
  // FIXME: This algorithm does not work on refined grids!
  {
#ifdef LIBMESH_ENABLE_UNIQUE_ID
    unique_id_type max_unique_id = mesh.parallel_max_unique_id();
#endif
 
    // For avoiding extraneous allocation when building side elements
    std::unique_ptr<const Elem> elem_side, subside_elem;
 
    for (auto & elem : mesh.element_ptr_range())
      {
        const ElemType etype = elem->type();
 
        // all_tri currently only works on coarse meshes
        if (elem->parent())
          libmesh_not_implemented_msg("Cannot convert a refined element into simplices\n");
 
        // The new elements we will split the original into. Reserving
        // for the maximum number of sub-elements created for each element
        std::vector<std::unique_ptr<Elem>> subelem(max_subelems);
 
        auto set_nodes = [&elem, &subelem]
          (const std::initializer_list<std::initializer_list<int>> & node_ids) {
          int i=0;
          for (auto row : node_ids)
            {
              int j=0;
              Elem * sub = subelem[i++].get();
              libmesh_assert(sub);
              for (auto node_id : row)
                sub->set_node(j++, elem->node_ptr(node_id));
            }
        };
 
        switch (etype)
          {
          case QUAD4:
            {
              subelem[0] = Elem::build(TRI3);
              subelem[1] = Elem::build(TRI3);
 
              // Check for possible edge swap
              if ((elem->point(0) - elem->point(2)).norm() <
                  (elem->point(1) - elem->point(3)).norm())
                set_nodes({{0,1,2},{0,2,3}});
              else
                set_nodes({{0,1,3},{1,2,3}});
 
              break;
            }
 
          case QUAD8:
            {
              if (elem->processor_id() != mesh.processor_id())
                added_new_ghost_point = true;
 
              subelem[0] = Elem::build(TRI6);
              subelem[1] = Elem::build(TRI6);
 
              // Add a new node at the center (vertex average) of the element.
              Node * new_node = mesh.add_point((mesh.point(elem->node_id(0)) +
                                                mesh.point(elem->node_id(1)) +
                                                mesh.point(elem->node_id(2)) +
                                                mesh.point(elem->node_id(3)))/4,
                                               DofObject::invalid_id,
                                               elem->processor_id());
 
              // Check for possible edge swap
              if ((elem->point(0) - elem->point(2)).norm() <
                  (elem->point(1) - elem->point(3)).norm())
                {
                  set_nodes({{0,1,2,4,5},{0,2,3,3,6,7}});
                  subelem[0]->set_node(5, new_node);
                  subelem[1]->set_node(3, new_node);
                }
              else
                {
                  set_nodes({{3,0,1,7,4},{1,2,3,5,6}});
                  subelem[0]->set_node(5, new_node);
                  subelem[1]->set_node(5, new_node);
                }
 
              break;
            }
 
          case QUAD9:
            {
              subelem[0] = Elem::build(TRI6);
              subelem[1] = Elem::build(TRI6);
 
              // Check for possible edge swap
              if ((elem->point(0) - elem->point(2)).norm() <
                  (elem->point(1) - elem->point(3)).norm())
                set_nodes({{0,1,2,4,5,8},{0,2,3,8,6,7}});
              else
                set_nodes({{0,1,3,4,8,7},{1,2,3,5,6,8}});
 
              break;
            }
 
          case HEX8:
            {
              BoundaryInfo & boundary_info = mesh.get_boundary_info();
 
              // Hexes all split into six tetrahedra
              subelem[0] = Elem::build(TET4);
              subelem[1] = Elem::build(TET4);
              subelem[2] = Elem::build(TET4);
              subelem[3] = Elem::build(TET4);
              subelem[4] = Elem::build(TET4);
              subelem[5] = Elem::build(TET4);
 
              // On faces, we choose the node with the highest
              // global id, and we split on the diagonal which
              // includes that node.  This ensures that (even in
              // parallel, even on distributed meshes) the same
              // diagonal split will be chosen for elements on either
              // side of the same quad face.
              const unsigned int highest_n = highest_vertex_on(elem);
 
              // opposing_node[n] is the local node number of the node
              // on the farthest corner of a hex8 from local node n
              static const std::array<unsigned int, 8> opposing_node =
                {6, 7, 4, 5, 2, 3, 0, 1};
 
              static const std::vector<std::vector<unsigned int>> sides_opposing_highest =
                {{2,3,5},{3,4,5},{1,4,5},{1,2,5},{0,2,3},{0,3,4},{0,1,4},{0,1,2}};
              static const std::vector<std::vector<unsigned int>> nodes_neighboring_highest =
                {{1,3,4},{0,2,5},{1,3,6},{0,2,7},{0,5,7},{1,4,6},{2,5,7},{3,4,6}};
 
              // Start by looking in three directions away from the
              // highest-id node.  In each direction there will be two
              // different possibilities for the split depending on
              // how the opposing face nodes are numbered.
              //
              // This is tricky enough that I'm not going to worry
              // about manually keeping tets oriented; we'll just call
              // orient() on each as we go.
 
              unsigned int next_subelem = 0;
              for (auto side : sides_opposing_highest[highest_n])
                {
                  const std::vector<unsigned int> nodes_on_side =
                    elem->nodes_on_side(side);
 
                  auto [dn, dn2] = split_diagonal(elem, nodes_on_side);
 
                  unsigned int split_on_neighbor = false;
                  for (auto n : nodes_neighboring_highest[highest_n])
                    if (dn == n || dn2 == n)
                      {
                        split_on_neighbor = true;
                        break;
                      }
 
                  // Add one or two elements for each opposing side,
                  // depending on whether the diagonal split there
                  // connects to the neighboring diagonal split or
                  // not.
                  if (split_on_neighbor)
                    {
                      subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
                      subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
                      subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
                      for (auto n : nodes_on_side)
                        if (n != dn && n != dn2)
                          {
                            subelem[next_subelem]->set_node(3, elem->node_ptr(n));
                            break;
                          }
                      subelem[next_subelem]->orient(&boundary_info);
                      ++next_subelem;
 
                      subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
                      subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
                      subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
                      for (auto n : reverse(nodes_on_side))
                        if (n != dn && n != dn2)
                          {
                            subelem[next_subelem]->set_node(3, elem->node_ptr(n));
                            break;
                          }
                      subelem[next_subelem]->orient(&boundary_info);
                      ++next_subelem;
                    }
                  else
                    {
                      subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
                      subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
                      subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
                      for (auto n : nodes_on_side)
                        for (auto n2 : nodes_neighboring_highest[highest_n])
                          if (n == n2)
                            {
                              subelem[next_subelem]->set_node(3, elem->node_ptr(n));
                              goto break_both_loops;
                            }
 
                      break_both_loops:
                      subelem[next_subelem]->orient(&boundary_info);
                      ++next_subelem;
                    }
                }
 
              // At this point we've created between 3 and 6 tets.
              // What's left to do depends on how many.
 
              // If we just chopped off three vertices into three
              // tets, then the best way to split this hex would be
              // the symmetric five-split.  Chop off the opposing
              // vertex too, and then the remaining interior is our
              // final tet.
              if (next_subelem == 3)
                {
                  subelem[next_subelem]->set_node(0, elem->node_ptr(opposing_nodes[highest_n][0]));
                  subelem[next_subelem]->set_node(1, elem->node_ptr(opposing_nodes[highest_n][1]));
                  subelem[next_subelem]->set_node(2, elem->node_ptr(opposing_nodes[highest_n][2]));
                  subelem[next_subelem]->set_node(3, elem->node_ptr(opposing_node[highest_n]));
                  subelem[next_subelem]->orient(&boundary_info);
                  ++next_subelem;
 
                  subelem[next_subelem]->set_node(0, elem->node_ptr(opposing_nodes[highest_n][0]));
                  subelem[next_subelem]->set_node(1, elem->node_ptr(opposing_nodes[highest_n][1]));
                  subelem[next_subelem]->set_node(2, elem->node_ptr(opposing_nodes[highest_n][2]));
                  subelem[next_subelem]->set_node(3, elem->node_ptr(highest_n));
                  subelem[next_subelem]->orient(&boundary_info);
                  ++next_subelem;
 
                  // We don't need the 6th tet after all
                  subelem[next_subelem].reset();
                  ++next_subelem;
                }
 
              // If we just chopped off one (or two) vertices into
              // tets, then the remaining gap is best (or only) filled
              // by pairing another tet with each.
              if (next_subelem == 4 ||
                  next_subelem == 5)
                {
                  for (auto side : sides_opposing_highest[highest_n])
                    {
                      const std::vector<unsigned int> nodes_on_side =
                        elem->nodes_on_side(side);
 
                      auto [dn, dn2] = split_diagonal(elem, nodes_on_side);
 
                      unsigned int split_on_neighbor = false;
                      for (auto n : nodes_neighboring_highest[highest_n])
                        if (dn == n || dn2 == n)
                          {
                            split_on_neighbor = true;
                            break;
                          }
 
                      // The two !split_on_neighbor sides are where we
                      // need the two remaining tets
                      if (!split_on_neighbor)
                        {
                          subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
                          subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
                          subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
                          subelem[next_subelem]->set_node(3, elem->node_ptr(opposing_node[highest_n]));
                          subelem[next_subelem]->orient(&boundary_info);
                          ++next_subelem;
                        }
                    }
                }
 
              // Whether we got there by creating six tets from the
              // first for loop or by patching up the split afterward,
              // we should have considered six tets (possibly
              // including one deleted one...) at this point.
              libmesh_assert(next_subelem == 6);
 
              break;
            }
 
          case PRISM6:
            {
              // Prisms all split into three tetrahedra
              subelem[0] = Elem::build(TET4);
              subelem[1] = Elem::build(TET4);
              subelem[2] = Elem::build(TET4);
 
              // Triangular faces are not split.
 
              // On quad faces, we choose the node with the highest
              // global id, and we split on the diagonal which
              // includes that node.  This ensures that (even in
              // parallel, even on distributed meshes) the same
              // diagonal split will be chosen for elements on either
              // side of the same quad face.  It also ensures that we
              // always have a mix of "clockwise" and
              // "counterclockwise" split faces (two of one and one
              // of the other on each prism; this is useful since the
              // alternative all-clockwise or all-counterclockwise
              // face splittings can't be turned into tets without
              // adding more nodes
 
              // Split on 0-4 diagonal
              if (split_first_diagonal(elem, 0,4, 1,3))
                {
                  // Split on 0-5 diagonal
                  if (split_first_diagonal(elem, 0,5, 2,3))
                    {
                      // Split on 1-5 diagonal
                      if (split_first_diagonal(elem, 1,5, 2,4))
                        set_nodes({{0,4,5,3},{0,4,1,5},{0,1,2,5}});
                      else // Split on 2-4 diagonal
                        {
                          libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
                          set_nodes({{0,4,5,3},{0,4,2,5},{0,1,2,4}});
                        }
                    }
                  else // Split on 2-3 diagonal
                    {
                      libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
 
                      // 0-4 and 2-3 split implies 2-4 split
                      libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
 
                      set_nodes({{0,4,2,3},{3,4,2,5},{0,1,2,4}});
                    }
                }
              else // Split on 1-3 diagonal
                {
                  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
 
                  // Split on 0-5 diagonal
                  if (split_first_diagonal(elem, 0,5, 2,3))
                    {
                      // 1-3 and 0-5 split implies 1-5 split
                      libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
 
                      set_nodes({{1,3,4,5},{1,0,3,5},{0,1,2,5}});
                    }
                  else // Split on 2-3 diagonal
                    {
                      libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
 
                      // Split on 1-5 diagonal
                      if (split_first_diagonal(elem, 1,5, 2,4))
                        set_nodes({{0,1,2,3},{3,1,2,5},{1,3,4,5}});
                      else // Split on 2-4 diagonal
                        {
                          libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
                          set_nodes({{0,1,2,3},{2,3,4,5},{3,1,2,4}});
                        }
                    }
                }
 
              break;
            }
 
          case PRISM20:
          case PRISM21:
            libmesh_experimental(); // We should upgrade this to TET14...
            libmesh_fallthrough();
          case PRISM18:
            {
              subelem[0] = Elem::build(TET10);
              subelem[1] = Elem::build(TET10);
              subelem[2] = Elem::build(TET10);
 
              // Split on 0-4 diagonal
              if (split_first_diagonal(elem, 0,4, 1,3))
                {
                  // Split on 0-5 diagonal
                  if (split_first_diagonal(elem, 0,5, 2,3))
                    {
                      // Split on 1-5 diagonal
                      if (split_first_diagonal(elem, 1,5, 2,4))
                        set_nodes({{0,4,5,3,15,13,17,9,12,14},
                                   {0,4,1,5,15,10,6,17,13,16},
                                   {0,1,2,5,6,7,8,17,16,11}});
                      else // Split on 2-4 diagonal
                        {
                          libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
 
                          set_nodes({{0,4,5,3,15,13,17,9,12,14},
                                     {0,4,2,5,15,16,8,17,13,11},
                                     {0,1,2,4,6,7,8,15,10,16}});
                        }
                    }
                  else // Split on 2-3 diagonal
                    {
                      libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
 
                      // 0-4 and 2-3 split implies 2-4 split
                      libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
 
                      set_nodes({{0,4,2,3,15,16,8,9,12,17},
                                 {3,4,2,5,12,16,17,14,13,11},
                                 {0,1,2,4,6,7,8,15,10,16}});
                    }
                }
              else // Split on 1-3 diagonal
                {
                  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
 
                  // Split on 0-5 diagonal
                  if (split_first_diagonal(elem, 0,5, 2,3))
                    {
                      // 1-3 and 0-5 split implies 1-5 split
                      libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
 
                      set_nodes({{1,3,4,5,15,12,10,16,14,13},
                                 {1,0,3,5,6,9,15,16,17,14},
                                 {0,1,2,5,6,7,8,17,16,11}});
                    }
                  else // Split on 2-3 diagonal
                    {
                      libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
 
                      // Split on 1-5 diagonal
                      if (split_first_diagonal(elem, 1,5, 2,4))
                        set_nodes({{0,1,2,3,6,7,8,9,15,17},
                                   {3,1,2,5,15,7,17,14,16,11},
                                   {1,3,4,5,15,12,10,16,14,13}});
                      else // Split on 2-4 diagonal
                        {
                          libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
 
                          set_nodes({{0,1,2,3,6,7,8,9,15,17},
                                     {2,3,4,5,17,12,16,11,14,13},
                                     {3,1,2,4,15,7,17,12,10,16}});
                        }
                    }
                }
 
              break;
            }
 
          case PYRAMID5:
            {
              // Pyramids all split into two tetrahedra
              subelem[0] = Elem::build(TET4);
              subelem[1] = Elem::build(TET4);
 
              // Choose how to split the quad face in a way that will
              // be consistent from possibly-different-type elements
              // splitting from the other side
              //
              // Split on 0-2 diagonal
              if (split_first_diagonal(elem, 0,2, 1,3))
                set_nodes({{0,1,2,4},{0,2,3,4}});
              // Split on 1-3 diagonal
              else
                {
                  libmesh_assert (split_first_diagonal(elem, 1,3, 0,2));
                  set_nodes({{0,1,3,4},{1,2,3,4}});
                }
 
              break;
            }
 
          case PYRAMID14:
            {
              // Pyramids all split into two tetrahedra
              subelem[0] = Elem::build(TET10);
              subelem[1] = Elem::build(TET10);
 
              // Choose how to split the quad face in a way that will
              // be consistent from possibly-different-type elements
              // splitting from the other side
              //
              // Split on 0-2 diagonal
              if (split_first_diagonal(elem, 0,2, 1,3))
                set_nodes({{0,1,2,4,5,6,13,9,10,11},
                           {0,2,3,4,13,7,8,9,11,12}});
              // Split on 1-3 diagonal
              else
                {
                  libmesh_assert (split_first_diagonal(elem, 1,3, 0,2));
                  set_nodes({{0,1,3,4,5,13,8,9,10,12},
                             {1,2,3,4,6,7,13,10,11,12}});
                }
 
              break;
            }
 
          case C0POLYGON:
            {
              // Split a C0Polygon into the triangles defined by its
              // current triangulation.  This relies on the user having
              // a valid triangulation (the constructor sets a default
              // one, and the user can refresh it via retriangulate()
              // after moving nodes).
              const C0Polygon * polygon = cast_ptr<const C0Polygon *>(elem);
              const unsigned int n_subtri = polygon->n_subtriangles();
              for (unsigned int t = 0; t != n_subtri; ++t)
                {
                  const std::array<int, 3> tri = polygon->subtriangle(t);
                  if (tri[0] < 0 || tri[1] < 0 || tri[2] < 0)
                    libmesh_not_implemented_msg
                      ("Cannot convert a C0Polygon whose triangulation\n"
                       "introduces special (non-vertex) points");
                  subelem[t] = Elem::build(TRI3);
                  subelem[t]->set_node(0, elem->node_ptr(tri[0]));
                  subelem[t]->set_node(1, elem->node_ptr(tri[1]));
                  subelem[t]->set_node(2, elem->node_ptr(tri[2]));
                }
 
              break;
            }
 
          case C0POLYHEDRON:
            {
              // Split a C0Polyhedron into the tetrahedra defined by its
              // current tetrahedralization.  If the polyhedron required
              // a mid-element node, the user is expected to have added
              // that node to the mesh during construction; we just
              // reference it via the polyhedron's node pointers.
              const C0Polyhedron * polyhedron =
                cast_ptr<const C0Polyhedron *>(elem);
              const unsigned int n_sub = polyhedron->n_subelements();
              for (unsigned int t = 0; t != n_sub; ++t)
                {
                  const std::array<int, 4> tet = polyhedron->subelement(t);
                  if (tet[0] < 0 || tet[1] < 0 || tet[2] < 0 || tet[3] < 0)
                    libmesh_not_implemented_msg
                      ("Cannot convert a C0Polyhedron whose triangulation\n"
                       "introduces special (non-vertex) points");
                  subelem[t] = Elem::build(TET4);
                  subelem[t]->set_node(0, elem->node_ptr(tet[0]));
                  subelem[t]->set_node(1, elem->node_ptr(tet[1]));
                  subelem[t]->set_node(2, elem->node_ptr(tet[2]));
                  subelem[t]->set_node(3, elem->node_ptr(tet[3]));
                }
              // There is a concern that two neighbor polyhedra might have
              // a triangulation of a side that does not match. But the
              // default triangulation is based on the side's triangulation
              // and the side element is supposed to be shared (that's why
              // shared pointers to polygons are used to build the polyhedra).
              // So the default one should work.
 
              break;
            }
 
            // No need to split elements that are already simplicial:
          case EDGE2:
          case EDGE3:
          case EDGE4:
          case TRI3:
          case TRI6:
          case TRI7:
          case TET4:
          case TET10:
          case TET14:
          case INFEDGE2:
            // No way to split infinite quad/prism elements, so
            // hopefully no need to
          case INFQUAD4:
          case INFQUAD6:
          case INFPRISM6:
          case INFPRISM12:
            continue;
            // If we're left with an unimplemented element we're
            // probably out of luck.  TODO: implement hex20, hex27,
            // pyramid13,...
          default:
            libmesh_not_implemented_msg
              ("Error, encountered unimplemented element "
               << Utility::enum_to_string<ElemType>(etype)
               << " in MeshTools::Modification::all_tri()...");
          } // end switch (etype)
 
        // Be sure the correct data is set for all subelems.
        const unsigned int nei = elem->n_extra_integers();
        for (unsigned int i=0; i != max_subelems; ++i)
          if (subelem[i]) {
            subelem[i]->processor_id() = elem->processor_id();
            subelem[i]->subdomain_id() = elem->subdomain_id();
 
            // Copy any extra element data.  Since the subelements
            // haven't been added to the mesh yet any allocation has
            // to be done manually.
            subelem[i]->add_extra_integers(nei);
            for (unsigned int ei=0; ei != nei; ++ei)
              subelem[ei]->set_extra_integer(ei, elem->get_extra_integer(ei));
 
 
            // Copy any mapping data.
            subelem[i]->set_mapping_type(elem->mapping_type());
            subelem[i]->set_mapping_data(elem->mapping_data());
          }
 
        // On a mesh with boundary data, we need to move that data to
        // the new elements.
 
        // On a mesh which is distributed, we need to move
        // remote_elem links to the new elements.
        bool mesh_is_serial = mesh.is_serial();
 
        if (mesh_has_boundary_data || !mesh_is_serial)
          {
            // Container to key boundary IDs handed back by the BoundaryInfo object.
            std::vector<boundary_id_type> bc_ids;
 
            for (auto sn : elem->side_index_range())
              {
                mesh.get_boundary_info().boundary_ids(elem, sn, bc_ids);
 
                if (bc_ids.empty() && elem->neighbor_ptr(sn) != remote_elem)
                  continue;
 
                // Make a sorted list of node ids for elem->side(sn)
                elem->build_side_ptr(elem_side, sn);
                std::vector<dof_id_type> elem_side_nodes(elem_side->n_nodes());
                for (unsigned int esn=0,
                     n_esn = cast_int<unsigned int>(elem_side_nodes.size());
                     esn != n_esn; ++esn)
                  elem_side_nodes[esn] = elem_side->node_id(esn);
                std::sort(elem_side_nodes.begin(), elem_side_nodes.end());
 
                for (unsigned int i=0; i != max_subelems; ++i)
                  if (subelem[i])
                    {
                      for (auto subside : subelem[i]->side_index_range())
                        {
                          subelem[i]->build_side_ptr(subside_elem, subside);
 
                          // Make a list of *vertex* node ids for this subside, see if they are all present
                          // in elem->side(sn).  Note 1: we can't just compare elem->key(sn) to
                          // subelem[i]->key(subside) in the Prism cases, since the new side is
                          // a different type.  Note 2: we only use vertex nodes since, in the future,
                          // a Hex20 or Prism15's QUAD8 face may be split into two Tri6 faces, and the
                          // original face will not contain the mid-edge node.
                          std::vector<dof_id_type> subside_nodes(subside_elem->n_vertices());
                          for (unsigned int ssn=0,
                               n_ssn = cast_int<unsigned int>(subside_nodes.size());
                               ssn != n_ssn; ++ssn)
                            subside_nodes[ssn] = subside_elem->node_id(ssn);
                          std::sort(subside_nodes.begin(), subside_nodes.end());
 
                          // std::includes returns true if every element of the second sorted range is
                          // contained in the first sorted range.
                          if (std::includes(elem_side_nodes.begin(), elem_side_nodes.end(),
                                            subside_nodes.begin(), subside_nodes.end()))
                            {
                              for (const auto & b_id : bc_ids)
                                if (b_id != BoundaryInfo::invalid_id)
                                  {
                                    new_bndry_ids.push_back(b_id);
                                    new_bndry_elements.push_back(subelem[i].get());
                                    new_bndry_sides.push_back(subside);
                                  }
 
                              // If the original element had a RemoteElem neighbor on side 'sn',
                              // then the subelem has one on side 'subside'.
                              if (elem->neighbor_ptr(sn) == remote_elem)
                                subelem[i]->set_neighbor(subside, const_cast<RemoteElem*>(remote_elem));
                            }
                        }
                    } // end for loop over subelem
              } // end for loop over sides
 
            // Remove the original element from the BoundaryInfo structure.
            mesh.get_boundary_info().remove(elem);
 
          } // end if (mesh_has_boundary_data)
 
        // Determine new IDs for the split elements which will be
        // the same on all processors, therefore keeping the Mesh
        // in sync.  Note: we offset the new IDs by max_orig_id to
        // avoid overwriting any of the original IDs.
        for (unsigned int i=0; i != max_subelems; ++i)
          if (subelem[i])
            {
              // Determine new IDs for the split elements which will be
              // the same on all processors, therefore keeping the Mesh
              // in sync.  Note: we offset the new IDs by the max of the
              // pre-existing ids to avoid conflicting with originals.
              subelem[i]->set_id( max_orig_id + max_subelems*elem->id() + i );
 
#ifdef LIBMESH_ENABLE_UNIQUE_ID
              subelem[i]->set_unique_id(max_unique_id + max_subelems*elem->unique_id() + i);
#endif
 
              // Prepare to add the newly-created simplices
              new_elements.push_back(std::move(subelem[i]));
            }
 
        // Delete the original element
        mesh.delete_elem(elem);
      } // End for loop over elements
  } // end scope
 
 
  // Now, iterate over the new elements vector, and add them each to
  // the Mesh.
  for (auto & elem : new_elements)
    mesh.add_elem(std::move(elem));
 
  if (mesh_has_boundary_data)
    {
      // If the old mesh had boundary data, the new mesh better have
      // some.  However, we can't assert that the size of
      // new_bndry_elements vector is > 0, since we may not have split
      // any elements actually on the boundary.  We also can't assert
      // that the original number of boundary sides is equal to the
      // sum of the boundary sides currently in the mesh and the
      // newly-added boundary sides, since in 3D, we may have split a
      // boundary QUAD into two boundary TRIs.  Therefore, we won't be
      // too picky about the actual number of BCs, and just assert that
      // there are some, somewhere.
#ifndef NDEBUG
      bool nbe_nonempty = new_bndry_elements.size();
      mesh.comm().max(nbe_nonempty);
      libmesh_assert(nbe_nonempty ||
                     mesh.get_boundary_info().n_boundary_conds()>0);
#endif
 
      // We should also be sure that the lengths of the new boundary data vectors
      // are all the same.
      libmesh_assert_equal_to (new_bndry_elements.size(), new_bndry_sides.size());
      libmesh_assert_equal_to (new_bndry_sides.size(), new_bndry_ids.size());
 
      // Add the new boundary info to the mesh
      for (auto s : index_range(new_bndry_elements))
        mesh.get_boundary_info().add_side(new_bndry_elements[s],
                                          new_bndry_sides[s],
                                          new_bndry_ids[s]);
    }
 
  // In a DistributedMesh any newly added ghost node ids may be
  // inconsistent, and unique_ids of newly added ghost nodes remain
  // unset.
  // make_nodes_parallel_consistent() will fix all this.
  if (!mesh.is_serial())
    {
      mesh.comm().max(added_new_ghost_point);
 
      if (added_new_ghost_point)
        MeshCommunication().make_nodes_parallel_consistent (mesh);
    }
 
  // Prepare the newly created mesh for use.
  mesh.prepare_for_use();
 
  // Let the new_elements and new_bndry_elements vectors go out of scope.
}

Referenced by OverlappingFunctorTest::checkCouplingFunctorTri(), OverlappingFunctorTest::checkCouplingFunctorTriUnifRef(), main(), AllTriTest::test_helper(), AllTriTest::test_helper_c0polyhedron(), AllRBBTest::test_sphere(), AllTriTest::testAllTriC0Polygon(), AllTriTest::testAllTriC0PolygonOctagon(), and libMesh::MeshTetInterface::volume_to_surface_mesh().

◆ 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 1921 of file mesh_modification.C.

{
  // This is just a shim around the member implementation, now
  mesh.get_boundary_info().renumber_id(old_id, new_id);
}

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

Referenced by MeshStitchTest::renameAndShift().

◆ 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 1931 of file mesh_modification.C.

{
  if (old_id == new_id)
    {
      // If the IDs are the same, this is a no-op.
      return;
    }
 
  Threads::parallel_for
    (mesh.element_stored_range(),
     [old_id, new_id](const ElemRange & range)
     {
       for (Elem * elem : range)
         if (elem->subdomain_id() == old_id)
           elem->subdomain_id() = new_id;
     });
 
  // We just invalidated mesh.get_subdomain_ids(), but it might not be
  // efficient to fix that here.
  mesh.unset_has_cached_elem_data();
}

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

◆ 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 148 of file mesh_modification.C.

{
  libmesh_assert (mesh.n_nodes());
  libmesh_assert (mesh.n_elem());
  libmesh_assert ((factor >= 0.) && (factor <= 1.));
 
  LOG_SCOPE("distort()", "MeshTools::Modification");
 
  // If we are not perturbing boundary nodes, make a
  // quickly-searchable list of node ids we can check against.
  std::unordered_set<dof_id_type> boundary_node_ids;
  if (!perturb_boundary)
    boundary_node_ids = MeshTools::find_boundary_nodes (mesh);
 
  // Now calculate the minimum distance to
  // neighboring nodes for each node.
  // hmin holds these distances.
  std::vector<float> hmin (mesh.max_node_id(),
                           std::numeric_limits<float>::max());
 
  for (const auto & elem : mesh.active_element_ptr_range())
    for (auto & n : elem->node_ref_range())
      hmin[n.id()] = std::min(hmin[n.id()],
                              static_cast<float>(elem->hmin()));
 
  // Now actually move the nodes
  {
    const unsigned int seed = 123456;
 
    // seed the random number generator.
    // We'll loop from 1 to n_nodes on every processor, even those
    // that don't have a particular node, so that the pseudorandom
    // numbers will be the same everywhere.
    std::srand(seed);
 
    // If the node is on the boundary or
    // the node is not used by any element (hmin[n]<1.e20)
    // then we should not move it.
    // [Note: Testing for (in)equality might be wrong
    // (different types, namely float and double)]
    for (auto n : make_range(mesh.max_node_id()))
      if ((perturb_boundary || !boundary_node_ids.count(n)) && hmin[n] < 1.e20)
        {
          // the direction, random but unit normalized
          Point dir (static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
                     (mesh.mesh_dimension() > 1) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.,
                     ((mesh.mesh_dimension() == 3) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.));
 
          dir(0) = (dir(0)-.5)*2.;
#if LIBMESH_DIM > 1
          if (mesh.mesh_dimension() > 1)
            dir(1) = (dir(1)-.5)*2.;
#endif
#if LIBMESH_DIM > 2
          if (mesh.mesh_dimension() == 3)
            dir(2) = (dir(2)-.5)*2.;
#endif
 
          dir = dir.unit();
 
          Node * node = mesh.query_node_ptr(n);
          if (!node)
            continue;
 
          (*node)(0) += dir(0)*factor*hmin[n];
#if LIBMESH_DIM > 1
          if (mesh.mesh_dimension() > 1)
            (*node)(1) += dir(1)*factor*hmin[n];
#endif
#if LIBMESH_DIM > 2
          if (mesh.mesh_dimension() == 3)
            (*node)(2) += dir(2)*factor*hmin[n];
#endif
        }
  }
 
  // We haven't changed any topology, but just changing geometry could
  // have invalidated a point locator.
  mesh.clear_point_locator();
}

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(), libMesh::Real, and libMesh::TypeVector< T >::unit().

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

◆ 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 1795 of file mesh_modification.C.

{
  libmesh_assert(mesh.is_prepared() || mesh.is_replicated());
 
  // Algorithm:
  // .) For each active element in the mesh: construct a
  //    copy which is the same in every way *except* it is
  //    a level 0 element.  Store the pointers to these in
  //    a separate vector. Save any boundary information as well.
  //    Delete the active element from the mesh.
  // .) Loop over all (remaining) elements in the mesh, delete them.
  // .) Add the level-0 copies back to the mesh
 
  // Temporary storage for new element pointers
  std::vector<std::unique_ptr<Elem>> new_elements;
 
  // BoundaryInfo Storage for element ids, sides, and BC ids
  std::vector<Elem *>              saved_boundary_elements;
  std::vector<boundary_id_type>   saved_bc_ids;
  std::vector<unsigned short int> saved_bc_sides;
 
  // Container to catch boundary ids passed back by BoundaryInfo
  std::vector<boundary_id_type> bc_ids;
 
  // Reserve a reasonable amt. of space for each
  new_elements.reserve(mesh.n_active_elem());
  saved_boundary_elements.reserve(mesh.get_boundary_info().n_boundary_conds());
  saved_bc_ids.reserve(mesh.get_boundary_info().n_boundary_conds());
  saved_bc_sides.reserve(mesh.get_boundary_info().n_boundary_conds());
 
  for (auto & elem : mesh.active_element_ptr_range())
    {
      // Make a new element of the same type
      auto copy = Elem::build(elem->type());
 
      // Set node pointers (they still point to nodes in the original mesh)
      for (auto n : elem->node_index_range())
        copy->set_node(n, elem->node_ptr(n));
 
      // Copy over ids
      copy->processor_id() = elem->processor_id();
      copy->subdomain_id() = elem->subdomain_id();
 
      // Retain the original element's ID(s) as well, otherwise
      // the Mesh may try to create them for you...
      copy->set_id( elem->id() );
#ifdef LIBMESH_ENABLE_UNIQUE_ID
      copy->set_unique_id(elem->unique_id());
#endif
 
      // This element could have boundary info or DistributedMesh
      // remote_elem links as well.  We need to save the (elem,
      // side, bc_id) triples and those links
      for (auto s : elem->side_index_range())
        {
          if (elem->neighbor_ptr(s) == remote_elem)
            copy->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
 
          mesh.get_boundary_info().boundary_ids(elem, s, bc_ids);
          for (const auto & bc_id : bc_ids)
            if (bc_id != BoundaryInfo::invalid_id)
              {
                saved_boundary_elements.push_back(copy.get());
                saved_bc_ids.push_back(bc_id);
                saved_bc_sides.push_back(s);
              }
        }
 
      // Copy any extra element data.  Since the copy hasn't been
      // added to the mesh yet any allocation has to be done manually.
      const unsigned int nei = elem->n_extra_integers();
      copy->add_extra_integers(nei);
      for (unsigned int i=0; i != nei; ++i)
        copy->set_extra_integer(i, elem->get_extra_integer(i));
 
      // Copy any mapping data.
      copy->set_mapping_type(elem->mapping_type());
      copy->set_mapping_data(elem->mapping_data());
 
      // We're done with this element
      mesh.delete_elem(elem);
 
      // But save the copy
      new_elements.push_back(std::move(copy));
    }
 
  // Make sure we saved the same number of boundary conditions
  // in each vector.
  libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
  libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
 
  // Loop again, delete any remaining elements
  for (auto & elem : mesh.element_ptr_range())
    mesh.delete_elem(elem);
 
  // Add the copied (now level-0) elements back to the mesh
  for (auto & new_elem : new_elements)
    {
      // Save the original ID, because the act of adding the Elem can
      // change new_elem's id!
      dof_id_type orig_id = new_elem->id();
 
      Elem * added_elem = mesh.add_elem(std::move(new_elem));
 
      // If the Elem, as it was re-added to the mesh, now has a
      // different ID (this is unlikely, so it's just an assert)
      // the boundary information will no longer be correct.
      libmesh_assert_equal_to (orig_id, added_elem->id());
 
      // Avoid compiler warnings in opt mode.
      libmesh_ignore(added_elem, orig_id);
    }
 
  // Finally, also add back the saved boundary information
  for (auto e : index_range(saved_boundary_elements))
    mesh.get_boundary_info().add_side(saved_boundary_elements[e],
                                      saved_bc_sides[e],
                                      saved_bc_ids[e]);
 
  // Trim unused and renumber nodes and elements
  mesh.prepare_for_use();
}

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().

◆ interpolate_surface()

void libMesh::MeshTools::Modification::interpolate_surface	(	MeshBase &	mesh,
		const Surface &	surface,
		std::set< std::size_t >	ids = `{}`,
		bool	use_boundary_nodes = `true`
	)

Move nodes in mesh to their closest points on the specified surface.

If use_boundary_nodes is true (the default), then ids is interpreted as a set of boundary side ids, and nodes are moved iff they are on element sides which have those boundary ids. We do not consider nodal boundary ids, which are for discrete points, which cannot be continuously moved without moving a neighborhood around them.

If use_boundary_nodes is false, then ids is interpreted as a set of subdomain ids, and nodes are moved iff they are on elements which are on those subdomains.

If ids is empty, then rather than "do nothing" we interpret that as "do everything". If use_boundary_nodes is true then we move all nodes on domain boundary sides (where elements have null neighbors), and if false then we move all nodes.

Definition at line 1695 of file mesh_modification.C.

{
  const bool is_serial = mesh.is_serial();
  const processor_id_type mesh_pid = mesh.processor_id();
 
  // We might have to move ghost nodes on a distributed mesh if their
  // owners don't see a requisite element or boundary they're on.
  std::unordered_set<dof_id_type> moved_ghost_nodes;
 
  auto move_node = [& moved_ghost_nodes, & surface, is_serial, mesh_pid]
                   (Node & node) {
    node = surface.closest_point(node);
 
    if (!is_serial && node.processor_id() != mesh_pid)
      moved_ghost_nodes.insert(node.id());
  };
 
  const bool no_ids = ids.empty();
  const BoundaryInfo & boundary_info = mesh.get_boundary_info();
 
  for (const auto & elem : mesh.active_element_ptr_range())
    {
      if (elem->mapping_type() != LAGRANGE_MAP)
        libmesh_not_implemented();
 
      if (use_boundary_nodes)
        {
          for (auto s : elem->side_index_range())
            {
              if (no_ids)
                {
                  // If we're not using boundary ids, we're
                  // interpolating all external and no internal
                  // boundaries
                  if (elem->neighbor_ptr(s))
                    continue;
                }
              else
                {
                  if (std::none_of(ids.begin(), ids.end(),
                                   [&boundary_info,elem,s](std::size_t bcid)
                                   {return boundary_info.has_boundary_id(elem, s, bcid);}))
                    continue;
                }
 
              for (auto n : elem->nodes_on_side(s))
                move_node(elem->node_ref(n));
            }
        }
      else
        {
          if (no_ids || ids.count(elem->subdomain_id()))
            for (Node & node : elem->node_ref_range())
              move_node(node);
        }
    }
 
  if (!is_serial)
    {
      std::map<processor_id_type, std::vector<dof_id_type>> moved_nodes_map;
      for (auto id : moved_ghost_nodes)
        {
          const Node & node = mesh.node_ref(id);
          moved_nodes_map[node.processor_id()].push_back(node.id());
        }
 
      auto action_functor =
        [& mesh, & surface]
        (processor_id_type /* pid */,
         const std::vector<dof_id_type> & my_moved_nodes)
        {
          for (auto id : my_moved_nodes)
            {
              Node & node = mesh.node_ref(id);
              node = surface.closest_point(node);
            }
        };
 
      // First get new node positions to their owners
      Parallel::push_parallel_vector_data
        (mesh.comm(), moved_nodes_map, action_functor);
 
      // Then get node positions to anyone else with them ghosted
      SyncNodalPositions sync_object(mesh);
      Parallel::sync_dofobject_data_by_id
        (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(),
         sync_object);
    }
 
  // We haven't changed any topology, but just changing geometry could
  // have invalidated a point locator.
  mesh.clear_point_locator();
}

References libMesh::MeshBase::clear_point_locator(), libMesh::Surface::closest_point(), libMesh::ParallelObject::comm(), libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::id(), libMesh::MeshBase::is_serial(), libMesh::LAGRANGE_MAP, libMesh::Elem::mapping_type(), mesh, libMesh::Elem::neighbor_ptr(), libMesh::MeshBase::node_ref(), libMesh::Elem::node_ref(), libMesh::Elem::node_ref_range(), libMesh::Elem::nodes_on_side(), libMesh::DofObject::processor_id(), libMesh::ParallelObject::processor_id(), TIMPI::push_parallel_vector_data(), libMesh::Elem::side_index_range(), libMesh::Elem::subdomain_id(), and libMesh::Parallel::sync_dofobject_data_by_id().

◆ 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 273 of file mesh_modification.C.

{
  LOG_SCOPE("orient_elements()", "MeshTools::Modification");
 
  // We don't yet support doing orient() on a parent element, which
  // would require us to consistently orient all its children and
  // give them different local child numbers.
  unsigned int n_levels = MeshTools::n_levels(mesh);
  if (n_levels > 1)
    libmesh_not_implemented_msg("orient_elements() does not support refined meshes");
 
  BoundaryInfo & boundary_info = mesh.get_boundary_info();
  for (auto elem : mesh.element_ptr_range())
    elem->orient(&boundary_info);
}

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

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

◆ 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 233 of file mesh_modification.C.

{
  LOG_SCOPE("permute_elements()", "MeshTools::Modification");
 
  // We don't yet support doing permute() on a parent element, which
  // would require us to consistently permute all its children and
  // give them different local child numbers.
  unsigned int n_levels = MeshTools::n_levels(mesh);
  if (n_levels > 1)
    libmesh_error();
 
  const unsigned int seed = 123456;
 
  // seed the random number generator.
  // We'll loop from 1 to max_elem_id on every processor, even those
  // that don't have a particular element, so that the pseudorandom
  // numbers will be the same everywhere.
  std::srand(seed);
 
 
  for (auto e_id : make_range(mesh.max_elem_id()))
    {
      int my_rand = std::rand();
 
      Elem * elem = mesh.query_elem_ptr(e_id);
 
      if (!elem)
        continue;
 
      const unsigned int max_permutation = elem->n_permutations();
      if (!max_permutation)
        continue;
 
      const unsigned int perm = my_rand % max_permutation;
 
      elem->permute(perm);
    }
}

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, build_nx, CaseName >::setUp().

◆ 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 291 of file mesh_modification.C.

{
  libmesh_assert (mesh.n_nodes());
  libmesh_assert (mesh.n_elem());
 
  LOG_SCOPE("redistribute()", "MeshTools::Modification");
 
  DenseVector<Real> output_vec(LIBMESH_DIM);
 
  // FIXME - we should thread this later.
  std::unique_ptr<FunctionBase<Real>> myfunc = mapfunc.clone();
 
  for (auto & node : mesh.node_ptr_range())
    {
      (*myfunc)(*node, output_vec);
 
      (*node)(0) = output_vec(0);
#if LIBMESH_DIM > 1
      (*node)(1) = output_vec(1);
#endif
#if LIBMESH_DIM > 2
      (*node)(2) = output_vec(2);
#endif
    }
 
  // If we just moved a mesh in or out out of the X axis or XY plane
  // then we might have changed its spatial_dimension()
  mesh.unset_has_cached_elem_data();
 
  // We haven't changed any topology, but just changing geometry could
  // have invalidated a point locator.
  mesh.clear_point_locator();
}

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

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

◆ rotate()

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

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 371 of file mesh_modification.C.

{
  // We won't change any topology, but just changing geometry could
  // invalidate a point locator.
  mesh.clear_point_locator();
 
#if LIBMESH_DIM == 3
  const auto R = RealTensorValue::intrinsic_rotation_matrix(phi, theta, psi);
 
  for (auto & node : mesh.node_ptr_range())
    {
      Point & pt = *node;
      pt = R * pt;
    }
 
  // If we just moved a mesh in or out out of the X axis or XY plane
  // then we might have changed its spatial_dimension()
  mesh.unset_has_cached_elem_data();
 
  return R;
 
#else
  libmesh_ignore(mesh, phi, theta, psi);
  libmesh_error_msg("MeshTools::Modification::rotate() requires libMesh to be compiled with LIBMESH_DIM==3");
  // We'll never get here
  return RealTensorValue();
#endif
}

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

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

◆ 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 404 of file mesh_modification.C.

{
  const Real x_scale = xs;
  Real y_scale       = ys;
  Real z_scale       = zs;
 
  if (ys == 0.)
    {
      libmesh_assert_equal_to (zs, 0.);
 
      y_scale = z_scale = x_scale;
    }
 
  // Scale the x coordinate in all dimensions
  for (auto & node : mesh.node_ptr_range())
    (*node)(0) *= x_scale;
 
  // Only scale the y coordinate in 2 and 3D
  if (LIBMESH_DIM < 2)
    return;
 
  for (auto & node : mesh.node_ptr_range())
    (*node)(1) *= y_scale;
 
  // Only scale the z coordinate in 3D
  if (LIBMESH_DIM < 3)
    return;
 
  for (auto & node : mesh.node_ptr_range())
    (*node)(2) *= z_scale;
 
  // If we just collapsed a manifold onto the X axis or XY plane
  // then we might have changed its spatial_dimension()
  mesh.unset_has_cached_elem_data();
 
  // We haven't changed any topology, but just changing geometry could
  // have invalidated a point locator.
  mesh.clear_point_locator();
}

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

Referenced by main().

◆ 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 1536 of file mesh_modification.C.

{
  libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
 
  /*
   * Create a quickly-searchable list of boundary nodes.
   */
  std::unordered_set<dof_id_type> boundary_node_ids =
    MeshTools::find_boundary_nodes (mesh);
 
  // For avoiding extraneous element side allocation
  ElemSideBuilder side_builder;
 
  for (unsigned int iter=0; iter<n_iterations; iter++)
    {
      /*
       * loop over the mesh refinement level
       */
      unsigned int n_levels = MeshTools::n_levels(mesh);
      for (unsigned int refinement_level=0; refinement_level != n_levels;
           refinement_level++)
        {
          // initialize the storage (have to do it on every level to get empty vectors
          std::vector<Point> new_positions;
          std::vector<Real>   weight;
          new_positions.resize(mesh.n_nodes());
          weight.resize(mesh.n_nodes());
 
          {
            // Loop over the elements to calculate new node positions
            for (const auto & elem : as_range(mesh.level_elements_begin(refinement_level),
                                              mesh.level_elements_end(refinement_level)))
              {
                /*
                 * We relax all nodes on level 0 first
                 * If the element is refined (level > 0), we interpolate the
                 * parents nodes with help of the embedding matrix
                 */
                if (refinement_level == 0)
                  {
                    for (auto s : elem->side_index_range())
                      {
                        /*
                         * Only operate on sides which are on the
                         * boundary or for which the current element's
                         * id is greater than its neighbor's.
                         * Sides get only built once.
                         */
                        if ((elem->neighbor_ptr(s) != nullptr) &&
                            (elem->id() > elem->neighbor_ptr(s)->id()))
                          {
                            const Elem & side = side_builder(*elem, s);
                            const Node & node0 = side.node_ref(0);
                            const Node & node1 = side.node_ref(1);
 
                            Real node_weight = 1.;
                            // calculate the weight of the nodes
                            if (power > 0)
                              {
                                Point diff = node0-node1;
                                node_weight = std::pow(diff.norm(), power);
                              }
 
                            const dof_id_type id0 = node0.id(), id1 = node1.id();
                            new_positions[id0].add_scaled( node1, node_weight );
                            new_positions[id1].add_scaled( node0, node_weight );
                            weight[id0] += node_weight;
                            weight[id1] += node_weight;
                          }
                      } // element neighbor loop
                  }
#ifdef LIBMESH_ENABLE_AMR
                else   // refinement_level > 0
                  {
                    /*
                     * Find the positions of the hanging nodes of refined elements.
                     * We do this by calculating their position based on the parent
                     * (one level less refined) element, and the embedding matrix
                     */
 
                    const Elem * parent = elem->parent();
 
                    /*
                     * find out which child I am
                     */
                    unsigned int c = parent->which_child_am_i(elem);
                    /*
                     *loop over the childs (that is, the current elements) nodes
                     */
                    for (auto nc : elem->node_index_range())
                      {
                        /*
                         * the new position of the node
                         */
                        Point point;
                        for (auto n : parent->node_index_range())
                          {
                            /*
                             * The value from the embedding matrix
                             */
                            const Real em_val = parent->embedding_matrix(c,nc,n);
 
                            if (em_val != 0.)
                              point.add_scaled (parent->point(n), em_val);
                          }
 
                        const dof_id_type id = elem->node_ptr(nc)->id();
                        new_positions[id] = point;
                        weight[id] = 1.;
                      }
                  } // if element refinement_level
#endif // #ifdef LIBMESH_ENABLE_AMR
 
              } // element loop
 
            /*
             * finally reposition the vertex nodes
             */
            for (auto nid : make_range(mesh.n_nodes()))
              if (!boundary_node_ids.count(nid) && weight[nid] > 0.)
                mesh.node_ref(nid) = new_positions[nid]/weight[nid];
          }
 
          // Now handle the additional second_order nodes by calculating
          // their position based on the vertex positions
          // we do a second loop over the level elements
          for (auto & elem : as_range(mesh.level_elements_begin(refinement_level),
                                      mesh.level_elements_end(refinement_level)))
            {
              const unsigned int son_begin = elem->n_vertices();
              const unsigned int son_end   = elem->n_nodes();
              for (unsigned int n=son_begin; n<son_end; n++)
                {
                  const unsigned int n_adjacent_vertices =
                    elem->n_second_order_adjacent_vertices(n);
 
                  Point point;
                  for (unsigned int v=0; v<n_adjacent_vertices; v++)
                    point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
 
                  const dof_id_type id = elem->node_ptr(n)->id();
                  mesh.node_ref(id) = point/n_adjacent_vertices;
                }
            }
        } // refinement_level loop
    } // end iteration
 
  // We haven't changed any topology, but just changing geometry could
  // have invalidated a point locator.
  mesh.clear_point_locator();
}

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::Elem::n_nodes(), libMesh::MeshBase::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::MeshBase::node_ref(), libMesh::Elem::node_ref(), libMesh::TypeVector< T >::norm(), libMesh::Elem::parent(), libMesh::Elem::point(), libMesh::Real, libMesh::Elem::second_order_adjacent_vertex(), libMesh::Elem::side_index_range(), libMesh::MeshTools::weight(), and libMesh::Elem::which_child_am_i().

◆ 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 328 of file mesh_modification.C.

{
  const Point p(xt, yt, zt);
 
  for (auto & node : mesh.node_ptr_range())
    *node += p;
 
  // If we just moved a mesh in or out out of the X axis or XY plane
  // then we might have changed its spatial_dimension()
  mesh.unset_has_cached_elem_data();
 
  // We haven't changed any topology, but just changing geometry could
  // have invalidated a point locator.
  mesh.clear_point_locator();
}

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

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

Functions