libMesh::LaplaceMeshSmoother Class Reference

#include <mesh_smoother_laplace.h>

Inheritance diagram for libMesh::LaplaceMeshSmoother:

Public Member Functions

 LaplaceMeshSmoother (UnstructuredMesh &mesh)
 
virtual ~LaplaceMeshSmoother ()
 
virtual void smooth () libmesh_override
 
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
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 44 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.

37  : MeshSmoother(mesh),
38  _initialized(false)
39 {
40 }
MeshBase & mesh
MeshSmoother(UnstructuredMesh &mesh)
Definition: mesh_smoother.h:46
virtual libMesh::LaplaceMeshSmoother::~LaplaceMeshSmoother ( )
inlinevirtual

Destructor.

Definition at line 57 of file mesh_smoother_laplace.h.

57 {}

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 296 of file mesh_smoother_laplace.C.

References _graph, libMesh::MeshSmoother::_mesh, libMesh::Parallel::Communicator::allgather(), libMesh::ParallelObject::comm(), libMesh::MeshBase::max_node_id(), and libMesh::ParallelObject::n_processors().

Referenced by init(), and smooth().

297 {
298  // The graph data structure is not well-suited for parallel communication,
299  // so copy the graph into a single vector defined by:
300  // NA A_0 A_1 ... A_{NA} | NB B_0 B_1 ... B_{NB} | NC C_0 C_1 ... C_{NC}
301  // where:
302  // * NA is the number of graph connections for node A
303  // * A_0, A_1, etc. are the IDs connected to node A
304  std::vector<dof_id_type> flat_graph;
305 
306  // Reserve at least enough space for each node to have zero entries
307  flat_graph.reserve(_graph.size());
308 
309  for (std::size_t i=0; i<_graph.size(); ++i)
310  {
311  // First push back the number of entries for this node
312  flat_graph.push_back (cast_int<dof_id_type>(_graph[i].size()));
313 
314  // Then push back all the IDs
315  for (std::size_t j=0; j<_graph[i].size(); ++j)
316  flat_graph.push_back(_graph[i][j]);
317  }
318 
319  // // A copy of the flat graph (for printing only, delete me later)
320  // std::vector<unsigned> copy_of_flat_graph(flat_graph);
321 
322  // Use the allgather routine to combine all the flat graphs on all processors
323  _mesh.comm().allgather(flat_graph);
324 
325  // Now reconstruct _graph from the allgathered flat_graph.
326 
327  // // (Delete me later, the copy is just for printing purposes.)
328  // std::vector<std::vector<unsigned > > copy_of_graph(_graph);
329 
330  // Make sure the old graph is cleared out
331  _graph.clear();
332  _graph.resize(_mesh.max_node_id());
333 
334  // Our current position in the allgather'd flat_graph
335  std::size_t cursor=0;
336 
337  // There are max_node_id * n_processors entries to read in total
338  for (processor_id_type p=0; p<_mesh.n_processors(); ++p)
339  for (dof_id_type node_ctr=0; node_ctr<_mesh.max_node_id(); ++node_ctr)
340  {
341  // Read the number of entries for this node, move cursor
342  std::size_t n_entries = flat_graph[cursor++];
343 
344  // Reserve space for that many more entries, then push back
345  _graph[node_ctr].reserve(_graph[node_ctr].size() + n_entries);
346 
347  // Read all graph connections for this node, move the cursor each time
348  // Note: there might be zero entries but that's fine
349  for (std::size_t i=0; i<n_entries; ++i)
350  _graph[node_ctr].push_back(flat_graph[cursor++]);
351  }
352 
353 
354  // // Print local graph to uniquely named file (debugging)
355  // {
356  // // Generate unique filename for this processor
357  // std::ostringstream oss;
358  // oss << "graph_filename_" << _mesh.processor_id() << ".txt";
359  // std::ofstream graph_stream(oss.str().c_str());
360  //
361  // // Print the local non-flat graph
362  // std::swap(_graph, copy_of_graph);
363  // print_graph(graph_stream);
364  //
365  // // Print the (local) flat graph for verification
366  // for (std::size_t i=0; i<copy_of_flat_graph.size(); ++i)
367  //graph_stream << copy_of_flat_graph[i] << " ";
368  // graph_stream << "\n";
369  //
370  // // Print the allgather'd grap for verification
371  // for (std::size_t i=0; i<flat_graph.size(); ++i)
372  //graph_stream << flat_graph[i] << " ";
373  // graph_stream << "\n";
374  //
375  // // Print the global non-flat graph
376  // std::swap(_graph, copy_of_graph);
377  // print_graph(graph_stream);
378  // }
379 } // allgather_graph()
std::vector< std::vector< dof_id_type > > _graph
processor_id_type n_processors() const
virtual dof_id_type max_node_id() const =0
uint8_t processor_id_type
Definition: id_types.h:99
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
const Parallel::Communicator & comm() const
uint8_t dof_id_type
Definition: id_types.h:64
void allgather(const T &send, std::vector< T > &recv) const
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 172 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_ptr(), end, libMesh::DofObject::id(), libmesh_nullptr, libMesh::MeshBase::max_node_id(), libMesh::MeshBase::mesh_dimension(), libMesh::Elem::n_neighbors(), libMesh::Elem::neighbor_ptr(), and side.

