libMesh::LaplaceMeshSmoother Class Reference

#include <mesh_smoother_laplace.h>

Inheritance diagram for libMesh::LaplaceMeshSmoother:

[legend]

List of all members.

Public Member Functions
	LaplaceMeshSmoother (UnstructuredMesh &mesh)
virtual	~LaplaceMeshSmoother ()
virtual void	smooth ()
void	smooth (unsigned int n_iterations)
void	init ()
void	print_graph (std::ostream &out=libMesh::out) const
Protected Attributes
UnstructuredMesh &	_mesh
Private Member Functions
void	allgather_graph ()
Private Attributes
bool	_initialized
std::vector< std::vector < dof_id_type > >	_graph

---

Detailed Description

This class defines the data structures necessary for Laplace smoothing. Note that this is a simple averaging smoother, which does NOT guarantee that points will be smoothed to valid locations, e.g. locations inside the boundary! This aspect could use work.

Author:

John W. Peterson

Date:

2002-2007

Definition at line 52 of file mesh_smoother_laplace.h.

---

Constructor & Destructor Documentation

libMesh::LaplaceMeshSmoother::LaplaceMeshSmoother

(

UnstructuredMesh &

mesh

)

[explicit]

Constructor. Sets the constant mesh reference in the protected data section of the class.

Definition at line 36 of file mesh_smoother_laplace.C.

    : MeshSmoother(mesh),
      _initialized(false)
  {
  }

virtual libMesh::LaplaceMeshSmoother::~LaplaceMeshSmoother

(

)

[inline, virtual]

Destructor.

Definition at line 65 of file mesh_smoother_laplace.h.

{}

---

Member Function Documentation

void libMesh::LaplaceMeshSmoother::allgather_graph

(

)

[private]

This function allgather's the (local) graph after it is computed on each processor by the init() function.

Definition at line 313 of file mesh_smoother_laplace.C.

References _graph, libMesh::MeshSmoother::_mesh, libMesh::Parallel::Communicator::allgather(), libMesh::CommWorld, libMesh::MeshBase::n_nodes(), and libMesh::n_processors().

Referenced by init().

  {
    // The graph data structure is not well-suited for parallel communication,
    // so copy the graph into a single vector defined by:
    // NA A_0 A_1 ... A_{NA} | NB B_0 B_1 ... B_{NB} | NC C_0 C_1 ... C_{NC}
    // where:
    // * NA is the number of graph connections for node A
    // * A_0, A_1, etc. are the IDs connected to node A
    std::vector<std::size_t> flat_graph;

    // Reserve at least enough space for each node to have zero entries
    flat_graph.reserve(_graph.size());

    for (std::size_t i=0; i<_graph.size(); ++i)
      {
        // First push back the number of entries for this node
        flat_graph.push_back (_graph[i].size());

        // Then push back all the IDs
        for (std::size_t j=0; j<_graph[i].size(); ++j)
          flat_graph.push_back(_graph[i][j]);
      }

    // // A copy of the flat graph (for printing only, delete me later)
    // std::vector<unsigned> copy_of_flat_graph(flat_graph);

    // Use the allgather routine to combine all the flat graphs on all processors
    CommWorld.allgather(flat_graph);

    // Now reconstruct _graph from the allgathered flat_graph.

    // // (Delete me later, the copy is just for printing purposes.)
    // std::vector<std::vector<unsigned > > copy_of_graph(_graph);

    // Make sure the old graph is cleared out
    _graph.clear();
    _graph.resize(_mesh.n_nodes());

    // Our current position in the allgather'd flat_graph
    std::size_t cursor=0;

    // There are n_nodes * n_processors vectors to read in total
    for (processor_id_type p=0; p<libMesh::n_processors(); ++p)
      for (dof_id_type node_ctr=0; node_ctr<_mesh.n_nodes(); ++node_ctr)
        {
          // Read the number of entries for this node, move cursor
          std::size_t n_entries = flat_graph[cursor++];

          // Reserve space for that many more entries, then push back
          _graph[node_ctr].reserve(_graph[node_ctr].size() + n_entries);

          // Read all graph connections for this node, move the cursor each time
          // Note: there might be zero entries but that's fine
          for (std::size_t i=0; i<n_entries; ++i)
            _graph[node_ctr].push_back(flat_graph[cursor++]);
        }


//    // Print local graph to uniquely named file (debugging)
//    {
//      // Generate unique filename for this processor
//      std::ostringstream oss;
//      oss << "graph_filename_" << libMesh::processor_id() << ".txt";
//      std::ofstream graph_stream(oss.str().c_str());
//
//      // Print the local non-flat graph
//      std::swap(_graph, copy_of_graph);
//      print_graph(graph_stream);
//
//      // Print the (local) flat graph for verification
//      for (std::size_t i=0; i<copy_of_flat_graph.size(); ++i)
//      graph_stream << copy_of_flat_graph[i] << " ";
//      graph_stream << "\n";
//
//      // Print the allgather'd grap for verification
//      for (std::size_t i=0; i<flat_graph.size(); ++i)
//      graph_stream << flat_graph[i] << " ";
//      graph_stream << "\n";
//
//      // Print the global non-flat graph
//      std::swap(_graph, copy_of_graph);
//      print_graph(graph_stream);
//    }
  } // allgather_graph()

void libMesh::LaplaceMeshSmoother::init

(

)

Initialization for the Laplace smoothing routine is basically identical to building an "L-graph" which is expensive. It's provided separately from the constructor since you may or may not want to build the L-graph on construction.

Definition at line 182 of file mesh_smoother_laplace.C.

References _graph, _initialized, libMesh::MeshSmoother::_mesh, libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), allgather_graph(), libMesh::Elem::build_side(), end, libMesh::err, libMesh::DofObject::id(), libMesh::MeshBase::mesh_dimension(), libMesh::Elem::n_neighbors(), libMesh::MeshBase::n_nodes(), libMesh::Elem::neighbor(), and side.

Referenced by smooth().

{
  switch (_mesh.mesh_dimension())
    {

      // TODO:[BSK] Fix this to work for refined meshes...  I think
      // the implementation was done quickly for Damien, who did not have
      // refined grids.  Fix it here and in the original Mesh member.

    case 2: // Stolen directly from build_L_graph in mesh_base.C
      {
        // Initialize space in the graph.  It is n_nodes long.  Each
        // node may be connected to an arbitrary number of other nodes
        // via edges.
        _graph.resize(_mesh.n_nodes());

        MeshBase::element_iterator       el  = _mesh.active_local_elements_begin();
        const MeshBase::element_iterator end = _mesh.active_local_elements_end();

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

            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));
                    _graph[side->node(0)].push_back(side->node(1));
                    _graph[side->node(1)].push_back(side->node(0));
                  }
              }
          }
        _initialized = true;
        break;
      } // case 2

    case 3: // Stolen blatantly from build_L_graph in mesh_base.C
      {
        // Initialize space in the graph.
        _graph.resize(_mesh.n_nodes());

        MeshBase::element_iterator       el  = _mesh.active_local_elements_begin();
        const MeshBase::element_iterator end = _mesh.active_local_elements_end();

        for (; el != end; ++el)
          {
            // Shortcut notation for simplicity
            const Elem* elem = *el;

            for (unsigned int f=0; f<elem->n_neighbors(); f++) // Loop over faces
              if ((elem->neighbor(f) == NULL) ||
                  (elem->id() > elem->neighbor(f)->id()))
                {
                  // We need a full (i.e. non-proxy) element for the face, since we will
                  // be looking at its sides as well!
                  AutoPtr<Elem> face = elem->build_side(f, /*proxy=*/false);

                  for (unsigned int s=0; s<face->n_neighbors(); s++) // Loop over face's edges
                    {
                      // Here we can use a proxy
                      AutoPtr<Elem> side = face->build_side(s);

                      // At this point, we just insert the node numbers
                      // again.  At the end we'll call sort and unique
                      // to make sure there are no duplicates
                      _graph[side->node(0)].push_back(side->node(1));
                      _graph[side->node(1)].push_back(side->node(0));
                    }
                }
          }

        _initialized = true;
        break;
      } // case 3

    default:
      {
        libMesh::err << "At this time it is not possible "
                      << "to smooth a dimension "
                      << _mesh.mesh_dimension()
                      << "mesh.  Aborting..."
                      << std::endl;
        libmesh_error();
      }
    }

  // Done building graph from local elements.  Let's now allgather the
  // graph so that it is available on all processors for the actual
  // smoothing operation?
  this->allgather_graph();

  // In 3D, it's possible for > 2 processor partitions to meet
  // at a single edge, while in 2D only 2 processor partitions
  // share an edge.  Therefore the allgather'd graph in 3D may
  // now have duplicate entries and we need to remove them so
  // they don't foul up the averaging algorithm employed by the
  // Laplace smoother.
  for (unsigned i=0; i<_graph.size(); ++i)
    {
      // The std::unique algorithm removes duplicate *consecutive* elements from a range,
      // so it only makes sense to call it on a sorted range...
      std::sort(_graph[i].begin(), _graph[i].end());
      _graph[i].erase(std::unique(_graph[i].begin(), _graph[i].end()), _graph[i].end());
    }

} // init()

