Poster of Linux kernelThe best gift for a Linux geek
MeshTools::Modification

MeshTools::Modification

Section: C Library Functions (3) Updated: Thu Apr 7 2011
Local index Up
 

NAME

MeshTools::Modification -  

SYNOPSIS


 

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 short int old_id, const short int new_id)
 

Detailed Description

Tools for Mesh modification.

Author:

Benjamin S. Kirk

Date:

2004

Version:

Revision:

3391

 

Function Documentation

 

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

References MeshBase::add_elem(), MeshBase::add_point(), MeshBase::boundary_info, MeshBase::delete_elem(), MeshBase::elements_begin(), MeshBase::elements_end(), Utility::enum_to_string< ElemType >(), BoundaryInfo::invalid_id, MeshBase::n_active_elem(), MeshBase::node(), MeshBase::prepare_for_use(), DofObject::processor_id(), libMeshEnums::QUAD4, libMeshEnums::QUAD8, libMeshEnums::QUAD9, Elem::set_node(), Elem::subdomain_id(), libMeshEnums::TRI3, and libMeshEnums::TRI6.

{
  // 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*mesh.n_active_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<short int> new_bndry_ids;
  
  // Iterate over the elements, splitting QUADS into
  // pairs of conforming triangles.
  {
    MeshBase::element_iterator       el  = mesh.elements_begin();
    const MeshBase::element_iterator end = mesh.elements_end(); 

    for (; el!=end; ++el)
      {
        const ElemType etype = (*el)->type();

        // all_tri currently only works on coarse meshes
        libmesh_assert ((*el)->parent() == NULL);

        // 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 (((*el)->point(0) - (*el)->point(2)).size() <
                  ((*el)->point(1) - (*el)->point(3)).size())
                {             
                  tri0->set_node(0) = (*el)->get_node(0);
                  tri0->set_node(1) = (*el)->get_node(1);
                  tri0->set_node(2) = (*el)->get_node(2);
                
                  tri1->set_node(0) = (*el)->get_node(0);
                  tri1->set_node(1) = (*el)->get_node(2);
                  tri1->set_node(2) = (*el)->get_node(3);
                }
            
              else
                {
                  edge_swap=true;
                  
                  tri0->set_node(0) = (*el)->get_node(0);
                  tri0->set_node(1) = (*el)->get_node(1);
                  tri0->set_node(2) = (*el)->get_node(3);
                
                  tri1->set_node(0) = (*el)->get_node(1);
                  tri1->set_node(1) = (*el)->get_node(2);
                  tri1->set_node(2) = (*el)->get_node(3);
                }

              
              break;
            }
      
          case QUAD8:
            {
              split_elem =  true;
              
              tri0 = new Tri6;
              tri1 = new Tri6;
          
              Node* new_node = mesh.add_point((mesh.node((*el)->node(0)) +
                                               mesh.node((*el)->node(1)) +
                                               mesh.node((*el)->node(2)) +
                                               mesh.node((*el)->node(3)))*.25
                                               );
          
              // Check for possible edge swap
              if (((*el)->point(0) - (*el)->point(2)).size() <
                  ((*el)->point(1) - (*el)->point(3)).size())
                {             
                  tri0->set_node(0) = (*el)->get_node(0);
                  tri0->set_node(1) = (*el)->get_node(1);
                  tri0->set_node(2) = (*el)->get_node(2);
                  tri0->set_node(3) = (*el)->get_node(4);
                  tri0->set_node(4) = (*el)->get_node(5);
                  tri0->set_node(5) = new_node;
              
                  tri1->set_node(0) = (*el)->get_node(0);
                  tri1->set_node(1) = (*el)->get_node(2);
                  tri1->set_node(2) = (*el)->get_node(3);
                  tri1->set_node(3) = new_node;
                  tri1->set_node(4) = (*el)->get_node(6);
                  tri1->set_node(5) = (*el)->get_node(7);

                }
          
              else
                {
                  edge_swap=true;
                  
                  tri0->set_node(0) = (*el)->get_node(3);
                  tri0->set_node(1) = (*el)->get_node(0);
                  tri0->set_node(2) = (*el)->get_node(1);
                  tri0->set_node(3) = (*el)->get_node(7);
                  tri0->set_node(4) = (*el)->get_node(4);
                  tri0->set_node(5) = new_node;
              
                  tri1->set_node(0) = (*el)->get_node(1);
                  tri1->set_node(1) = (*el)->get_node(2);
                  tri1->set_node(2) = (*el)->get_node(3);
                  tri1->set_node(3) = (*el)->get_node(5);
                  tri1->set_node(4) = (*el)->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 (((*el)->point(0) - (*el)->point(2)).size() <
                  ((*el)->point(1) - (*el)->point(3)).size())
                {             
                  tri0->set_node(0) = (*el)->get_node(0);
                  tri0->set_node(1) = (*el)->get_node(1);
                  tri0->set_node(2) = (*el)->get_node(2);
                  tri0->set_node(3) = (*el)->get_node(4);
                  tri0->set_node(4) = (*el)->get_node(5);
                  tri0->set_node(5) = (*el)->get_node(8);
              
                  tri1->set_node(0) = (*el)->get_node(0);
                  tri1->set_node(1) = (*el)->get_node(2);
                  tri1->set_node(2) = (*el)->get_node(3);
                  tri1->set_node(3) = (*el)->get_node(8);
                  tri1->set_node(4) = (*el)->get_node(6);
                  tri1->set_node(5) = (*el)->get_node(7);
                }

              else
                {
                  edge_swap=true;
                  
                  tri0->set_node(0) = (*el)->get_node(0);
                  tri0->set_node(1) = (*el)->get_node(1);
                  tri0->set_node(2) = (*el)->get_node(3);
                  tri0->set_node(3) = (*el)->get_node(4);
                  tri0->set_node(4) = (*el)->get_node(8);
                  tri0->set_node(5) = (*el)->get_node(7);
              
                  tri1->set_node(0) = (*el)->get_node(1);
                  tri1->set_node(1) = (*el)->get_node(2);
                  tri1->set_node(2) = (*el)->get_node(3);
                  tri1->set_node(3) = (*el)->get_node(5);
                  tri1->set_node(4) = (*el)->get_node(6);
                  tri1->set_node(5) = (*el)->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:
            {
              std::cerr << '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() = (*el)->processor_id();
            tri0->subdomain_id() = (*el)->subdomain_id();
            tri1->processor_id() = (*el)->processor_id();
            tri1->subdomain_id() = (*el)->subdomain_id();
           
            if (mesh_has_boundary_data)
              {
                for (unsigned int sn=0; sn<(*el)->n_sides(); ++sn)
                  {
                    short int b_id = mesh.boundary_info->boundary_id(*el, sn);

                    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:
                                {
                                  std::cerr << '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:
                                {
                                  std::cerr << '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 sides

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

              } // end if (mesh_has_boundary_data)
            
            // 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(*el);
          } // end if (split_elem)
      } // End for loop over elements
  } 

  
    // 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 (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 (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 (new_bndry_elements.size() == new_bndry_sides.size());
      libmesh_assert (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();

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

void MeshTools::Modification::change_boundary_id (MeshBase &mesh, const short intold_id, const short intnew_id)any boundary ids that are currently old_id, changes them to new_id

Definition at line 1279 of file mesh_modification.C.

References MeshBase::boundary_info, MeshBase::level_elements_begin(), MeshBase::level_elements_end(), and 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)
        if (mesh.boundary_info->boundary_id(elem, s) == old_id)
          {
            mesh.boundary_info->remove_side(elem, s, old_id);
            mesh.boundary_info->add_side(elem, s, new_id);
          }
    }
}
 

void MeshTools::Modification::distort (MeshBase &mesh, const Realfactor, const boolperturb_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 43 of file mesh_modification.C.

References MeshBase::active_elements_begin(), MeshBase::active_elements_end(), MeshTools::find_boundary_nodes(), MeshBase::mesh_dimension(), std::min(), MeshBase::n_elem(), MeshBase::n_nodes(), MeshBase::node(), and 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(), 1.e20);

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

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


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

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

References MeshBase::active_elements_begin(), MeshBase::active_elements_end(), MeshBase::add_elem(), MeshBase::boundary_info, Elem::build(), MeshBase::delete_elem(), MeshBase::elements_begin(), MeshBase::elements_end(), Elem::get_node(), BoundaryInfo::invalid_id, MeshBase::n_active_elem(), Elem::n_nodes(), Elem::n_sides(), Elem::neighbor(), MeshBase::prepare_for_use(), DofObject::processor_id(), Elem::set_node(), Elem::subdomain_id(), and 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<short int>          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 point to nodes in the original mesh)
        // The copy's element ID will be set upon adding it to the 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();

        // This element could have boundary info as well.  We need
        // to save the (elem, side, bc_id) triples
        for (unsigned int s=0; s<elem->n_sides(); s++)
            if (elem->neighbor(s) == NULL)
              {
                short int bc_id = mesh.boundary_info->boundary_id (elem,s);

                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 (saved_boundary_elements.size() == saved_bc_ids.size());
    libmesh_assert (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)
      mesh.add_elem(*it);
  }

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

void MeshTools::Modification::rotate (MeshBase &mesh, const Realphi, const Realtheta = 0., const Realpsi = 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 159 of file mesh_modification.C.

References MeshBase::mesh_dimension(), MeshBase::n_nodes(), MeshBase::node(), and libMesh::pi.

{
  libmesh_assert (mesh.mesh_dimension() != 1);

  const Real pi = std::acos(-1.);
  const Real  p = -phi/180.*pi;
  const Real  t = -theta/180.*pi;
  const Real  s = -psi/180.*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.
  for (unsigned int n=0; n<mesh.n_nodes(); n++)
    {
      const Point p = mesh.node(n);
      const Real  x = p(0);
      const Real  y = p(1);
      const Real  z = p(2);
      mesh.node(n) = 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 MeshTools::Modification::scale (MeshBase &mesh, const Realxs, const Realys = 0., const Realzs = 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 191 of file mesh_modification.C.

References MeshBase::n_nodes(), MeshBase::node(), and MeshBase::spatial_dimension().

Referenced by DenseVector< T >::operator*=(), and DenseMatrix< T >::operator*=().

{
  const Real x_scale = xs;
  Real y_scale       = ys;
  Real z_scale       = zs;
  
  if (ys == 0.)
    {
      libmesh_assert (zs == 0.);

      y_scale = z_scale = x_scale;
    }

  // Scale the x coordinate in all dimensions
  for (unsigned int n=0; n<mesh.n_nodes(); n++)
    mesh.node(n)(0) *= x_scale;


  // Only scale the y coordinate in 2 and 3D
  if (mesh.spatial_dimension() > 1)
    {

      for (unsigned int n=0; n<mesh.n_nodes(); n++)
        mesh.node(n)(1) *= y_scale;

      // Only scale the z coordinate in 3D
      if (mesh.spatial_dimension() == 3)
        {
          for (unsigned int n=0; n<mesh.n_nodes(); n++)
            mesh.node(n)(2) *= z_scale;
        }
    }
}
 

void MeshTools::Modification::smooth (MeshBase &mesh, unsigned intn_iterations, Realpower)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 992 of file mesh_modification.C.

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

{
  libmesh_assert (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 unsigned int 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 unsigned int 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 unsigned int id = elem->get_node(n)->id();
                    mesh.node(id) = point/n_adjacent_vertices;
                  }
              }
          }
      
        } // refinement_level loop

    } // end iteration
}
 

void MeshTools::Modification::translate (MeshBase &mesh, const Realxt = 0., const Realyt = 0., const Realzt = 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 125 of file mesh_modification.C.

References MeshBase::n_nodes(), and MeshBase::node().

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

  for (unsigned int n=0; n<mesh.n_nodes(); n++)
    mesh.node(n) += p;
}
 

Author

Generated automatically by Doxygen for libMesh from the source code.


 

Index

NAME
SYNOPSIS
Functions
Detailed Description
Function Documentation
void 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.
void MeshTools::Modification::change_boundary_id (MeshBase &mesh, const short intold_id, const short intnew_id)any boundary ids that are currently old_id, changes them to new_id
void MeshTools::Modification::distort (MeshBase &mesh, const Realfactor, const boolperturb_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.
void 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.
void MeshTools::Modification::rotate (MeshBase &mesh, const Realphi, const Realtheta = 0., const Realpsi = 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)
void MeshTools::Modification::scale (MeshBase &mesh, const Realxs, const Realys = 0., const Realzs = 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.
void MeshTools::Modification::smooth (MeshBase &mesh, unsigned intn_iterations, Realpower)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.
void MeshTools::Modification::translate (MeshBase &mesh, const Realxt = 0., const Realyt = 0., const Realzt = 0.)Translates the mesh. The grid points are translated in the x direction by xt, in the y direction by yt, etc...
Author

This document was created by man2html, using the manual pages.
Time: 21:50:11 GMT, April 16, 2011