Referenced by smooth().

173 {
174  switch (_mesh.mesh_dimension())
175  {
176 
177  // TODO:[BSK] Fix this to work for refined meshes... I think
178  // the implementation was done quickly for Damien, who did not have
179  // refined grids. Fix it here and in the original Mesh member.
180 
181  case 2: // Stolen directly from build_L_graph in mesh_base.C
182  {
183  // Initialize space in the graph. It is indexed by node id.
184  // Each node may be connected to an arbitrary number of other
185  // nodes via edges.
186  _graph.resize(_mesh.max_node_id());
187 
188  MeshBase::element_iterator el = _mesh.active_local_elements_begin();
189  const MeshBase::element_iterator end = _mesh.active_local_elements_end();
190 
191  for (; el != end; ++el)
192  {
193  // Constant handle for the element
194  const Elem * elem = *el;
195 
196  for (unsigned int s=0; s<elem->n_neighbors(); s++)
197  {
198  // Only operate on sides which are on the
199  // boundary or for which the current element's
200  // id is greater than its neighbor's.
201  // Sides get only built once.
202  if ((elem->neighbor_ptr(s) == libmesh_nullptr) ||
203  (elem->id() > elem->neighbor_ptr(s)->id()))
204  {
205  UniquePtr<const Elem> side(elem->build_side_ptr(s));
206  _graph[side->node_id(0)].push_back(side->node_id(1));
207  _graph[side->node_id(1)].push_back(side->node_id(0));
208  }
209  }
210  }
211  _initialized = true;
212  break;
213  } // case 2
214 
215  case 3: // Stolen blatantly from build_L_graph in mesh_base.C
216  {
217  // Initialize space in the graph.
218  _graph.resize(_mesh.max_node_id());
219 
220  MeshBase::element_iterator el = _mesh.active_local_elements_begin();
221  const MeshBase::element_iterator end = _mesh.active_local_elements_end();
222 
223  for (; el != end; ++el)
224  {
225  // Shortcut notation for simplicity
226  const Elem * elem = *el;
227 
228  for (unsigned int f=0; f<elem->n_neighbors(); f++) // Loop over faces
229  if ((elem->neighbor_ptr(f) == libmesh_nullptr) ||
230  (elem->id() > elem->neighbor_ptr(f)->id()))
231  {
232  // We need a full (i.e. non-proxy) element for the face, since we will
233  // be looking at its sides as well!
234  UniquePtr<const Elem> face = elem->build_side_ptr(f, /*proxy=*/false);
235 
236  for (unsigned int s=0; s<face->n_neighbors(); s++) // Loop over face's edges
237  {
238  // Here we can use a proxy
239  UniquePtr<const Elem> side = face->build_side_ptr(s);
240 
241  // At this point, we just insert the node numbers
242  // again. At the end we'll call sort and unique
243  // to make sure there are no duplicates
244  _graph[side->node_id(0)].push_back(side->node_id(1));
245  _graph[side->node_id(1)].push_back(side->node_id(0));
246  }
247  }
248  }
249 
250  _initialized = true;
251  break;
252  } // case 3
253 
254  default:
255  libmesh_error_msg("At this time it is not possible to smooth a dimension " << _mesh.mesh_dimension() << "mesh. Aborting...");
256  }
257 
258  // Done building graph from local elements. Let's now allgather the
259  // graph so that it is available on all processors for the actual
260  // smoothing operation?
261  this->allgather_graph();
262 
263  // In 3D, it's possible for > 2 processor partitions to meet
264  // at a single edge, while in 2D only 2 processor partitions
265  // share an edge. Therefore the allgather'd graph in 3D may
266  // now have duplicate entries and we need to remove them so
267  // they don't foul up the averaging algorithm employed by the
268  // Laplace smoother.
269  for (std::size_t i=0; i<_graph.size(); ++i)
270  {
271  // The std::unique algorithm removes duplicate *consecutive* elements from a range,
272  // so it only makes sense to call it on a sorted range...
273  std::sort(_graph[i].begin(), _graph[i].end());
274  _graph[i].erase(std::unique(_graph[i].begin(), _graph[i].end()), _graph[i].end());
275  }
276 
277 } // init()
std::vector< std::vector< dof_id_type > > _graph
unsigned short int side
Definition: xdr_io.C:49
virtual dof_id_type max_node_id() const =0
const class libmesh_nullptr_t libmesh_nullptr
IterBase * end
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
virtual element_iterator active_local_elements_begin()=0
unsigned int mesh_dimension() const
Definition: mesh_base.C:147
virtual element_iterator active_local_elements_end()=0
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 282 of file mesh_smoother_laplace.C.

