libMesh::MeshTools::Modification Namespace Reference

Functions
void	distort (MeshBase &mesh, const Real factor, const bool perturb_boundary=false)
void	translate (MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.)
void	rotate (MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
void	scale (MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
void	all_tri (MeshBase &mesh)
void	smooth (MeshBase &, unsigned int, Real)
void	flatten (MeshBase &mesh)
void	change_boundary_id (MeshBase &mesh, const boundary_id_type old_id, const boundary_id_type new_id)
void	change_subdomain_id (MeshBase &mesh, const subdomain_id_type old_id, const subdomain_id_type new_id)

---

Detailed Description

Tools for Mesh modification.

Author:

Benjamin S. Kirk

Date:

2004

---

Function Documentation

void libMesh::MeshTools::Modification::all_tri

(

MeshBase &

mesh

)

Converts the 2D quadrilateral elements of a Mesh into triangular elements. Note: Only works for 2D elements! 3D elements are ignored. Note: Probably won't do the right thing for meshes which have been refined previously.

Definition at line 686 of file mesh_modification.C.

References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::MeshBase::boundary_info, libMesh::MeshBase::delete_elem(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), end, libMesh::Utility::enum_to_string< ElemType >(), libMesh::err, libMesh::Elem::get_node(), libMesh::DofObject::id(), libMesh::BoundaryInfo::invalid_id, libMesh::MeshBase::n_elem(), libMesh::Elem::n_sides(), libMesh::Elem::neighbor(), libMesh::Elem::node(), libMesh::MeshBase::node(), libMesh::Elem::parent(), libMesh::Elem::point(), libMesh::MeshBase::prepare_for_use(), libMesh::DofObject::processor_id(), libMeshEnums::QUAD4, libMeshEnums::QUAD8, libMeshEnums::QUAD9, libMesh::remote_elem, libMesh::DofObject::set_id(), libMesh::Elem::set_neighbor(), libMesh::Elem::set_node(), libMesh::Elem::subdomain_id(), libMeshEnums::TRI3, libMeshEnums::TRI6, and libMesh::Elem::type().

{
  // The number of elements in the original mesh before any additions
  // or deletions.
  const dof_id_type n_orig_elem = mesh.n_elem();

  // 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 iterators.
  std::vector<Elem*> new_elements;
  new_elements.reserve (2*n_orig_elem);

  // If the original mesh has boundary data, we carry that over
  // to the new mesh with triangular elements.
  const bool mesh_has_boundary_data = (mesh.boundary_info->n_boundary_ids() > 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;

  // Iterate over the elements, splitting QUADS into
  // pairs of conforming triangles.
  // FIXME: This algorithm does not work on refined grids!
  {
    MeshBase::element_iterator       el  = mesh.elements_begin();
    const MeshBase::element_iterator end = mesh.elements_end();

    for (; el!=end; ++el)
      {
        Elem* elem = *el;

        const ElemType etype = elem->type();

        // all_tri currently only works on coarse meshes
        libmesh_assert (!elem->parent());

        // We split the quads using the shorter of the two diagonals
        // to maintain the best angle properties.
        bool edge_swap = false;

        // True if we actually split the current element.
        bool split_elem = false;

        // The two new triangular elements we will split the quad into.
        Elem* tri0 = NULL;
        Elem* tri1 = NULL;


        switch (etype)
          {
          case QUAD4:
            {
              split_elem = true;

              tri0 = new Tri3;
              tri1 = new Tri3;

              // Check for possible edge swap
              if ((elem->point(0) - elem->point(2)).size() <
                  (elem->point(1) - elem->point(3)).size())
                {
                  tri0->set_node(0) = elem->get_node(0);
                  tri0->set_node(1) = elem->get_node(1);
                  tri0->set_node(2) = elem->get_node(2);

                  tri1->set_node(0) = elem->get_node(0);
                  tri1->set_node(1) = elem->get_node(2);
                  tri1->set_node(2) = elem->get_node(3);
                }

              else
                {
                  edge_swap=true;

                  tri0->set_node(0) = elem->get_node(0);
                  tri0->set_node(1) = elem->get_node(1);
                  tri0->set_node(2) = elem->get_node(3);

                  tri1->set_node(0) = elem->get_node(1);
                  tri1->set_node(1) = elem->get_node(2);
                  tri1->set_node(2) = elem->get_node(3);
                }


              break;
            }

          case QUAD8:
            {
              split_elem =  true;

              tri0 = new Tri6;
              tri1 = new Tri6;

              Node* new_node = mesh.add_point( (mesh.node(elem->node(0)) +
                                                mesh.node(elem->node(1)) +
                                                mesh.node(elem->node(2)) +
                                                mesh.node(elem->node(3)) / 4)
                                               );

              // Check for possible edge swap
              if ((elem->point(0) - elem->point(2)).size() <
                  (elem->point(1) - elem->point(3)).size())
                {
                  tri0->set_node(0) = elem->get_node(0);
                  tri0->set_node(1) = elem->get_node(1);
                  tri0->set_node(2) = elem->get_node(2);
                  tri0->set_node(3) = elem->get_node(4);
                  tri0->set_node(4) = elem->get_node(5);
                  tri0->set_node(5) = new_node;

                  tri1->set_node(0) = elem->get_node(0);
                  tri1->set_node(1) = elem->get_node(2);
                  tri1->set_node(2) = elem->get_node(3);
                  tri1->set_node(3) = new_node;
                  tri1->set_node(4) = elem->get_node(6);
                  tri1->set_node(5) = elem->get_node(7);

                }

              else
                {
                  edge_swap=true;

                  tri0->set_node(0) = elem->get_node(3);
                  tri0->set_node(1) = elem->get_node(0);
                  tri0->set_node(2) = elem->get_node(1);
                  tri0->set_node(3) = elem->get_node(7);
                  tri0->set_node(4) = elem->get_node(4);
                  tri0->set_node(5) = new_node;

                  tri1->set_node(0) = elem->get_node(1);
                  tri1->set_node(1) = elem->get_node(2);
                  tri1->set_node(2) = elem->get_node(3);
                  tri1->set_node(3) = elem->get_node(5);
                  tri1->set_node(4) = elem->get_node(6);
                  tri1->set_node(5) = new_node;
                }

              break;
            }

          case QUAD9:
            {
              split_elem =  true;

              tri0 = new Tri6;
              tri1 = new Tri6;

              // Check for possible edge swap
              if ((elem->point(0) - elem->point(2)).size() <
                  (elem->point(1) - elem->point(3)).size())
                {
                  tri0->set_node(0) = elem->get_node(0);
                  tri0->set_node(1) = elem->get_node(1);
                  tri0->set_node(2) = elem->get_node(2);
                  tri0->set_node(3) = elem->get_node(4);
                  tri0->set_node(4) = elem->get_node(5);
                  tri0->set_node(5) = elem->get_node(8);

                  tri1->set_node(0) = elem->get_node(0);
                  tri1->set_node(1) = elem->get_node(2);
                  tri1->set_node(2) = elem->get_node(3);
                  tri1->set_node(3) = elem->get_node(8);
                  tri1->set_node(4) = elem->get_node(6);
                  tri1->set_node(5) = elem->get_node(7);
                }

              else
                {
                  edge_swap=true;

                  tri0->set_node(0) = elem->get_node(0);
                  tri0->set_node(1) = elem->get_node(1);
                  tri0->set_node(2) = elem->get_node(3);
                  tri0->set_node(3) = elem->get_node(4);
                  tri0->set_node(4) = elem->get_node(8);
                  tri0->set_node(5) = elem->get_node(7);

                  tri1->set_node(0) = elem->get_node(1);
                  tri1->set_node(1) = elem->get_node(2);
                  tri1->set_node(2) = elem->get_node(3);
                  tri1->set_node(3) = elem->get_node(5);
                  tri1->set_node(4) = elem->get_node(6);
                  tri1->set_node(5) = elem->get_node(8);
                }

              break;
            }
          // No need to split elements that are already triangles
          case TRI3:
          case TRI6:
            continue;
          // Try to ignore non-2D elements for now
          default:
            {
              libMesh::err << "Warning, encountered non-2D element "
                            << Utility::enum_to_string<ElemType>(etype)
                            << " in MeshTools::Modification::all_tri(), hope that's OK..."
                            << std::endl;
            }
          } // end switch (etype)



        if (split_elem)
          {
            // Be sure the correct ID's are also set for tri0 and
            // tri1.
            tri0->processor_id() = elem->processor_id();
            tri0->subdomain_id() = elem->subdomain_id();
            tri1->processor_id() = elem->processor_id();
            tri1->subdomain_id() = elem->subdomain_id();

            if (mesh_has_boundary_data)
              {
                for (unsigned int sn=0; sn<elem->n_sides(); ++sn)
                  {
                    const std::vector<boundary_id_type>& bc_ids = mesh.boundary_info->boundary_ids(*el, sn);
                    for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
                      {
                        const boundary_id_type b_id = *id_it;

                        if (b_id != BoundaryInfo::invalid_id)
                          {
                            // Add the boundary ID to the list of new boundary ids
                            new_bndry_ids.push_back(b_id);

                            // Convert the boundary side information of the old element to
                            // boundary side information for the new element.
                            if (!edge_swap)
                              {
                                switch (sn)
                                  {
                                  case 0:
                                    {
                                      // New boundary side is Tri 0, side 0
                                      new_bndry_elements.push_back(tri0);
                                      new_bndry_sides.push_back(0);
                                      break;
                                    }
                                  case 1:
                                    {
                                      // New boundary side is Tri 0, side 1
                                      new_bndry_elements.push_back(tri0);
                                      new_bndry_sides.push_back(1);
                                      break;
                                    }
                                  case 2:
                                    {
                                      // New boundary side is Tri 1, side 1
                                      new_bndry_elements.push_back(tri1);
                                      new_bndry_sides.push_back(1);
                                      break;
                                    }
                                  case 3:
                                    {
                                      // New boundary side is Tri 1, side 2
                                      new_bndry_elements.push_back(tri1);
                                      new_bndry_sides.push_back(2);
                                      break;
                                    }

                                  default:
                                    {
                                      libMesh::err << "Quad4/8/9 cannot have more than 4 sides." << std::endl;
                                      libmesh_error();
                                    }
                                  }
                              }

                            else // edge_swap==true
                              {
                                switch (sn)
                                  {
                                  case 0:
                                    {
                                      // New boundary side is Tri 0, side 0
                                      new_bndry_elements.push_back(tri0);
                                      new_bndry_sides.push_back(0);
                                      break;
                                    }
                                  case 1:
                                    {
                                      // New boundary side is Tri 1, side 0
                                      new_bndry_elements.push_back(tri1);
                                      new_bndry_sides.push_back(0);
                                      break;
                                    }
                                  case 2:
                                    {
                                      // New boundary side is Tri 1, side 1
                                      new_bndry_elements.push_back(tri1);
                                      new_bndry_sides.push_back(1);
                                      break;
                                    }
                                  case 3:
                                    {
                                      // New boundary side is Tri 0, side 2
                                      new_bndry_elements.push_back(tri0);
                                      new_bndry_sides.push_back(2);
                                      break;
                                    }

                                  default:
                                    {
                                      libMesh::err << "Quad4/8/9 cannot have more than 4 sides." << std::endl;
                                      libmesh_error();
                                    }
                                  }
                              } // end edge_swap==true
                          } // end if (b_id != BoundaryInfo::invalid_id)
                      } // end for loop over boundary IDs
                  } // end for loop over sides

                // Remove the original element from the BoundaryInfo structure.
                mesh.boundary_info->remove(elem);

              } // end if (mesh_has_boundary_data)


            // On a distributed mesh, we need to preserve remote_elem
            // links, since prepare_for_use can't reconstruct them for
            // us.
            for (unsigned int sn=0; sn<elem->n_sides(); ++sn)
              {
                if (elem->neighbor(sn) == remote_elem)
                  {
                    // Create a remote_elem link on one of the new
                    // elements corresponding to the link from the old
                    // element.
                    if (!edge_swap)
                      {
                        switch (sn)
                              {
                              case 0:
                                {
                                  // New remote side is Tri 0, side 0
                                  tri0->set_neighbor(0, const_cast<RemoteElem*>(remote_elem));
                                  break;
                                }
                              case 1:
                                {
                                  // New remote side is Tri 0, side 1
                                  tri0->set_neighbor(1, const_cast<RemoteElem*>(remote_elem));
                                  break;
                                }
                              case 2:
                                {
                                  // New remote side is Tri 1, side 1
                                  tri1->set_neighbor(1, const_cast<RemoteElem*>(remote_elem));
                                  break;
                                }
                              case 3:
                                {
                                  // New remote side is Tri 1, side 2
                                  tri1->set_neighbor(2, const_cast<RemoteElem*>(remote_elem));
                                  break;
                                }

                              default:
                                {
                                  libMesh::err << "Quad4/8/9 cannot have more than 4 sides." << std::endl;
                                  libmesh_error();
                                }
                              }
                          }

                        else // edge_swap==true
                          {
                            switch (sn)
                              {
                              case 0:
                                {
                                  // New remote side is Tri 0, side 0
                                  tri0->set_neighbor(0, const_cast<RemoteElem*>(remote_elem));
                                  break;
                                }
                              case 1:
                                {
                                  // New remote side is Tri 1, side 0
                                  tri1->set_neighbor(0, const_cast<RemoteElem*>(remote_elem));
                                  break;
                                }
                              case 2:
                                {
                                  // New remote side is Tri 1, side 1
                                  tri1->set_neighbor(1, const_cast<RemoteElem*>(remote_elem));
                                  break;
                                }
                              case 3:
                                {
                                  // New remote side is Tri 0, side 2
                                  tri0->set_neighbor(2, const_cast<RemoteElem*>(remote_elem));
                                  break;
                                }

                              default:
                                {
                                  libMesh::err << "Quad4/8/9 cannot have more than 4 sides." << std::endl;
                                  libmesh_error();
                                }
                              }
                          } // end edge_swap==true
                      } // end if (elem->neighbor(sn) == remote_elem)
              } // end for loop over sides

            // 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 n_orig_elem to
            // avoid overwriting any of the original IDs, this assumes
            // they were contiguously-numbered to begin with...
            tri0->set_id( n_orig_elem + 2*elem->id() + 0 );
            tri1->set_id( n_orig_elem + 2*elem->id() + 1 );

            // Add the newly-created triangles to the temporary vector of new elements.
            new_elements.push_back(tri0);
            new_elements.push_back(tri1);

            // Delete the original element
            mesh.delete_elem(elem);
          } // end if (split_elem)
      } // End for loop over elements
  } // end scope


  // Now, iterate over the new elements vector, and add them each to
  // the Mesh.
  {
    std::vector<Elem*>::iterator el        = new_elements.begin();
    const std::vector<Elem*>::iterator end = new_elements.end();
    for (; el != end; ++el)
      mesh.add_elem(*el);
  }

  if (mesh_has_boundary_data)
    {
      // By this time, we should have removed all of the original boundary sides
      // - except on a hybrid mesh, where we can't "start from a blank slate"! - RHS
      // libmesh_assert_equal_to (mesh.boundary_info->n_boundary_conds(), 0);

      // Clear the boundary info, to be sure and start from a blank slate.
      // mesh.boundary_info->clear();

      // If the old mesh had boundary data, the new mesh better have some.
      libmesh_assert_greater (new_bndry_elements.size(), 0);

      // 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 (unsigned int s=0; s<new_bndry_elements.size(); ++s)
        mesh.boundary_info->add_side(new_bndry_elements[s],
                                     new_bndry_sides[s],
                                     new_bndry_ids[s]);
    }


  // Prepare the newly created mesh for use.
  mesh.prepare_for_use(/*skip_renumber =*/ false);

  // Let the new_elements and new_bndry_elements vectors go out of scope.
}

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

References libMesh::MeshBase::boundary_info, libMesh::MeshBase::level_elements_begin(), libMesh::MeshBase::level_elements_end(), and libMesh::Elem::n_sides().

{
  // Only level-0 elements store BCs.  Loop over them.
  MeshBase::element_iterator           el = mesh.level_elements_begin(0);
  const MeshBase::element_iterator end_el = mesh.level_elements_end(0);
  for (; el != end_el; ++el)
    {
      Elem *elem = *el;
      unsigned int n_sides = elem->n_sides();
      for (unsigned int s=0; s != n_sides; ++s)
        {
          const std::vector<boundary_id_type>& old_ids = mesh.boundary_info->boundary_ids(elem, s);
          if (std::find(old_ids.begin(), old_ids.end(), old_id) != old_ids.end())
            {
              std::vector<boundary_id_type> new_ids(old_ids);
              std::replace(new_ids.begin(), new_ids.end(), old_id, new_id);
              mesh.boundary_info->remove_side(elem, s);
              mesh.boundary_info->add_side(elem, s, new_ids);
            }
        }
    }
}

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

References libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), and libMesh::Elem::subdomain_id().

{
  MeshBase::element_iterator           el = mesh.elements_begin();
  const MeshBase::element_iterator end_el = mesh.elements_end();

  for (; el != end_el; ++el)
    {
      Elem *elem = *el;

      if (elem->subdomain_id() == old_id)
        elem->subdomain_id() = new_id;
    }
}

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

References libMesh::MeshBase::active_elements_begin(), libMesh::MeshBase::active_elements_end(), end, libMesh::MeshTools::find_boundary_nodes(), std::max(), libMesh::MeshBase::mesh_dimension(), std::min(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::MeshBase::node_ptr(), and libMesh::TypeVector< T >::unit().

{
  libmesh_assert (mesh.n_nodes());
  libmesh_assert (mesh.n_elem());
  libmesh_assert ((factor >= 0.) && (factor <= 1.));

  START_LOG("distort()", "MeshTools::Modification");



  // First find nodes on the boundary and flag them
  // so that we don't move them
  // on_boundary holds false (not on boundary) and true (on boundary)
  std::vector<bool> on_boundary (mesh.n_nodes(), false);

  if (!perturb_boundary) MeshTools::find_boundary_nodes (mesh, on_boundary);

  // Now calculate the minimum distance to
  // neighboring nodes for each node.
  // hmin holds these distances.
  std::vector<float> hmin (mesh.n_nodes(),
                           std::numeric_limits<float>::max());

  MeshBase::element_iterator       el  = mesh.active_elements_begin();
  const MeshBase::element_iterator end = mesh.active_elements_end();

  for (; el!=end; ++el)
    for (unsigned int n=0; n<(*el)->n_nodes(); n++)
      hmin[(*el)->node(n)] = std::min(hmin[(*el)->node(n)],
                                      static_cast<float>((*el)->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 (unsigned int n=0; n<mesh.n_nodes(); n++)
      if (!on_boundary[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 (mesh.mesh_dimension() > 1)
            dir(1) = (dir(1)-.5)*2.;
          if (mesh.mesh_dimension() == 3)
            dir(2) = (dir(2)-.5)*2.;

          dir = dir.unit();

          Node *node = mesh.node_ptr(n);
          if (!node)
            continue;

          (*node)(0) += dir(0)*factor*hmin[n];
          if (mesh.mesh_dimension() > 1)
            (*node)(1) += dir(1)*factor*hmin[n];
          if (mesh.mesh_dimension() == 3)
            (*node)(2) += dir(2)*factor*hmin[n];
        }
  }


  // All done
  STOP_LOG("distort()", "MeshTools::Modification");
}

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

References libMesh::MeshBase::active_elements_begin(), libMesh::MeshBase::active_elements_end(), libMesh::MeshBase::add_elem(), bc_id, libMesh::MeshBase::boundary_info, libMesh::Elem::build(), libMesh::MeshBase::delete_elem(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), end, libMesh::Elem::get_node(), libMesh::DofObject::id(), libMesh::BoundaryInfo::invalid_id, libMesh::MeshBase::n_active_elem(), libMesh::Elem::n_nodes(), libMesh::Elem::n_sides(), libMesh::Elem::neighbor(), libMesh::MeshBase::prepare_for_use(), libMesh::DofObject::processor_id(), libMesh::remote_elem, libMesh::DofObject::set_id(), libMesh::Elem::set_neighbor(), libMesh::Elem::set_node(), libMesh::Elem::subdomain_id(), and libMesh::Elem::type().

{
  // 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<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;

  // Reserve a reasonable amt. of space for each
  new_elements.reserve(mesh.n_active_elem());
  saved_boundary_elements.reserve(mesh.boundary_info->n_boundary_conds());
  saved_bc_ids.reserve(mesh.boundary_info->n_boundary_conds());
  saved_bc_sides.reserve(mesh.boundary_info->n_boundary_conds());
  {
    MeshBase::element_iterator       it  = mesh.active_elements_begin();
    const MeshBase::element_iterator end = mesh.active_elements_end();

    for (; it != end; ++it)
      {
        Elem* elem = *it;

        // Make a new element of the same type
        Elem* copy = Elem::build(elem->type()).release();

        // Set node pointers (they still point to nodes in the original mesh)
        for(unsigned int n=0; n<elem->n_nodes(); n++)
          copy->set_node(n) = elem->get_node(n);

        // Copy over ids
        copy->processor_id() = elem->processor_id();
        copy->subdomain_id() = elem->subdomain_id();

        // Retain the original element's ID as well, otherwise ParallelMesh will
        // try to create one for you...
        copy->set_id( elem->id() );

        // This element could have boundary info or ParallelMesh
        // remote_elem links as well.  We need to save the (elem,
        // side, bc_id) triples and those links
        for (unsigned int s=0; s<elem->n_sides(); s++)
          {
            if (elem->neighbor(s) == remote_elem)
              copy->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));

            const std::vector<boundary_id_type>& bc_ids = mesh.boundary_info->boundary_ids(elem,s);
            for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
              {
                const boundary_id_type bc_id = *id_it;

                if (bc_id != BoundaryInfo::invalid_id)
                  {
                    saved_boundary_elements.push_back(copy);
                    saved_bc_ids.push_back(bc_id);
                    saved_bc_sides.push_back(s);
                  }
              }
          }


        // We're done with this element
        mesh.delete_elem(elem);

        // But save the copy
        new_elements.push_back(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
  {
    MeshBase::element_iterator       it  = mesh.elements_begin();
    const MeshBase::element_iterator end = mesh.elements_end();

    for (; it != end; ++it)
      mesh.delete_elem( *it );
  }


  // Add the copied (now level-0) elements back to the mesh
  {
    for (std::vector<Elem*>::iterator it = new_elements.begin();
         it != new_elements.end();
         ++it)
      {
        dof_id_type orig_id = (*it)->id();

        Elem* added_elem = mesh.add_elem(*it);

        dof_id_type added_id = added_elem->id();

        // 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_id);
      }
  }

  // Finally, also add back the saved boundary information
  for (unsigned int e=0; e<saved_boundary_elements.size(); ++e)
    mesh.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(/*skip_renumber =*/ false);
}

void 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)

Definition at line 174 of file mesh_modification.C.

References libMesh::MeshBase::mesh_dimension(), libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::pi, and libMesh::Real.

{
  libmesh_assert_not_equal_to (mesh.mesh_dimension(), 1);

  const Real  p = -phi/180.*libMesh::pi;
  const Real  t = -theta/180.*libMesh::pi;
  const Real  s = -psi/180.*libMesh::pi;
  const Real sp = std::sin(p), cp = std::cos(p);
  const Real st = std::sin(t), ct = std::cos(t);
  const Real ss = std::sin(s), cs = std::cos(s);

  // We follow the convention described at http://mathworld.wolfram.com/EulerAngles.html
  // (equations 6-14 give the entries of the composite transformation matrix).
  // The rotations are performed sequentially about the z, x, and z axes, in that order.
  // A positive angle yields a counter-clockwise rotation about the axis in question.
  const MeshBase::node_iterator nd_end = mesh.nodes_end();

  for (MeshBase::node_iterator nd = mesh.nodes_begin();
       nd != nd_end; ++nd)
    {
      const Point pt = **nd;
      const Real  x  = pt(0);
      const Real  y  = pt(1);
      const Real  z  = pt(2);
      **nd = Point(( cp*cs-sp*ct*ss)*x + ( sp*cs+cp*ct*ss)*y + (st*ss)*z,
                   (-cp*ss-sp*ct*cs)*x + (-sp*ss+cp*ct*cs)*y + (st*cs)*z,
                   ( sp*st)*x          + (-cp*st)*y          + (ct)*z   );
    }
}

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

References libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::Real, and libMesh::MeshBase::spatial_dimension().

{
  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
  const MeshBase::node_iterator nd_end = mesh.nodes_end();

  for (MeshBase::node_iterator nd = mesh.nodes_begin();
       nd != nd_end; ++nd)
    (**nd)(0) *= x_scale;


  // Only scale the y coordinate in 2 and 3D
  if (mesh.spatial_dimension() < 2)
    return;

  for (MeshBase::node_iterator nd = mesh.nodes_begin();
       nd != nd_end; ++nd)
    (**nd)(1) *= y_scale;

  // Only scale the z coordinate in 3D
  if (mesh.spatial_dimension() < 3)
    return;

  for (MeshBase::node_iterator nd = mesh.nodes_begin();
       nd != nd_end; ++nd)
    (**nd)(2) *= z_scale;
}

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 postition 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 Lüthi (luthi@gi.alaska.edu)

Date:

2005

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

Definition at line 1155 of file mesh_modification.C.

References libMesh::TypeVector< T >::add(), libMesh::TypeVector< T >::add_scaled(), libMesh::Elem::build_side(), libMesh::Elem::child(), libMesh::Elem::embedding_matrix(), end, libMesh::MeshTools::find_boundary_nodes(), libMesh::Elem::get_node(), libMesh::DofObject::id(), libMesh::MeshBase::level_elements_begin(), libMesh::MeshBase::level_elements_end(), libMesh::MeshBase::mesh_dimension(), libMesh::Elem::n_children(), libMesh::MeshTools::n_levels(), libMesh::Elem::n_neighbors(), libMesh::Elem::n_nodes(), libMesh::MeshBase::n_nodes(), libMesh::Elem::n_second_order_adjacent_vertices(), libMesh::Elem::n_vertices(), libMesh::Elem::neighbor(), libMesh::MeshBase::node(), libMesh::Elem::parent(), libMesh::Elem::point(), std::pow(), libMesh::Real, libMesh::Elem::second_order_adjacent_vertex(), side, libMesh::TypeVector< T >::size(), and libMesh::MeshTools::weight().

{
  libmesh_assert_equal_to (mesh.mesh_dimension(), 2);

  /*
   * find the boundary nodes
   */
  std::vector<bool>  on_boundary;
  MeshTools::find_boundary_nodes(mesh, on_boundary);

  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
             */
            MeshBase::element_iterator       el  = mesh.level_elements_begin(refinement_level);
            const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level);

            for (; el != end; ++el)
              {
                /*
                 * Constant handle for the element
                 */
                const Elem* elem = *el;

                /*
                 * 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 (unsigned int s=0; s<elem->n_neighbors(); s++)
                      {
                        /*
                         * 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(s) != NULL) &&
                            (elem->id() > elem->neighbor(s)->id()) )
                          {
                            AutoPtr<Elem> side(elem->build_side(s));

                            Node* node0 = side->get_node(0);
                            Node* node1 = side->get_node(1);

                            Real node_weight = 1.;
                            // calculate the weight of the nodes
                            if (power > 0)
                              {
                                Point diff = (*node0)-(*node1);
                                node_weight = std::pow( diff.size(), 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
                     */
                    for (unsigned int c=0; c < parent->n_children(); c++)
                      {
                        if (parent->child(c) == elem)
                          {
                            /*
                             *loop over the childs (that is, the current elements) nodes
                             */
                            for (unsigned int nc=0; nc < elem->n_nodes(); nc++)
                              {
                                /*
                                 * the new position of the node
                                 */
                                Point point;
                                for (unsigned int n=0; n<parent->n_nodes(); n++)
                                  {
                                    /*
                                     * The value from the embedding matrix
                                     */
                                    const float 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->get_node(nc)->id();
                                new_positions[id] = point;
                                weight[id] = 1.;
                              }

                          } // if parent->child == elem
                      } // for parent->n_children
                  } // if element refinement_level
#endif // #ifdef LIBMESH_ENABLE_AMR

              } // element loop

            /*
             * finally reposition the vertex nodes
             */
            for (unsigned int nid=0; nid<mesh.n_nodes(); ++nid)
              if (!on_boundary[nid] && weight[nid] > 0.)
                mesh.node(nid) = new_positions[nid]/weight[nid];
          }

          {
            /*
             * Now handle the additional second_order nodes by calculating
             * their position based on the vertex postitions
             * we do a second loop over the level elements
             */
            MeshBase::element_iterator       el  = mesh.level_elements_begin(refinement_level);
            const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level);

            for (; el != end; ++el)
              {
                /*
                 * Constant handle for the element
                 */
                const Elem* elem = *el;
                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->get_node(n)->id();
                    mesh.node(id) = point/n_adjacent_vertices;
                  }
              }
          }

        } // refinement_level loop

    } // end iteration
}

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

References libMesh::MeshBase::nodes_begin(), and libMesh::MeshBase::nodes_end().

{
  const Point p(xt, yt, zt);

  const MeshBase::node_iterator nd_end = mesh.nodes_end();

  for (MeshBase::node_iterator nd = mesh.nodes_begin();
       nd != nd_end; ++nd)
    **nd += p;
}