 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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 }
1.9.1.