mesh_smoother_laplace.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 #include <algorithm> // for std::copy, std::sort
22 
23 
24 // Local includes
26 #include "libmesh/mesh_tools.h"
27 #include "libmesh/elem.h"
29 #include "libmesh/parallel.h"
30 #include "libmesh/parallel_ghost_sync.h" // sync_dofobject_data_by_id()
31 #include "libmesh/parallel_algebra.h" // StandardType<Point>
32 
33 namespace libMesh
34 {
35  // LaplaceMeshSmoother member functions
37  : MeshSmoother(mesh),
38  _initialized(false)
39  {
40  }
41 
42 
43 
44 
45 void LaplaceMeshSmoother::smooth(unsigned int n_iterations)
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.n_nodes())
58  {
59  libMesh::err << "MeshTools::find_boundary_nodes() returned incorrect length vector!" << std::endl;
60  libmesh_error();
61  }
62 
63  // We can only update the nodes after all new positions were
64  // determined. We store the new positions here
65  std::vector<Point> new_positions;
66 
67  for (unsigned int n=0; n<n_iterations; n++)
68  {
69  new_positions.resize(_mesh.n_nodes());
70 
71  {
74  for (; it != it_end; ++it)
75  {
76  Node* node = *it;
77 
78  if (node == NULL)
79  {
80  libMesh::err << "[" << _mesh.processor_id() << "]: Node iterator returned NULL pointer." << std::endl;
81  libmesh_error();
82  }
83 
84  // leave the boundary intact
85  // Only relocate the nodes which are vertices of an element
86  // All other entries of _graph (the secondary nodes) are empty
87  if (!on_boundary[node->id()] && (_graph[node->id()].size() > 0))
88  {
89  Point avg_position(0.,0.,0.);
90 
91  for (unsigned j=0; j<_graph[node->id()].size(); ++j)
92  {
93  // Will these nodal positions always be available
94  // or will they refer to remote nodes? To be
95  // careful, we grab a pointer and test it against
96  // NULL.
97  Node* connected_node = _mesh.node_ptr(_graph[node->id()][j]);
98 
99  if (connected_node == NULL)
100  {
101  libMesh::err << "Error! Libmesh returned NULL pointer for node " << _graph[connected_node->id()][j] << std::endl;
102  libmesh_error();
103  }
104 
105  avg_position.add( *connected_node );
106  } // end for(j)
107 
108  // Compute the average, store in the new_positions vector
109  new_positions[node->id()] = avg_position / static_cast<Real>(_graph[node->id()].size());
110  } // end if
111  } // end for
112  } // end scope
113 
114 
115  // now update the node positions (local node positions only)
116  {
119  for (; it != it_end; ++it)
120  {
121  Node* node = *it;
122 
123  if (!on_boundary[node->id()] && (_graph[node->id()].size() > 0))
124  {
125  // Should call Point::op=
126  // libMesh::out << "Setting node id " << node->id() << " to position " << new_positions[node->id()];
127  _mesh.node(node->id()) = new_positions[node->id()];
128  }
129  } // end for
130  } // end scope
131 
132  // Now the nodes which are ghosts on this processor may have been moved on
133  // the processors which own them. So we need to synchronize with our neighbors
134  // and get the most up-to-date positions for the ghosts.
135  SyncNodalPositions sync_object(_mesh);
137  (_mesh.comm(), _mesh.nodes_begin(), _mesh.nodes_end(), sync_object);
138 
139  } // end for n_iterations
140 
141  // finally adjust the second order nodes (those located between vertices)
142  // these nodes will be located between their adjacent nodes
143  // do this element-wise
146 
147  for (; el != end; ++el)
148  {
149  // Constant handle for the element
150  const Elem* elem = *el;
151 
152  // get the second order nodes (son)
153  // their element indices start at n_vertices and go to n_nodes
154  const unsigned int son_begin = elem->n_vertices();
155  const unsigned int son_end = elem->n_nodes();
156 
157  // loop over all second order nodes (son)
158  for (unsigned int son=son_begin; son<son_end; son++)
159  {
160  // Don't smooth second-order nodes which are on the boundary
161  if (!on_boundary[elem->node(son)])
162  {
163  const unsigned int n_adjacent_vertices =
165 
166  // calculate the new position which is the average of the
167  // position of the adjacent vertices
168  Point avg_position(0,0,0);
169  for (unsigned int v=0; v<n_adjacent_vertices; v++)
170  avg_position +=
171  _mesh.point( elem->node( elem->second_order_adjacent_vertex(son,v) ) );
172 
173  _mesh.node(elem->node(son)) = avg_position / n_adjacent_vertices;
174  }
175  }
176  }
177 }
178 
179 
180 
181 
182 
184 {
185  switch (_mesh.mesh_dimension())
186  {
187 
188  // TODO:[BSK] Fix this to work for refined meshes... I think
189  // the implementation was done quickly for Damien, who did not have
190  // refined grids. Fix it here and in the original Mesh member.
191 
192  case 2: // Stolen directly from build_L_graph in mesh_base.C
193  {
194  // Initialize space in the graph. It is n_nodes long. Each
195  // node may be connected to an arbitrary number of other nodes
196  // via edges.
197  _graph.resize(_mesh.n_nodes());
198 
201 
202  for (; el != end; ++el)
203  {
204  // Constant handle for the element
205  const Elem* elem = *el;
206 
207  for (unsigned int s=0; s<elem->n_neighbors(); s++)
208  {
209  // Only operate on sides which are on the
210  // boundary or for which the current element's
211  // id is greater than its neighbor's.
212  // Sides get only built once.
213  if ((elem->neighbor(s) == NULL) ||
214  (elem->id() > elem->neighbor(s)->id()))
215  {
216  AutoPtr<Elem> side(elem->build_side(s));
217  _graph[side->node(0)].push_back(side->node(1));
218  _graph[side->node(1)].push_back(side->node(0));
219  }
220  }
221  }
222  _initialized = true;
223  break;
224  } // case 2
225 
226  case 3: // Stolen blatantly from build_L_graph in mesh_base.C
227  {
228  // Initialize space in the graph.
229  _graph.resize(_mesh.n_nodes());
230 
233 
234  for (; el != end; ++el)
235  {
236  // Shortcut notation for simplicity
237  const Elem* elem = *el;
238 
239  for (unsigned int f=0; f<elem->n_neighbors(); f++) // Loop over faces
240  if ((elem->neighbor(f) == NULL) ||
241  (elem->id() > elem->neighbor(f)->id()))
242  {
243  // We need a full (i.e. non-proxy) element for the face, since we will
244  // be looking at its sides as well!
245  AutoPtr<Elem> face = elem->build_side(f, /*proxy=*/false);
246 
247  for (unsigned int s=0; s<face->n_neighbors(); s++) // Loop over face's edges
248  {
249  // Here we can use a proxy
250  AutoPtr<Elem> side = face->build_side(s);
251 
252  // At this point, we just insert the node numbers
253  // again. At the end we'll call sort and unique
254  // to make sure there are no duplicates
255  _graph[side->node(0)].push_back(side->node(1));
256  _graph[side->node(1)].push_back(side->node(0));
257  }
258  }
259  }
260 
261  _initialized = true;
262  break;
263  } // case 3
264 
265  default:
266  {
267  libMesh::err << "At this time it is not possible "
268  << "to smooth a dimension "
269  << _mesh.mesh_dimension()
270  << "mesh. Aborting..."
271  << std::endl;
272  libmesh_error();
273  }
274  }
275 
276  // Done building graph from local elements. Let's now allgather the
277  // graph so that it is available on all processors for the actual
278  // smoothing operation?
279  this->allgather_graph();
280 
281  // In 3D, it's possible for > 2 processor partitions to meet
282  // at a single edge, while in 2D only 2 processor partitions
283  // share an edge. Therefore the allgather'd graph in 3D may
284  // now have duplicate entries and we need to remove them so
285  // they don't foul up the averaging algorithm employed by the
286  // Laplace smoother.
287  for (unsigned i=0; i<_graph.size(); ++i)
288  {
289  // The std::unique algorithm removes duplicate *consecutive* elements from a range,
290  // so it only makes sense to call it on a sorted range...
291  std::sort(_graph[i].begin(), _graph[i].end());
292  _graph[i].erase(std::unique(_graph[i].begin(), _graph[i].end()), _graph[i].end());
293  }
294 
295 } // init()
296 
297 
298 
299 
300 void LaplaceMeshSmoother::print_graph(std::ostream& out_stream) const
301 {
302  for (unsigned int i=0; i<_graph.size(); ++i)
303  {
304  out_stream << i << ": ";
305  std::copy(_graph[i].begin(),
306  _graph[i].end(),
307  std::ostream_iterator<unsigned>(out_stream, " "));
308  out_stream << std::endl;
309  }
310 }
311 
312 
313 
315  {
316  // The graph data structure is not well-suited for parallel communication,
317  // so copy the graph into a single vector defined by:
318  // NA A_0 A_1 ... A_{NA} | NB B_0 B_1 ... B_{NB} | NC C_0 C_1 ... C_{NC}
319  // where:
320  // * NA is the number of graph connections for node A
321  // * A_0, A_1, etc. are the IDs connected to node A
322  std::vector<std::size_t> flat_graph;
323 
324  // Reserve at least enough space for each node to have zero entries
325  flat_graph.reserve(_graph.size());
326 
327  for (std::size_t i=0; i<_graph.size(); ++i)
328  {
329  // First push back the number of entries for this node
330  flat_graph.push_back (_graph[i].size());
331 
332  // Then push back all the IDs
333  for (std::size_t j=0; j<_graph[i].size(); ++j)
334  flat_graph.push_back(_graph[i][j]);
335  }
336 
337  // // A copy of the flat graph (for printing only, delete me later)
338  // std::vector<unsigned> copy_of_flat_graph(flat_graph);
339 
340  // Use the allgather routine to combine all the flat graphs on all processors
341  _mesh.comm().allgather(flat_graph);
342 
343  // Now reconstruct _graph from the allgathered flat_graph.
344 
345  // // (Delete me later, the copy is just for printing purposes.)
346  // std::vector<std::vector<unsigned > > copy_of_graph(_graph);
347 
348  // Make sure the old graph is cleared out
349  _graph.clear();
350  _graph.resize(_mesh.n_nodes());
351 
352  // Our current position in the allgather'd flat_graph
353  std::size_t cursor=0;
354 
355  // There are n_nodes * n_processors vectors to read in total
356  for (processor_id_type p=0; p<_mesh.n_processors(); ++p)
357  for (dof_id_type node_ctr=0; node_ctr<_mesh.n_nodes(); ++node_ctr)
358  {
359  // Read the number of entries for this node, move cursor
360  std::size_t n_entries = flat_graph[cursor++];
361 
362  // Reserve space for that many more entries, then push back
363  _graph[node_ctr].reserve(_graph[node_ctr].size() + n_entries);
364 
365  // Read all graph connections for this node, move the cursor each time
366  // Note: there might be zero entries but that's fine
367  for (std::size_t i=0; i<n_entries; ++i)
368  _graph[node_ctr].push_back(flat_graph[cursor++]);
369  }
370 
371 
372 // // Print local graph to uniquely named file (debugging)
373 // {
374 // // Generate unique filename for this processor
375 // std::ostringstream oss;
376 // oss << "graph_filename_" << _mesh.processor_id() << ".txt";
377 // std::ofstream graph_stream(oss.str().c_str());
378 //
379 // // Print the local non-flat graph
380 // std::swap(_graph, copy_of_graph);
381 // print_graph(graph_stream);
382 //
383 // // Print the (local) flat graph for verification
384 // for (std::size_t i=0; i<copy_of_flat_graph.size(); ++i)
385 // graph_stream << copy_of_flat_graph[i] << " ";
386 // graph_stream << "\n";
387 //
388 // // Print the allgather'd grap for verification
389 // for (std::size_t i=0; i<flat_graph.size(); ++i)
390 // graph_stream << flat_graph[i] << " ";
391 // graph_stream << "\n";
392 //
393 // // Print the global non-flat graph
394 // std::swap(_graph, copy_of_graph);
395 // print_graph(graph_stream);
396 // }
397  } // allgather_graph()
398 
399 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo