libMesh
mesh_smoother_laplace.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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
25 #include "libmesh/mesh_smoother_laplace.h"
26 #include "libmesh/mesh_tools.h"
27 #include "libmesh/elem.h"
28 #include "libmesh/unstructured_mesh.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.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  for (auto & node : _mesh.local_node_ptr_range())
69  {
70  if (node == libmesh_nullptr)
71  libmesh_error_msg("[" << _mesh.processor_id() << "]: Node iterator returned NULL pointer.");
72 
73  // leave the boundary intact
74  // Only relocate the nodes which are vertices of an element
75  // All other entries of _graph (the secondary nodes) are empty
76  if (!on_boundary[node->id()] && (_graph[node->id()].size() > 0))
77  {
78  Point avg_position(0.,0.,0.);
79 
80  for (std::size_t j=0; j<_graph[node->id()].size(); ++j)
81  {
82  // Will these nodal positions always be available
83  // or will they refer to remote nodes? This will
84  // fail an assertion in the latter case, which
85  // shouldn't occur if DistributedMesh is working
86  // correctly.
87  const Point & connected_node = _mesh.point(_graph[node->id()][j]);
88 
89  avg_position.add( connected_node );
90  } // end for (j)
91 
92  // Compute the average, store in the new_positions vector
93  new_positions[node->id()] = avg_position / static_cast<Real>(_graph[node->id()].size());
94  } // end if
95  } // end for
96 
97 
98  // now update the node positions (local node positions only)
99  for (auto & node : _mesh.local_node_ptr_range())
100  if (!on_boundary[node->id()] && (_graph[node->id()].size() > 0))
101  {
102  // Should call Point::op=
103  // libMesh::out << "Setting node id " << node->id() << " to position " << new_positions[node->id()];
104  _mesh.node_ref(node->id()) = new_positions[node->id()];
105  }
106 
107  // Now the nodes which are ghosts on this processor may have been moved on
108  // the processors which own them. So we need to synchronize with our neighbors
109  // and get the most up-to-date positions for the ghosts.
110  SyncNodalPositions sync_object(_mesh);
112  (_mesh.comm(), _mesh.nodes_begin(), _mesh.nodes_end(), sync_object);
113 
114  } // end for n_iterations
115 
116  // finally adjust the second order nodes (those located between vertices)
117  // these nodes will be located between their adjacent nodes
118  // do this element-wise
119  for (auto & elem : _mesh.active_element_ptr_range())
120  {
121  // get the second order nodes (son)
122  // their element indices start at n_vertices and go to n_nodes
123  const unsigned int son_begin = elem->n_vertices();
124  const unsigned int son_end = elem->n_nodes();
125 
126  // loop over all second order nodes (son)
127  for (unsigned int son=son_begin; son<son_end; son++)
128  {
129  // Don't smooth second-order nodes which are on the boundary
130  if (!on_boundary[elem->node_id(son)])
131  {
132  const unsigned int n_adjacent_vertices =
133  elem->n_second_order_adjacent_vertices(son);
134 
135  // calculate the new position which is the average of the
136  // position of the adjacent vertices
137  Point avg_position(0,0,0);
138  for (unsigned int v=0; v<n_adjacent_vertices; v++)
139  avg_position +=
140  _mesh.point( elem->node_id( elem->second_order_adjacent_vertex(son,v) ) );
141 
142  _mesh.node_ref(elem->node_id(son)) = avg_position / n_adjacent_vertices;
143  }
144  }
145  }
146 }
147 
148 
149 
150 
151 
153 {
154  switch (_mesh.mesh_dimension())
155  {
156 
157  // TODO:[BSK] Fix this to work for refined meshes... I think
158  // the implementation was done quickly for Damien, who did not have
159  // refined grids. Fix it here and in the original Mesh member.
160 
161  case 2: // Stolen directly from build_L_graph in mesh_base.C
162  {
163  // Initialize space in the graph. It is indexed by node id.
164  // Each node may be connected to an arbitrary number of other
165  // nodes via edges.
166  _graph.resize(_mesh.max_node_id());
167 
168  for (auto & elem : _mesh.active_local_element_ptr_range())
169  for (auto s : elem->side_index_range())
170  {
171  // Only operate on sides which are on the
172  // boundary or for which the current element's
173  // id is greater than its neighbor's.
174  // Sides get only built once.
175  if ((elem->neighbor_ptr(s) == libmesh_nullptr) ||
176  (elem->id() > elem->neighbor_ptr(s)->id()))
177  {
178  UniquePtr<const Elem> side(elem->build_side_ptr(s));
179  _graph[side->node_id(0)].push_back(side->node_id(1));
180  _graph[side->node_id(1)].push_back(side->node_id(0));
181  }
182  }
183  _initialized = true;
184  break;
185  } // case 2
186 
187  case 3: // Stolen blatantly from build_L_graph in mesh_base.C
188  {
189  // Initialize space in the graph.
190  _graph.resize(_mesh.max_node_id());
191 
192  for (auto & elem : _mesh.active_local_element_ptr_range())
193  for (auto f : elem->side_index_range()) // Loop over faces
194  if ((elem->neighbor_ptr(f) == libmesh_nullptr) ||
195  (elem->id() > elem->neighbor_ptr(f)->id()))
196  {
197  // We need a full (i.e. non-proxy) element for the face, since we will
198  // be looking at its sides as well!
199  UniquePtr<const Elem> face = elem->build_side_ptr(f, /*proxy=*/false);
200 
201  for (auto s : face->side_index_range()) // Loop over face's edges
202  {
203  // Here we can use a proxy
204  UniquePtr<const Elem> side = face->build_side_ptr(s);
205 
206  // At this point, we just insert the node numbers
207  // again. At the end we'll call sort and unique
208  // to make sure there are no duplicates
209  _graph[side->node_id(0)].push_back(side->node_id(1));
210  _graph[side->node_id(1)].push_back(side->node_id(0));
211  }
212  }
213 
214  _initialized = true;
215  break;
216  } // case 3
217 
218  default:
219  libmesh_error_msg("At this time it is not possible to smooth a dimension " << _mesh.mesh_dimension() << "mesh. Aborting...");
220  }
221 
222  // Done building graph from local elements. Let's now allgather the
223  // graph so that it is available on all processors for the actual
224  // smoothing operation?
225  this->allgather_graph();
226 
227  // In 3D, it's possible for > 2 processor partitions to meet
228  // at a single edge, while in 2D only 2 processor partitions
229  // share an edge. Therefore the allgather'd graph in 3D may
230  // now have duplicate entries and we need to remove them so
231  // they don't foul up the averaging algorithm employed by the
232  // Laplace smoother.
233  for (std::size_t i=0; i<_graph.size(); ++i)
234  {
235  // The std::unique algorithm removes duplicate *consecutive* elements from a range,
236  // so it only makes sense to call it on a sorted range...
237  std::sort(_graph[i].begin(), _graph[i].end());
238  _graph[i].erase(std::unique(_graph[i].begin(), _graph[i].end()), _graph[i].end());
239  }
240 
241 } // init()
242 
243 
244 
245 
246 void LaplaceMeshSmoother::print_graph(std::ostream & out_stream) const
247 {
248  for (std::size_t i=0; i<_graph.size(); ++i)
249  {
250  out_stream << i << ": ";
251  std::copy(_graph[i].begin(),
252  _graph[i].end(),
253  std::ostream_iterator<unsigned>(out_stream, " "));
254  out_stream << std::endl;
255  }
256 }
257 
258 
259 
261 {
262  // The graph data structure is not well-suited for parallel communication,
263  // so copy the graph into a single vector defined by:
264  // NA A_0 A_1 ... A_{NA} | NB B_0 B_1 ... B_{NB} | NC C_0 C_1 ... C_{NC}
265  // where:
266  // * NA is the number of graph connections for node A
267  // * A_0, A_1, etc. are the IDs connected to node A
268  std::vector<dof_id_type> flat_graph;
269 
270  // Reserve at least enough space for each node to have zero entries
271  flat_graph.reserve(_graph.size());
272 
273  for (std::size_t i=0; i<_graph.size(); ++i)
274  {
275  // First push back the number of entries for this node
276  flat_graph.push_back (cast_int<dof_id_type>(_graph[i].size()));
277 
278  // Then push back all the IDs
279  for (std::size_t j=0; j<_graph[i].size(); ++j)
280  flat_graph.push_back(_graph[i][j]);
281  }
282 
283  // // A copy of the flat graph (for printing only, delete me later)
284  // std::vector<unsigned> copy_of_flat_graph(flat_graph);
285 
286  // Use the allgather routine to combine all the flat graphs on all processors
287  _mesh.comm().allgather(flat_graph);
288 
289  // Now reconstruct _graph from the allgathered flat_graph.
290 
291  // // (Delete me later, the copy is just for printing purposes.)
292  // std::vector<std::vector<unsigned >> copy_of_graph(_graph);
293 
294  // Make sure the old graph is cleared out
295  _graph.clear();
296  _graph.resize(_mesh.max_node_id());
297 
298  // Our current position in the allgather'd flat_graph
299  std::size_t cursor=0;
300 
301  // There are max_node_id * n_processors entries to read in total
302  for (processor_id_type p=0; p<_mesh.n_processors(); ++p)
303  for (dof_id_type node_ctr=0; node_ctr<_mesh.max_node_id(); ++node_ctr)
304  {
305  // Read the number of entries for this node, move cursor
306  std::size_t n_entries = flat_graph[cursor++];
307 
308  // Reserve space for that many more entries, then push back
309  _graph[node_ctr].reserve(_graph[node_ctr].size() + n_entries);
310 
311  // Read all graph connections for this node, move the cursor each time
312  // Note: there might be zero entries but that's fine
313  for (std::size_t i=0; i<n_entries; ++i)
314  _graph[node_ctr].push_back(flat_graph[cursor++]);
315  }
316 
317 
318  // // Print local graph to uniquely named file (debugging)
319  // {
320  // // Generate unique filename for this processor
321  // std::ostringstream oss;
322  // oss << "graph_filename_" << _mesh.processor_id() << ".txt";
323  // std::ofstream graph_stream(oss.str().c_str());
324  //
325  // // Print the local non-flat graph
326  // std::swap(_graph, copy_of_graph);
327  // print_graph(graph_stream);
328  //
329  // // Print the (local) flat graph for verification
330  // for (std::size_t i=0; i<copy_of_flat_graph.size(); ++i)
331  //graph_stream << copy_of_flat_graph[i] << " ";
332  // graph_stream << "\n";
333  //
334  // // Print the allgather'd grap for verification
335  // for (std::size_t i=0; i<flat_graph.size(); ++i)
336  //graph_stream << flat_graph[i] << " ";
337  // graph_stream << "\n";
338  //
339  // // Print the global non-flat graph
340  // std::swap(_graph, copy_of_graph);
341  // print_graph(graph_stream);
342  // }
343 } // allgather_graph()
344 
345 } // namespace libMesh
std::vector< std::vector< dof_id_type > > _graph
Data structure for holding the L-graph.
bool _initialized
True if the L-graph has been created, false otherwise.
virtual const Point & point(const dof_id_type i) const =0
virtual SimpleRange< node_iterator > local_node_ptr_range()=0
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:420
unsigned short int side
Definition: xdr_io.C:49
processor_id_type n_processors() const
virtual dof_id_type max_node_id() const =0
MeshBase & mesh
uint8_t processor_id_type
Definition: id_types.h:99
const class libmesh_nullptr_t libmesh_nullptr
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
void init()
Initialization for the Laplace smoothing routine is basically identical to building an "L-graph" whic...
UnstructuredMesh & _mesh
Definition: mesh_smoother.h:61
The libMesh namespace provides an interface to certain functionality in the library.
void add(const TypeVector< T2 > &)
Add to this vector without creating a temporary.
Definition: type_vector.h:600
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual node_iterator nodes_begin()=0
Iterate over all the nodes in the Mesh.
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
The UnstructuredMesh class is derived from the MeshBase class.
void sync_dofobject_data_by_id(const Communicator &comm, const Iterator &range_begin, const Iterator &range_end, SyncFunctor &sync)
Request data about a range of ghost dofobjects uniquely identified by their id.
LaplaceMeshSmoother(UnstructuredMesh &mesh)
Constructor.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual node_iterator nodes_end()=0
void allgather_graph()
This function allgather&#39;s the (local) graph after it is computed on each processor by the init() func...
const Parallel::Communicator & comm() const
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
void print_graph(std::ostream &out=libMesh::out) const
Mainly for debugging, this function will print out the connectivity graph which has been created...
This class provides the necessary interface for mesh smoothing.
Definition: mesh_smoother.h:38
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
void find_boundary_nodes(const MeshBase &mesh, std::vector< bool > &on_boundary)
Calling this function on a 2D mesh will convert all the elements to triangles.
Definition: mesh_tools.C:290
virtual void smooth() libmesh_override
Redefinition of the smooth function from the base class.
uint8_t dof_id_type
Definition: id_types.h:64
processor_id_type processor_id() const
void allgather(const T &send, std::vector< T > &recv) const
Take a vector of length this->size(), and fill in recv[processor_id] = the value of send on that proc...