References _graph, and end.

Referenced by smooth().

283 {
284  for (std::size_t i=0; i<_graph.size(); ++i)
285  {
286  out_stream << i << ": ";
287  std::copy(_graph[i].begin(),
288  _graph[i].end(),
289  std::ostream_iterator<unsigned>(out_stream, " "));
290  out_stream << std::endl;
291  }
292 }
std::vector< std::vector< dof_id_type > > _graph
IterBase * end
virtual void libMesh::LaplaceMeshSmoother::smooth ( )
inlinevirtual

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 65 of file mesh_smoother_laplace.h.

References allgather_graph(), init(), libMesh::out, print_graph(), and smooth().

Referenced by libMesh::MeshTools::Generation::build_sphere(), smooth(), and libMesh::TriangleInterface::triangulate().

65 { this->smooth(1); }
virtual void smooth() libmesh_override
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(), libMesh::ParallelObject::comm(), end, libMesh::MeshTools::find_boundary_nodes(), libMesh::DofObject::id(), init(), libmesh_nullptr, libMesh::MeshBase::local_nodes_begin(), libMesh::MeshBase::local_nodes_end(), libMesh::MeshBase::max_node_id(), libMesh::Elem::n_nodes(), libMesh::Elem::n_second_order_adjacent_vertices(), libMesh::Elem::n_vertices(), libMesh::Elem::node_id(), libMesh::MeshBase::node_ref(), libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::MeshBase::point(), libMesh::ParallelObject::processor_id(), libMesh::Real, libMesh::Elem::second_order_adjacent_vertex(), and libMesh::Parallel::sync_dofobject_data_by_id().