void libMesh::LaplaceMeshSmoother::print_graph

(

std::ostream &

out = libMesh::out

)

const

Mainly for debugging, this function will print out the connectivity graph which has been created.

Definition at line 299 of file mesh_smoother_laplace.C.

References _graph, and end.

{
  for (unsigned int i=0; i<_graph.size(); ++i)
    {
      out_stream << i << ": ";
      std::copy(_graph[i].begin(),
                _graph[i].end(),
                std::ostream_iterator<unsigned>(out_stream, " "));
      out_stream << std::endl;
    }
}

void libMesh::LaplaceMeshSmoother::smooth

(

unsigned int

n_iterations

)

The actual smoothing function, gets called whenever the user specifies an actual number of smoothing iterations.

Definition at line 45 of file mesh_smoother_laplace.C.

References _graph, _initialized, libMesh::MeshSmoother::_mesh, libMesh::MeshBase::active_elements_begin(), libMesh::MeshBase::active_elements_end(), libMesh::TypeVector< T >::add(), end, libMesh::err, libMesh::MeshTools::find_boundary_nodes(), libMesh::DofObject::id(), init(), libMesh::MeshBase::local_nodes_begin(), libMesh::MeshBase::local_nodes_end(), libMesh::Elem::n_nodes(), libMesh::MeshBase::n_nodes(), libMesh::Elem::n_second_order_adjacent_vertices(), libMesh::Elem::n_vertices(), libMesh::Elem::node(), libMesh::MeshBase::node(), libMesh::MeshBase::node_ptr(), libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::MeshBase::point(), libMesh::processor_id(), libMesh::Real, libMesh::Elem::second_order_adjacent_vertex(), and libMesh::Parallel::sync_dofobject_data_by_id().

{
  if (!_initialized)
    this->init();

  // Don't smooth the nodes on the boundary...
  // this would change the mesh geometry which
  // is probably not something we want!
  std::vector<bool> on_boundary;
  MeshTools::find_boundary_nodes(_mesh, on_boundary);

  // Ensure that the find_boundary_nodes() function returned a properly-sized vector
  if (on_boundary.size() != _mesh.n_nodes())
    {
      libMesh::err << "MeshTools::find_boundary_nodes() returned incorrect length vector!" << std::endl;
      libmesh_error();
    }

  // We can only update the nodes after all new positions were
  // determined. We store the new positions here
  std::vector<Point> new_positions;

  for (unsigned int n=0; n<n_iterations; n++)
    {
      new_positions.resize(_mesh.n_nodes());

      {
        MeshBase::node_iterator       it     = _mesh.local_nodes_begin();
        const MeshBase::node_iterator it_end = _mesh.local_nodes_end();
        for (; it != it_end; ++it)
          {
            Node* node = *it;

            if (node == NULL)
              {
                libMesh::err << "[" << libMesh::processor_id() << "]: Node iterator returned NULL pointer." << std::endl;
                libmesh_error();
              }

            // leave the boundary intact
            // Only relocate the nodes which are vertices of an element
            // All other entries of _graph (the secondary nodes) are empty
            if (!on_boundary[node->id()] && (_graph[node->id()].size() > 0))
              {
                Point avg_position(0.,0.,0.);

                for (unsigned j=0; j<_graph[node->id()].size(); ++j)
                  {
                    // Will these nodal positions always be available
                    // or will they refer to remote nodes?  To be
                    // careful, we grab a pointer and test it against
                    // NULL.
                    Node* connected_node = _mesh.node_ptr(_graph[node->id()][j]);

                    if (connected_node == NULL)
                      {
                        libMesh::err << "Error! Libmesh returned NULL pointer for node " << _graph[connected_node->id()][j] << std::endl;
                        libmesh_error();
                      }

                    avg_position.add( *connected_node );
                  } // end for(j)

                // Compute the average, store in the new_positions vector
                new_positions[node->id()] = avg_position / static_cast<Real>(_graph[node->id()].size());
              } // end if
          } // end for
      } // end scope


      // now update the node positions (local node positions only)
      {
        MeshBase::node_iterator it           = _mesh.local_nodes_begin();
        const MeshBase::node_iterator it_end = _mesh.local_nodes_end();
        for (; it != it_end; ++it)
          {
            Node* node = *it;

            if (!on_boundary[node->id()] && (_graph[node->id()].size() > 0))
              {
                // Should call Point::op=
                // libMesh::out << "Setting node id " << node->id() << " to position " << new_positions[node->id()];
                _mesh.node(node->id()) = new_positions[node->id()];
              }
          } // end for
      } // end scope

      // Now the nodes which are ghosts on this processor may have been moved on
      // the processors which own them.  So we need to synchronize with our neighbors
      // and get the most up-to-date positions for the ghosts.
      SyncNodalPositions sync_object(_mesh);
      Parallel::sync_dofobject_data_by_id(_mesh.nodes_begin(), _mesh.nodes_end(), sync_object);

    } // end for n_iterations

  // finally adjust the second order nodes (those located between vertices)
  // these nodes will be located between their adjacent nodes
  // do this element-wise
  MeshBase::element_iterator       el  = _mesh.active_elements_begin();
  const MeshBase::element_iterator end = _mesh.active_elements_end();

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

      // get the second order nodes (son)
      // their element indices start at n_vertices and go to n_nodes
      const unsigned int son_begin = elem->n_vertices();
      const unsigned int son_end   = elem->n_nodes();

      // loop over all second order nodes (son)
      for (unsigned int son=son_begin; son<son_end; son++)
        {
          // Don't smooth second-order nodes which are on the boundary
          if (!on_boundary[elem->node(son)])
            {
              const unsigned int n_adjacent_vertices =
                elem->n_second_order_adjacent_vertices(son);

              // calculate the new position which is the average of the
              // position of the adjacent vertices
              Point avg_position(0,0,0);
              for (unsigned int v=0; v<n_adjacent_vertices; v++)
                avg_position +=
                  _mesh.point( elem->node( elem->second_order_adjacent_vertex(son,v) ) );

              _mesh.node(elem->node(son)) = avg_position / n_adjacent_vertices;
            }
        }
    }
}

virtual void libMesh::LaplaceMeshSmoother::smooth

(

)

[inline, virtual]

Redefinition of the smooth function from the base class. All this does is call the smooth function in this class which takes an int, using a default value of 1.

Implements libMesh::MeshSmoother.

Definition at line 73 of file mesh_smoother_laplace.h.

References smooth().

Referenced by smooth().

{ this->smooth(1); }

---

Member Data Documentation

std::vector<std::vector<dof_id_type> > libMesh::LaplaceMeshSmoother::_graph [private]

Data structure for holding the L-graph

Definition at line 112 of file mesh_smoother_laplace.h.

Referenced by allgather_graph(), init(), print_graph(), and smooth().

bool libMesh::LaplaceMeshSmoother::_initialized [private]

True if the L-graph has been created, false otherwise.

Definition at line 107 of file mesh_smoother_laplace.h.

Referenced by init(), and smooth().

UnstructuredMesh& libMesh::MeshSmoother::_mesh [protected, inherited]

Definition at line 71 of file mesh_smoother.h.

Referenced by libMesh::VariationalMeshSmoother::adjust_adapt_data(), allgather_graph(), init(), libMesh::VariationalMeshSmoother::readgr(), libMesh::VariationalMeshSmoother::smooth(), smooth(), and libMesh::VariationalMeshSmoother::writegr().

---

The documentation for this class was generated from the following files:

• mesh_smoother_laplace.h
• mesh_smoother_laplace.C