Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
LloydRelaxation.cpp
Go to the documentation of this file.
1 /** @example LloydRelaxation.cpp \n
2  * \brief Perform Lloyd relaxation on a mesh and its dual \n
3  * <b>To run</b>: mpiexec -np <np> LloydRelaxation [filename]\n
4  *
5  * Briefly, Lloyd relaxation is a technique to smooth out a mesh. The centroid of each cell is
6  * computed from its vertex positions, then vertices are placed at the average of their connected
7  * cells' centroids.
8  *
9  * In the parallel algorithm, an extra ghost layer of cells is exchanged. This allows us to compute
10  * the centroids for boundary cells on each processor where they appear; this eliminates the need
11  * for one round of data exchange (for those centroids) between processors. New vertex positions
12  * must be sent from owning processors to processors sharing those vertices. Convergence is
13  * measured as the maximum distance moved by any vertex.
14  *
15  * In this implementation, a fixed number of iterations is performed. The final mesh is output to
16  * 'lloydfinal.h5m' in the current directory (H5M format must be used since the file is written in
17  * parallel).
18  */
19 
20 #include "moab/Core.hpp"
21 #ifdef MOAB_HAVE_MPI
22 #include "moab/ParallelComm.hpp"
23 #include "MBParallelConventions.h"
24 #endif
25 #include "moab/Skinner.hpp"
26 #include "moab/CN.hpp"
27 #include "moab/CartVect.hpp"
28 #include <iostream>
29 #include <sstream>
30 
31 using namespace moab;
32 using namespace std;
33 
34 string test_file_name = string( MESH_DIR ) + string( "/surfrandomtris-4part.h5m" );
35 
36 ErrorCode perform_lloyd_relaxation( Interface* mb, Range& verts, Range& cells, Tag fixed, int num_its, int report_its );
37 
38 int main( int argc, char** argv )
39 {
40  int num_its = 10;
41  int report_its = 1;
42 
43 #ifdef MOAB_HAVE_MPI
44  MPI_Init( &argc, &argv );
45 #endif
46 
47  // Need option handling here for input filename
48  if( argc > 1 )
49  {
50  // User has input a mesh file
51  test_file_name = argv[1];
52  }
53 
54  // Get MOAB instance
55  Interface* mb = new( std::nothrow ) Core;
56  if( NULL == mb ) return 1;
57 
58 #ifdef MOAB_HAVE_MPI
59  // Get the ParallelComm instance
60  ParallelComm* pcomm = new ParallelComm( mb, MPI_COMM_WORLD );
61  int nprocs = pcomm->size();
62 #else
63  int nprocs = 1;
64 #endif
65  string options;
66  if( nprocs > 1 ) // If reading in parallel, need to tell it how
67  options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"
68  "PARALLEL_GHOSTS=2.0.1;DEBUG_IO=0;DEBUG_PIO=0";
69 
70  // Load the test file with specified options
71  ErrorCode rval = mb->load_file( test_file_name.c_str(), 0, options.c_str() );MB_CHK_ERR( rval );
72 
73  // Make tag to specify fixed vertices, since it's input to the algorithm; use a default value of
74  // non-fixed so we only need to set the fixed tag for skin vertices
75  Tag fixed;
76  int def_val = 0;
77  rval = mb->tag_get_handle( "fixed", 1, MB_TYPE_INTEGER, fixed, MB_TAG_CREAT | MB_TAG_DENSE, &def_val );MB_CHK_ERR( rval );
78 
79  // Get all vertices and faces
80  Range verts, faces, skin_verts;
81  rval = mb->get_entities_by_type( 0, MBVERTEX, verts );MB_CHK_ERR( rval );
82  rval = mb->get_entities_by_dimension( 0, 2, faces );MB_CHK_ERR( rval );
83 
84  // Get the skin vertices of those faces and mark them as fixed; we don't want to fix the
85  // vertices on a part boundary, but since we exchanged a layer of ghost faces, those vertices
86  // aren't on the skin locally ok to mark non-owned skin vertices too, I won't move those anyway
87  // use MOAB's skinner class to find the skin
88  Skinner skinner( mb );
89  rval = skinner.find_skin( 0, faces, true, skin_verts );MB_CHK_ERR( rval ); // 'true' param indicates we want vertices back, not faces
90 
91  vector< int > fix_tag( skin_verts.size(), 1 ); // Initialized to 1 to indicate fixed
92  rval = mb->tag_set_data( fixed, skin_verts, &fix_tag[0] );MB_CHK_ERR( rval );
93 
94  // Now perform the Lloyd relaxation
95  rval = perform_lloyd_relaxation( mb, verts, faces, fixed, num_its, report_its );MB_CHK_ERR( rval );
96 
97  // Delete fixed tag, since we created it here
98  rval = mb->tag_delete( fixed );MB_CHK_ERR( rval );
99 
100  // Output file, using parallel write
101 
102 #ifdef MOAB_HAVE_MPI
103  options = "PARALLEL=WRITE_PART";
104 #endif
105 
106  rval = mb->write_file( "lloydfinal.h5m", NULL, options.c_str() );MB_CHK_ERR( rval );
107 
108  // Delete MOAB instance
109  delete mb;
110 
111 #ifdef MOAB_HAVE_MPI
112  MPI_Finalize();
113 #endif
114 
115  return 0;
116 }
117 
118 ErrorCode perform_lloyd_relaxation( Interface* mb, Range& verts, Range& faces, Tag fixed, int num_its, int report_its )
119 {
120  ErrorCode rval;
121 
122 #ifdef MOAB_HAVE_MPI
123  ParallelComm* pcomm = ParallelComm::get_pcomm( mb, 0 );
124  int nprocs = pcomm->size();
125 #else
126  int nprocs = 1;
127 #endif
128 
129  // Perform Lloyd relaxation:
130  // 1. Setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
131 
132  // Get all verts coords into tag; don't need to worry about filtering out fixed verts,
133  // we'll just be setting to their fixed coords
134  vector< double > vcentroids( 3 * verts.size() );
135  rval = mb->get_coords( verts, &vcentroids[0] );MB_CHK_ERR( rval );
136 
137  Tag centroid;
138  rval = mb->tag_get_handle( "centroid", 3, MB_TYPE_DOUBLE, centroid, MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_ERR( rval );
139  rval = mb->tag_set_data( centroid, verts, &vcentroids[0] );MB_CHK_ERR( rval );
140 
141  // Filter verts down to owned ones and get fixed tag for them
142  Range owned_verts, shared_owned_verts;
143  if( nprocs > 1 )
144  {
145 #ifdef MOAB_HAVE_MPI
146  rval = pcomm->filter_pstatus( verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts );MB_CHK_ERR( rval );
147 #endif
148  }
149  else
150  owned_verts = verts;
151  vector< int > fix_tag( owned_verts.size() );
152  rval = mb->tag_get_data( fixed, owned_verts, &fix_tag[0] );MB_CHK_ERR( rval );
153 
154  // Now fill vcentroids array with positions of just owned vertices, since those are the ones
155  // we're actually computing
156  vcentroids.resize( 3 * owned_verts.size() );
157  rval = mb->tag_get_data( centroid, owned_verts, &vcentroids[0] );MB_CHK_ERR( rval );
158 
159 #ifdef MOAB_HAVE_MPI
160  // Get shared owned verts, for exchanging tags
161  rval = pcomm->get_shared_entities( -1, shared_owned_verts, 0, false, true );MB_CHK_ERR( rval );
162  // Workaround: if no shared owned verts, put a non-shared one in the list, to prevent exchanging
163  // tags for all shared entities
164  if( shared_owned_verts.empty() ) shared_owned_verts.insert( *verts.begin() );
165 #endif
166 
167  // Some declarations for later iterations
168  vector< double > fcentroids( 3 * faces.size() ); // fcentroids for face centroids
169  vector< double > ctag( 3 * CN::MAX_NODES_PER_ELEMENT ); // Temporary coordinate storage for verts bounding a face
170  const EntityHandle* conn; // const ptr & size to face connectivity
171  int nconn;
172  Range::iterator fit, vit; // For iterating over faces, verts
173  int f, v; // For indexing into centroid vectors
174  vector< EntityHandle > adj_faces; // Used in vertex iteration
175 
176  // 2. For num_its iterations:
177  for( int nit = 0; nit < num_its; nit++ )
178  {
179  double mxdelta = 0.0;
180 
181  // 2a. For each face: centroid = sum(vertex centroids)/num_verts_in_cell
182  for( fit = faces.begin(), f = 0; fit != faces.end(); ++fit, f++ )
183  {
184  // Get verts for this face
185  rval = mb->get_connectivity( *fit, conn, nconn );MB_CHK_ERR( rval );
186  // Get centroid tags for those verts
187  rval = mb->tag_get_data( centroid, conn, nconn, &ctag[0] );MB_CHK_ERR( rval );
188  fcentroids[3 * f + 0] = fcentroids[3 * f + 1] = fcentroids[3 * f + 2] = 0.0;
189  for( v = 0; v < nconn; v++ )
190  {
191  fcentroids[3 * f + 0] += ctag[3 * v + 0];
192  fcentroids[3 * f + 1] += ctag[3 * v + 1];
193  fcentroids[3 * f + 2] += ctag[3 * v + 2];
194  }
195  for( v = 0; v < 3; v++ )
196  fcentroids[3 * f + v] /= nconn;
197  }
198  rval = mb->tag_set_data( centroid, faces, &fcentroids[0] );MB_CHK_ERR( rval );
199 
200  // 2b. For each owned vertex:
201  for( vit = owned_verts.begin(), v = 0; vit != owned_verts.end(); ++vit, v++ )
202  {
203  // if !fixed
204  if( fix_tag[v] ) continue;
205  // vertex centroid = sum(cell centroids)/ncells
206  adj_faces.clear();
207  rval = mb->get_adjacencies( &( *vit ), 1, 2, false, adj_faces );MB_CHK_ERR( rval );
208  rval = mb->tag_get_data( centroid, &adj_faces[0], adj_faces.size(), &fcentroids[0] );MB_CHK_ERR( rval );
209  double vnew[] = { 0.0, 0.0, 0.0 };
210  for( f = 0; f < (int)adj_faces.size(); f++ )
211  {
212  vnew[0] += fcentroids[3 * f + 0];
213  vnew[1] += fcentroids[3 * f + 1];
214  vnew[2] += fcentroids[3 * f + 2];
215  }
216  for( f = 0; f < 3; f++ )
217  vnew[f] /= adj_faces.size();
218  double delta = ( CartVect( vnew ) - CartVect( &vcentroids[3 * v] ) ).length();
219  mxdelta = std::max( delta, mxdelta );
220  for( f = 0; f < 3; f++ )
221  vcentroids[3 * v + f] = vnew[f];
222  }
223 
224  // Set the centroid tag; having them only in vcentroids array isn't enough, as vertex
225  // centroids are accessed randomly in loop over faces
226  rval = mb->tag_set_data( centroid, owned_verts, &vcentroids[0] );MB_CHK_ERR( rval );
227 
228  // 2c. Exchange tags on owned verts
229  if( nprocs > 1 )
230  {
231 #ifdef MOAB_HAVE_MPI
232  rval = pcomm->exchange_tags( centroid, shared_owned_verts );MB_CHK_ERR( rval );
233 #endif
234  }
235 
236  if( !( nit % report_its ) )
237  {
238  // Global reduce for maximum delta, then report it
239  double global_max = mxdelta;
240  int myrank = 0;
241 #ifdef MOAB_HAVE_MPI
242  if( nprocs > 1 ) MPI_Reduce( &mxdelta, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm() );
243  myrank = pcomm->rank();
244 #endif
245  if( 1 == nprocs || 0 == myrank ) cout << "Max delta = " << global_max << endl;
246  }
247  }
248 
249  // Write the tag back onto vertex coordinates
250  rval = mb->set_coords( owned_verts, &vcentroids[0] );MB_CHK_ERR( rval );
251 
252  // Delete the centroid tag, since we don't need it anymore
253  rval = mb->tag_delete( centroid );MB_CHK_ERR( rval );
254 
255  return MB_SUCCESS;
256 }