46 {
47  if (!_initialized)
48  this->init();
49 
50  // Don't smooth the nodes on the boundary...
51  // this would change the mesh geometry which
52  // is probably not something we want!
53  std::vector<bool> on_boundary;
55 
56  // Ensure that the find_boundary_nodes() function returned a properly-sized vector
57  if (on_boundary.size() != _mesh.max_node_id())
58  libmesh_error_msg("MeshTools::find_boundary_nodes() returned incorrect length vector!");
59 
60  // We can only update the nodes after all new positions were
61  // determined. We store the new positions here
62  std::vector<Point> new_positions;
63 
64  for (unsigned int n=0; n<n_iterations; n++)
65  {
66  new_positions.resize(_mesh.max_node_id());
67 
68  {
69  MeshBase::node_iterator it = _mesh.local_nodes_begin();
70  const MeshBase::node_iterator it_end = _mesh.local_nodes_end();
71  for (; it != it_end; ++it)
72  {
73  Node * node = *it;
74 
75  if (node == libmesh_nullptr)
76  libmesh_error_msg("[" << _mesh.processor_id() << "]: Node iterator returned NULL pointer.");
77 
78  // leave the boundary intact
79  // Only relocate the nodes which are vertices of an element
80  // All other entries of _graph (the secondary nodes) are empty
81  if (!on_boundary[node->id()] && (_graph[node->id()].size() > 0))
82  {
83  Point avg_position(0.,0.,0.);
84 
85  for (std::size_t j=0; j<_graph[node->id()].size(); ++j)
86  {
87  // Will these nodal positions always be available
88  // or will they refer to remote nodes? This will
89  // fail an assertion in the latter case, which
90  // shouldn't occur if DistributedMesh is working
91  // correctly.
92  const Point & connected_node = _mesh.point(_graph[node->id()][j]);
93 
94  avg_position.add( connected_node );
95  } // end for (j)
96 
97  // Compute the average, store in the new_positions vector
98  new_positions[node->id()] = avg_position / static_cast<Real>(_graph[node->id()].size());
99  } // end if
100  } // end for
101  } // end scope
102 
103 
104  // now update the node positions (local node positions only)
105  {
106  MeshBase::node_iterator it = _mesh.local_nodes_begin();
107  const MeshBase::node_iterator it_end = _mesh.local_nodes_end();
108  for (; it != it_end; ++it)
109  {
110  Node * node = *it;
111 
112  if (!on_boundary[node->id()] && (_graph[node->id()].size() > 0))
113  {
114  // Should call Point::op=
115  // libMesh::out << "Setting node id " << node->id() << " to position " << new_positions[node->id()];
116  _mesh.node_ref(node->id()) = new_positions[node->id()];
117  }
118  } // end for
119  } // end scope
120 
121  // Now the nodes which are ghosts on this processor may have been moved on
122  // the processors which own them. So we need to synchronize with our neighbors
123  // and get the most up-to-date positions for the ghosts.
124  SyncNodalPositions sync_object(_mesh);
126  (_mesh.comm(), _mesh.nodes_begin(), _mesh.nodes_end(), sync_object);
127 
128  } // end for n_iterations
129 
130  // finally adjust the second order nodes (those located between vertices)
131  // these nodes will be located between their adjacent nodes
132  // do this element-wise
133  MeshBase::element_iterator el = _mesh.active_elements_begin();
134  const MeshBase::element_iterator end = _mesh.active_elements_end();
135 
136  for (; el != end; ++el)
137  {
138  // Constant handle for the element
139  const Elem * elem = *el;
140 
141  // get the second order nodes (son)
142  // their element indices start at n_vertices and go to n_nodes
143  const unsigned int son_begin = elem->n_vertices();
144  const unsigned int son_end = elem->n_nodes();
145 
146  // loop over all second order nodes (son)
147  for (unsigned int son=son_begin; son<son_end; son++)
148  {
149  // Don't smooth second-order nodes which are on the boundary
150  if (!on_boundary[elem->node_id(son)])
151  {
152  const unsigned int n_adjacent_vertices =
153  elem->n_second_order_adjacent_vertices(son);
154 
155  // calculate the new position which is the average of the
156  // position of the adjacent vertices
157  Point avg_position(0,0,0);
158  for (unsigned int v=0; v<n_adjacent_vertices; v++)
159  avg_position +=
160  _mesh.point( elem->node_id( elem->second_order_adjacent_vertex(son,v) ) );
161 
162  _mesh.node_ref(elem->node_id(son)) = avg_position / n_adjacent_vertices;
163  }
164  }
165  }
166 }
std::vector< std::vector< dof_id_type > > _graph
virtual const Point & point(const dof_id_type i) const =0
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:417
virtual dof_id_type max_node_id() const =0
const class libmesh_nullptr_t libmesh_nullptr
IterBase * end
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
void add(const TypeVector< T2 > &)
Definition: type_vector.h:596
virtual node_iterator nodes_begin()=0
virtual node_iterator local_nodes_end()=0
void sync_dofobject_data_by_id(const Communicator &comm, const Iterator &range_begin, const Iterator &range_end, SyncFunctor &sync)
virtual element_iterator active_elements_begin()=0
virtual element_iterator active_elements_end()=0
virtual node_iterator local_nodes_begin()=0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual node_iterator nodes_end()=0
const Parallel::Communicator & comm() const
void find_boundary_nodes(const MeshBase &mesh, std::vector< bool > &on_boundary)
Definition: mesh_tools.C:297
processor_id_type processor_id() const

Member Data Documentation

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

Data structure for holding the L-graph

Definition at line 104 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 99 of file mesh_smoother_laplace.h.

Referenced by init(), and smooth().


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