Branch data Line data Source code
1 : : #include "moab/ParallelMergeMesh.hpp"
2 : : #include "moab/Core.hpp"
3 : : #include "moab/CartVect.hpp"
4 : : #include "moab/BoundBox.hpp"
5 : : #include "moab/Skinner.hpp"
6 : : #include "moab/MergeMesh.hpp"
7 : : #include "moab/CN.hpp"
8 : : #include "float.h"
9 : : #include <algorithm>
10 : :
11 : : #ifdef MOAB_HAVE_MPI
12 : : #include "moab_mpi.h"
13 : : #endif
14 : :
15 : : namespace moab
16 : : {
17 : :
18 : : // Constructor
19 : : /*Get Merge Data and tolerance*/
20 [ # # ][ # # ]: 0 : ParallelMergeMesh::ParallelMergeMesh( ParallelComm* pc, const double epsilon ) : myPcomm( pc ), myEps( epsilon )
[ # # ]
21 : : {
22 [ # # ]: 0 : myMB = pc->get_moab();
23 [ # # ]: 0 : mySkinEnts.resize( 4 );
24 : 0 : }
25 : :
26 : : // Have a wrapper function on the actual merge to avoid memory leaks
27 : : // Merges elements within a proximity of epsilon
28 : 0 : ErrorCode ParallelMergeMesh::merge( EntityHandle levelset, bool skip_local_merge, int dim )
29 : : {
30 [ # # ][ # # ]: 0 : ErrorCode rval = PerformMerge( levelset, skip_local_merge, dim );MB_CHK_ERR( rval );
31 : 0 : CleanUp();
32 : 0 : return rval;
33 : : }
34 : :
35 : : // Perform the merge
36 : 0 : ErrorCode ParallelMergeMesh::PerformMerge( EntityHandle levelset, bool skip_local_merge, int dim )
37 : : {
38 : : // Get the mesh dimension
39 : : ErrorCode rval;
40 [ # # ]: 0 : if( dim < 0 )
41 : : {
42 [ # # ][ # # ]: 0 : rval = myMB->get_dimension( dim );MB_CHK_ERR( rval );
[ # # ][ # # ]
43 : : }
44 : :
45 : : // Get the local skin elements
46 [ # # ]: 0 : rval = PopulateMySkinEnts( levelset, dim, skip_local_merge );
47 : : // If there is only 1 proc, we can return now
48 [ # # ][ # # ]: 0 : if( rval != MB_SUCCESS || myPcomm->size() == 1 ) { return rval; }
[ # # ][ # # ]
49 : :
50 : : // Determine the global bounding box
51 : : double gbox[6];
52 [ # # ][ # # ]: 0 : rval = GetGlobalBox( gbox );MB_CHK_ERR( rval );
[ # # ][ # # ]
53 : :
54 : : /* Assemble The Destination Tuples */
55 : : // Get a list of tuples which contain (toProc, handle, x,y,z)
56 [ # # ][ # # ]: 0 : myTup.initialize( 1, 0, 1, 3, mySkinEnts[0].size() );
[ # # ]
57 [ # # ][ # # ]: 0 : rval = PopulateMyTup( gbox );MB_CHK_ERR( rval );
[ # # ][ # # ]
58 : :
59 : : /* Gather-Scatter Tuple
60 : : -tup comes out as (remoteProc,handle,x,y,z) */
61 [ # # ][ # # ]: 0 : myCD.initialize( myPcomm->comm() );
62 : :
63 : : // 1 represents dynamic tuple, 0 represents index of the processor to send to
64 [ # # ]: 0 : myCD.gs_transfer( 1, myTup, 0 );
65 : :
66 : : /* Sort By X,Y,Z
67 : : -Utilizes a custom quick sort incorporating epsilon*/
68 [ # # ]: 0 : SortTuplesByReal( myTup, myEps );
69 : :
70 : : // Initialize another tuple list for matches
71 [ # # ][ # # ]: 0 : myMatches.initialize( 2, 0, 2, 0, mySkinEnts[0].size() );
[ # # ]
72 : :
73 : : // ID the matching tuples
74 [ # # ][ # # ]: 0 : rval = PopulateMyMatches();MB_CHK_ERR( rval );
[ # # ][ # # ]
75 : :
76 : : // We can free up the tuple myTup now
77 [ # # ]: 0 : myTup.reset();
78 : :
79 : : /*Gather-Scatter Again*/
80 : : // 1 represents dynamic list, 0 represents proc index to send tuple to
81 [ # # ]: 0 : myCD.gs_transfer( 1, myMatches, 0 );
82 : : // We can free up the crystal router now
83 [ # # ]: 0 : myCD.reset();
84 : :
85 : : // Sort the matches tuple list
86 [ # # ]: 0 : SortMyMatches();
87 : :
88 : : // Tag the shared elements
89 [ # # ][ # # ]: 0 : rval = TagSharedElements( dim );MB_CHK_ERR( rval );
[ # # ][ # # ]
90 : :
91 : : // Free up the matches tuples
92 [ # # ]: 0 : myMatches.reset();
93 : 0 : return rval;
94 : : }
95 : :
96 : : // Sets mySkinEnts with all of the skin entities on the processor
97 : 0 : ErrorCode ParallelMergeMesh::PopulateMySkinEnts( const EntityHandle meshset, int dim, bool skip_local_merge )
98 : : {
99 : : /*Merge Mesh Locally*/
100 : : // Get all dim dimensional entities
101 [ # # ]: 0 : Range ents;
102 [ # # ][ # # ]: 0 : ErrorCode rval = myMB->get_entities_by_dimension( meshset, dim, ents );MB_CHK_ERR( rval );
[ # # ][ # # ]
103 : :
104 [ # # ][ # # ]: 0 : if( ents.empty() && dim == 3 )
[ # # ][ # # ]
105 : : {
106 : 0 : dim--;
107 [ # # ][ # # ]: 0 : rval = myMB->get_entities_by_dimension( meshset, dim, ents );MB_CHK_ERR( rval ); // maybe dimension 2
[ # # ][ # # ]
108 : : }
109 : :
110 : : // Merge Mesh Locally
111 [ # # ]: 0 : if( !skip_local_merge )
112 : : {
113 [ # # ]: 0 : MergeMesh merger( myMB, false );
114 [ # # ]: 0 : merger.merge_entities( ents, myEps );
115 : : // We can return if there is only 1 proc
116 [ # # ][ # # ]: 0 : if( rval != MB_SUCCESS || myPcomm->size() == 1 ) { return rval; }
[ # # ][ # # ]
117 : :
118 : : // Rebuild the ents range
119 [ # # ]: 0 : ents.clear();
120 [ # # ][ # # ]: 0 : rval = myMB->get_entities_by_dimension( meshset, dim, ents );MB_CHK_ERR( rval );
[ # # ][ # # ]
[ # # ]
121 : : }
122 : :
123 : : /*Get Skin
124 : : -Get Range of all dimensional entities
125 : : -skinEnts[i] is the skin entities of dimension i*/
126 [ # # ]: 0 : Skinner skinner( myMB );
127 [ # # ]: 0 : for( int skin_dim = dim; skin_dim >= 0; skin_dim-- )
128 : : {
129 [ # # ][ # # ]: 0 : rval = skinner.find_skin( meshset, ents, skin_dim, mySkinEnts[skin_dim] );MB_CHK_ERR( rval );
[ # # ][ # # ]
[ # # ]
130 : : }
131 : 0 : return MB_SUCCESS;
132 : : }
133 : :
134 : : // Determine the global assembly box
135 : 0 : ErrorCode ParallelMergeMesh::GetGlobalBox( double* gbox )
136 : : {
137 : : ErrorCode rval;
138 : :
139 : : /*Get Bounding Box*/
140 [ # # ]: 0 : BoundBox box;
141 [ # # ][ # # ]: 0 : if( mySkinEnts[0].size() != 0 )
[ # # ]
142 : : {
143 [ # # ][ # # ]: 0 : rval = box.update( *myMB, mySkinEnts[0] );MB_CHK_ERR( rval );
[ # # ][ # # ]
[ # # ]
144 : : }
145 : :
146 : : // Invert the max
147 [ # # ]: 0 : box.bMax *= -1;
148 : :
149 : : /*Communicate to all processors*/
150 [ # # ]: 0 : MPI_Allreduce( (void*)&box, gbox, 6, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
151 : :
152 : : /*Assemble Global Bounding Box*/
153 : : // Flip the max back
154 [ # # ]: 0 : for( int i = 3; i < 6; i++ )
155 : : {
156 : 0 : gbox[i] *= -1;
157 : : }
158 : 0 : return MB_SUCCESS;
159 : : }
160 : :
161 : : // Assemble the tuples with their processor destination
162 : 0 : ErrorCode ParallelMergeMesh::PopulateMyTup( double* gbox )
163 : : {
164 : : /*Figure out how do partition the global box*/
165 : : double lengths[3];
166 : : int parts[3];
167 [ # # ][ # # ]: 0 : ErrorCode rval = PartitionGlobalBox( gbox, lengths, parts );MB_CHK_ERR( rval );
[ # # ][ # # ]
168 : :
169 : : /* Get Skin Coordinates, Vertices */
170 [ # # ][ # # ]: 0 : double* x = new double[mySkinEnts[0].size()];
[ # # ][ # # ]
171 [ # # ][ # # ]: 0 : double* y = new double[mySkinEnts[0].size()];
[ # # ][ # # ]
172 [ # # ][ # # ]: 0 : double* z = new double[mySkinEnts[0].size()];
[ # # ][ # # ]
173 [ # # ][ # # ]: 0 : rval = myMB->get_coords( mySkinEnts[0], x, y, z );
174 [ # # ]: 0 : if( rval != MB_SUCCESS )
175 : : {
176 : : // Prevent Memory Leak
177 [ # # ]: 0 : delete[] x;
178 [ # # ]: 0 : delete[] y;
179 [ # # ]: 0 : delete[] z;
180 : 0 : return rval;
181 : : }
182 : :
183 : : // Initialize variable to be used in the loops
184 [ # # ]: 0 : std::vector< int > toProcs;
185 : : int xPart, yPart, zPart, xEps, yEps, zEps, baseProc;
186 : 0 : unsigned long long tup_i = 0, tup_ul = 0, tup_r = 0, count = 0;
187 : : // These are boolean to determine if the vertex is on close enough to a given border
188 : : bool xDup, yDup, zDup;
189 [ # # ]: 0 : bool canWrite = myTup.get_writeEnabled();
190 [ # # ][ # # ]: 0 : if( !canWrite ) myTup.enableWriteAccess();
191 : : // Go through each vertex
192 [ # # ][ # # ]: 0 : for( Range::iterator it = mySkinEnts[0].begin(); it != mySkinEnts[0].end(); ++it )
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
193 : : {
194 : 0 : xDup = false;
195 : 0 : yDup = false;
196 : 0 : zDup = false;
197 : : // Figure out which x,y,z partition the element is in.
198 : 0 : xPart = static_cast< int >( floor( ( x[count] - gbox[0] ) / lengths[0] ) );
199 [ # # ]: 0 : xPart = ( xPart < parts[0] ? xPart : parts[0] - 1 ); // Make sure it stays within the bounds
200 : :
201 : 0 : yPart = static_cast< int >( floor( ( y[count] - gbox[1] ) / lengths[1] ) );
202 [ # # ]: 0 : yPart = ( yPart < parts[1] ? yPart : parts[1] - 1 ); // Make sure it stays within the bounds
203 : :
204 : 0 : zPart = static_cast< int >( floor( ( z[count] - gbox[2] ) / lengths[2] ) );
205 [ # # ]: 0 : zPart = ( zPart < parts[2] ? zPart : parts[2] - 1 ); // Make sure it stays within the bounds
206 : :
207 : : // Figure out the partition with the addition of Epsilon
208 : 0 : xEps = static_cast< int >( floor( ( x[count] - gbox[0] + myEps ) / lengths[0] ) );
209 : 0 : yEps = static_cast< int >( floor( ( y[count] - gbox[1] + myEps ) / lengths[1] ) );
210 : 0 : zEps = static_cast< int >( floor( ( z[count] - gbox[2] + myEps ) / lengths[2] ) );
211 : :
212 : : // Figure out if the vertex needs to be sent to multiple procs
213 [ # # ][ # # ]: 0 : xDup = ( xPart != xEps && xEps < parts[0] );
214 [ # # ][ # # ]: 0 : yDup = ( yPart != yEps && yEps < parts[1] );
215 [ # # ][ # # ]: 0 : zDup = ( zPart != zEps && zEps < parts[2] );
216 : :
217 : : // Add appropriate processors to the vector
218 : 0 : baseProc = xPart + yPart * parts[0] + zPart * parts[0] * parts[1];
219 [ # # ]: 0 : toProcs.push_back( baseProc );
220 [ # # ]: 0 : if( xDup )
221 : : {
222 [ # # ]: 0 : toProcs.push_back( baseProc + 1 ); // Get partition to the right
223 : : }
224 [ # # ]: 0 : if( yDup )
225 : : {
226 : : // Partition up 1
227 [ # # ]: 0 : toProcs.push_back( baseProc + parts[0] );
228 : : }
229 [ # # ]: 0 : if( zDup )
230 : : {
231 : : // Partition above 1
232 [ # # ]: 0 : toProcs.push_back( baseProc + parts[0] * parts[1] );
233 : : }
234 [ # # ][ # # ]: 0 : if( xDup && yDup )
235 : : {
236 : : // Partition up 1 and right 1
237 [ # # ]: 0 : toProcs.push_back( baseProc + parts[0] + 1 );
238 : : }
239 [ # # ][ # # ]: 0 : if( xDup && zDup )
240 : : {
241 : : // Partition right 1 and above 1
242 [ # # ]: 0 : toProcs.push_back( baseProc + parts[0] * parts[1] + 1 );
243 : : }
244 [ # # ][ # # ]: 0 : if( yDup && zDup )
245 : : {
246 : : // Partition up 1 and above 1
247 [ # # ]: 0 : toProcs.push_back( baseProc + parts[0] * parts[1] + parts[0] );
248 : : }
249 [ # # ][ # # ]: 0 : if( xDup && yDup && zDup )
[ # # ]
250 : : {
251 : : // Partition right 1, up 1, and above 1
252 [ # # ]: 0 : toProcs.push_back( baseProc + parts[0] * parts[1] + parts[0] + 1 );
253 : : }
254 : : // Grow the tuple list if necessary
255 [ # # ][ # # ]: 0 : while( myTup.get_n() + toProcs.size() >= myTup.get_max() )
[ # # ]
256 : : {
257 [ # # ][ # # ]: 0 : myTup.resize( myTup.get_max() ? myTup.get_max() + myTup.get_max() / 2 + 1 : 2 );
[ # # ][ # # ]
[ # # ]
258 : : }
259 : :
260 : : // Add each proc as a tuple
261 [ # # ][ # # ]: 0 : for( std::vector< int >::iterator proc = toProcs.begin(); proc != toProcs.end(); ++proc )
[ # # ]
262 : : {
263 [ # # ]: 0 : myTup.vi_wr[tup_i++] = *proc;
264 [ # # ]: 0 : myTup.vul_wr[tup_ul++] = *it;
265 : 0 : myTup.vr_wr[tup_r++] = x[count];
266 : 0 : myTup.vr_wr[tup_r++] = y[count];
267 : 0 : myTup.vr_wr[tup_r++] = z[count];
268 [ # # ]: 0 : myTup.inc_n();
269 : : }
270 : 0 : count++;
271 : 0 : toProcs.clear();
272 : : }
273 [ # # ]: 0 : delete[] x;
274 [ # # ]: 0 : delete[] y;
275 [ # # ]: 0 : delete[] z;
276 [ # # ][ # # ]: 0 : if( !canWrite ) myTup.disableWriteAccess();
277 : 0 : return MB_SUCCESS;
278 : : }
279 : :
280 : : // Partition the global box by the number of procs
281 : 0 : ErrorCode ParallelMergeMesh::PartitionGlobalBox( double* gbox, double* lengths, int* parts )
282 : : {
283 : : // Determine the length of each side
284 : 0 : double xLen = gbox[3] - gbox[0];
285 : 0 : double yLen = gbox[4] - gbox[1];
286 : 0 : double zLen = gbox[5] - gbox[2];
287 : 0 : unsigned numProcs = myPcomm->size();
288 : :
289 : : // Partition sides from the longest to shortest lengths
290 : : // If x is the longest side
291 [ # # ][ # # ]: 0 : if( xLen >= yLen && xLen >= zLen )
292 : : {
293 : 0 : parts[0] = PartitionSide( xLen, yLen * zLen, numProcs, true );
294 : 0 : numProcs /= parts[0];
295 : : // If y is second longest
296 [ # # ]: 0 : if( yLen >= zLen )
297 : : {
298 : 0 : parts[1] = PartitionSide( yLen, zLen, numProcs, false );
299 : 0 : parts[2] = numProcs / parts[1];
300 : : }
301 : : // If z is the longer
302 : : else
303 : : {
304 : 0 : parts[2] = PartitionSide( zLen, yLen, numProcs, false );
305 : 0 : parts[1] = numProcs / parts[2];
306 : : }
307 : : }
308 : : // If y is the longest side
309 [ # # ]: 0 : else if( yLen >= zLen )
310 : : {
311 : 0 : parts[1] = PartitionSide( yLen, xLen * zLen, numProcs, true );
312 : 0 : numProcs /= parts[1];
313 : : // If x is the second longest
314 [ # # ]: 0 : if( xLen >= zLen )
315 : : {
316 : 0 : parts[0] = PartitionSide( xLen, zLen, numProcs, false );
317 : 0 : parts[2] = numProcs / parts[0];
318 : : }
319 : : // If z is the second longest
320 : : else
321 : : {
322 : 0 : parts[2] = PartitionSide( zLen, xLen, numProcs, false );
323 : 0 : parts[0] = numProcs / parts[2];
324 : : }
325 : : }
326 : : // If z is the longest side
327 : : else
328 : : {
329 : 0 : parts[2] = PartitionSide( zLen, xLen * yLen, numProcs, true );
330 : 0 : numProcs /= parts[2];
331 : : // If x is the second longest
332 [ # # ]: 0 : if( xLen >= yLen )
333 : : {
334 : 0 : parts[0] = PartitionSide( xLen, yLen, numProcs, false );
335 : 0 : parts[1] = numProcs / parts[0];
336 : : }
337 : : // If y is the second longest
338 : : else
339 : : {
340 : 0 : parts[1] = PartitionSide( yLen, xLen, numProcs, false );
341 : 0 : parts[0] = numProcs / parts[1];
342 : : }
343 : : }
344 : :
345 : : // Divide up each side to give the lengths
346 : 0 : lengths[0] = xLen / (double)parts[0];
347 : 0 : lengths[1] = yLen / (double)parts[1];
348 : 0 : lengths[2] = zLen / (double)parts[2];
349 : 0 : return MB_SUCCESS;
350 : : }
351 : :
352 : : // Partition a side based on the length ratios
353 : 0 : int ParallelMergeMesh::PartitionSide( double sideLen, double restLen, unsigned numProcs, bool altRatio )
354 : : {
355 : : // If theres only 1 processor, then just return 1
356 [ # # ]: 0 : if( numProcs == 1 ) { return 1; }
357 : : // Initialize with the ratio of 1 proc
358 : 0 : double ratio = -DBL_MAX;
359 : 0 : unsigned factor = 1;
360 : : // We need to be able to save the last ratio and factor (for comparison)
361 : 0 : double oldRatio = ratio;
362 : 0 : double oldFactor = 1;
363 : :
364 : : // This is the ratio were shooting for
365 : 0 : double goalRatio = sideLen / restLen;
366 : :
367 : : // Calculate the divisor and numerator power
368 : : // This avoid if statements in the loop and is useful since both calculations are similar
369 : : double divisor, p;
370 [ # # ]: 0 : if( altRatio )
371 : : {
372 : 0 : divisor = (double)numProcs * sideLen;
373 : 0 : p = 3;
374 : : }
375 : : else
376 : : {
377 : 0 : divisor = (double)numProcs;
378 : 0 : p = 2;
379 : : }
380 : :
381 : : // Find each possible factor
382 [ # # ]: 0 : for( unsigned i = 2; i <= numProcs / 2; i++ )
383 : : {
384 : : // If it is a factor...
385 [ # # ]: 0 : if( numProcs % i == 0 )
386 : : {
387 : : // We need to save the past factor
388 : 0 : oldRatio = ratio;
389 : 0 : oldFactor = factor;
390 : : // There are 2 different ways to calculate the ratio:
391 : : // Comparing 1 side to 2 sides: (i*i*i)/(numProcs*x)
392 : : // Justification: We have a ratio x:y:z (side Lengths) == a:b:c (procs). So a=kx,
393 : : // b=ky, c=kz. Also, abc=n (numProcs) => bc = n/a. Also, a=kx => k=a/x => 1/k=x/a And so
394 : : // x/(yz) == (kx)/(kyz) == (kx)/(kykz(1/k)) == a/(bc(x/a)) == a/((n/a)(x/a)) == a^3/(nx).
395 : : // Comparing 1 side to 1 side: (i*i)/numprocs
396 : : // Justification: i/(n/i) == i^2/n
397 : 0 : ratio = pow( (double)i, p ) / divisor;
398 : 0 : factor = i;
399 : : // Once we have passed the goal ratio, we can break since we'll only move away from the
400 : : // goal ratio
401 [ # # ]: 0 : if( ratio >= goalRatio ) { break; }
402 : : }
403 : : }
404 : : // If we haven't reached the goal ratio yet, check out factor = numProcs
405 [ # # ]: 0 : if( ratio < goalRatio )
406 : : {
407 : 0 : oldRatio = ratio;
408 : 0 : oldFactor = factor;
409 : 0 : factor = numProcs;
410 : 0 : ratio = pow( (double)numProcs, p ) / divisor;
411 : : }
412 : :
413 : : // Figure out if our oldRatio is better than ratio
414 [ # # ]: 0 : if( fabs( ratio - goalRatio ) > fabs( oldRatio - goalRatio ) ) { factor = oldFactor; }
415 : : // Return our findings
416 : 0 : return factor;
417 : : }
418 : :
419 : : // Id the tuples that are matching
420 : 0 : ErrorCode ParallelMergeMesh::PopulateMyMatches()
421 : : {
422 : : // Counters for accessing tuples more efficiently
423 : 0 : unsigned long i = 0, mat_i = 0, mat_ul = 0, j = 0, tup_r = 0;
424 : 0 : double eps2 = myEps * myEps;
425 : :
426 : : uint tup_mi, tup_ml, tup_mul, tup_mr;
427 [ # # ]: 0 : myTup.getTupleSize( tup_mi, tup_ml, tup_mul, tup_mr );
428 : :
429 [ # # ]: 0 : bool canWrite = myMatches.get_writeEnabled();
430 [ # # ][ # # ]: 0 : if( !canWrite ) myMatches.enableWriteAccess();
431 : :
432 [ # # ][ # # ]: 0 : while( ( i + 1 ) < myTup.get_n() )
433 : : {
434 : : // Proximity Comparison
435 : 0 : double xi = myTup.vr_rd[tup_r], yi = myTup.vr_rd[tup_r + 1], zi = myTup.vr_rd[tup_r + 2];
436 : :
437 : 0 : bool done = false;
438 [ # # ]: 0 : while( !done )
439 : : {
440 : 0 : j++;
441 : 0 : tup_r += tup_mr;
442 [ # # ][ # # ]: 0 : if( j >= myTup.get_n() ) { break; }
443 [ # # ]: 0 : CartVect cv( myTup.vr_rd[tup_r] - xi, myTup.vr_rd[tup_r + 1] - yi, myTup.vr_rd[tup_r + 2] - zi );
444 [ # # ][ # # ]: 0 : if( cv.length_squared() > eps2 ) { done = true; }
445 : : }
446 : : // Allocate the tuple list before adding matches
447 [ # # ][ # # ]: 0 : while( myMatches.get_n() + ( j - i ) * ( j - i - 1 ) >= myMatches.get_max() )
[ # # ]
448 : : {
449 [ # # ][ # # ]: 0 : myMatches.resize( myMatches.get_max() ? myMatches.get_max() + myMatches.get_max() / 2 + 1 : 2 );
[ # # ][ # # ]
[ # # ]
450 : : }
451 : :
452 : : // We now know that tuples [i to j) exclusive match.
453 : : // If n tuples match, n*(n-1) match tuples will be made
454 : : // tuples are of the form (proc1,proc2,handle1,handle2)
455 [ # # ]: 0 : if( i + 1 < j )
456 : : {
457 : 0 : int kproc = i * tup_mi;
458 : 0 : unsigned long khand = i * tup_mul;
459 [ # # ]: 0 : for( unsigned long k = i; k < j; k++ )
460 : : {
461 : 0 : int lproc = kproc + tup_mi;
462 : 0 : unsigned long lhand = khand + tup_mul;
463 [ # # ]: 0 : for( unsigned long l = k + 1; l < j; l++ )
464 : : {
465 : 0 : myMatches.vi_wr[mat_i++] = myTup.vi_rd[kproc]; // proc1
466 : 0 : myMatches.vi_wr[mat_i++] = myTup.vi_rd[lproc]; // proc2
467 : 0 : myMatches.vul_wr[mat_ul++] = myTup.vul_rd[khand]; // handle1
468 : 0 : myMatches.vul_wr[mat_ul++] = myTup.vul_rd[lhand]; // handle2
469 [ # # ]: 0 : myMatches.inc_n();
470 : :
471 : 0 : myMatches.vi_wr[mat_i++] = myTup.vi_rd[lproc]; // proc1
472 : 0 : myMatches.vi_wr[mat_i++] = myTup.vi_rd[kproc]; // proc2
473 : 0 : myMatches.vul_wr[mat_ul++] = myTup.vul_rd[lhand]; // handle1
474 : 0 : myMatches.vul_wr[mat_ul++] = myTup.vul_rd[khand]; // handle2
475 [ # # ]: 0 : myMatches.inc_n();
476 : 0 : lproc += tup_mi;
477 : 0 : lhand += tup_mul;
478 : : }
479 : 0 : kproc += tup_mi;
480 : 0 : khand += tup_mul;
481 : : } // End for(int k...
482 : : }
483 : 0 : i = j;
484 : : } // End while(i+1<tup.n)
485 : :
486 [ # # ][ # # ]: 0 : if( !canWrite ) myMatches.disableWriteAccess();
487 : 0 : return MB_SUCCESS;
488 : : }
489 : :
490 : : // Sort the matching tuples so that vertices can be tagged accurately
491 : 0 : ErrorCode ParallelMergeMesh::SortMyMatches()
492 : : {
493 [ # # ][ # # ]: 0 : TupleList::buffer buf( mySkinEnts[0].size() );
[ # # ]
494 : : // Sorts are necessary to check for doubles
495 : : // Sort by remote handle
496 [ # # ]: 0 : myMatches.sort( 3, &buf );
497 : : // Sort by matching proc
498 [ # # ]: 0 : myMatches.sort( 1, &buf );
499 : : // Sort by local handle
500 [ # # ]: 0 : myMatches.sort( 2, &buf );
501 [ # # ]: 0 : buf.reset();
502 : 0 : return MB_SUCCESS;
503 : : }
504 : :
505 : : // Tag the shared elements using existing PComm functionality
506 : 0 : ErrorCode ParallelMergeMesh::TagSharedElements( int dim )
507 : : {
508 : : // Manipulate the matches list to tag vertices and entities
509 : : // Set up proc ents
510 [ # # ]: 0 : Range proc_ents;
511 : : ErrorCode rval;
512 : :
513 : : // get the entities in the partition sets
514 [ # # ][ # # ]: 0 : for( Range::iterator rit = myPcomm->partitionSets.begin(); rit != myPcomm->partitionSets.end(); ++rit )
[ # # ][ # # ]
[ # # ]
515 : : {
516 [ # # ]: 0 : Range tmp_ents;
517 [ # # ][ # # ]: 0 : rval = myMB->get_entities_by_handle( *rit, tmp_ents, true );
518 [ # # ]: 0 : if( MB_SUCCESS != rval ) { return rval; }
519 [ # # ][ # # ]: 0 : proc_ents.merge( tmp_ents );
520 : 0 : }
521 [ # # ][ # # ]: 0 : if( myMB->dimension_from_handle( *proc_ents.rbegin() ) != myMB->dimension_from_handle( *proc_ents.begin() ) )
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
522 : : {
523 [ # # ]: 0 : Range::iterator lower = proc_ents.lower_bound( CN::TypeDimensionMap[0].first ),
524 [ # # ]: 0 : upper = proc_ents.upper_bound( CN::TypeDimensionMap[dim - 1].second );
525 [ # # ]: 0 : proc_ents.erase( lower, upper );
526 : : }
527 : :
528 : : // This vector doesn't appear to be used but its in resolve_shared_ents
529 : 0 : int maxp = -1;
530 [ # # ]: 0 : std::vector< int > sharing_procs( MAX_SHARING_PROCS );
531 [ # # ]: 0 : std::fill( sharing_procs.begin(), sharing_procs.end(), maxp );
532 : :
533 : : // get ents shared by 1 or n procs
534 [ # # ]: 0 : std::map< std::vector< int >, std::vector< EntityHandle > > proc_nranges;
535 [ # # ]: 0 : Range proc_verts;
536 [ # # ]: 0 : rval = myMB->get_adjacencies( proc_ents, 0, false, proc_verts, Interface::UNION );
537 [ # # ]: 0 : if( rval != MB_SUCCESS ) { return rval; }
538 : :
539 [ # # ]: 0 : rval = myPcomm->tag_shared_verts( myMatches, proc_nranges, proc_verts );
540 [ # # ]: 0 : if( rval != MB_SUCCESS ) { return rval; }
541 : :
542 : : // get entities shared by 1 or n procs
543 [ # # ][ # # ]: 0 : rval = myPcomm->get_proc_nvecs( dim, dim - 1, &mySkinEnts[0], proc_nranges );
544 [ # # ]: 0 : if( rval != MB_SUCCESS ) { return rval; }
545 : :
546 : : // create the sets for each interface; store them as tags on
547 : : // the interface instance
548 [ # # ]: 0 : Range iface_sets;
549 [ # # ]: 0 : rval = myPcomm->create_interface_sets( proc_nranges );
550 [ # # ]: 0 : if( rval != MB_SUCCESS ) { return rval; }
551 : : // establish comm procs and buffers for them
552 [ # # ]: 0 : std::set< unsigned int > procs;
553 [ # # ]: 0 : rval = myPcomm->get_interface_procs( procs, true );
554 [ # # ]: 0 : if( rval != MB_SUCCESS ) { return rval; }
555 : :
556 : : // resolve shared entity remote handles; implemented in ghost cell exchange
557 : : // code because it's so similar
558 [ # # ]: 0 : rval = myPcomm->exchange_ghost_cells( -1, -1, 0, true, true );
559 [ # # ]: 0 : if( rval != MB_SUCCESS ) { return rval; }
560 : : // now build parent/child links for interface sets
561 [ # # ]: 0 : rval = myPcomm->create_iface_pc_links();
562 : 0 : return rval;
563 : : }
564 : :
565 : : // Make sure to free up any allocated data
566 : : // Need to avoid a double free
567 : 0 : void ParallelMergeMesh::CleanUp()
568 : : {
569 : : // The reset operation is now safe and avoids a double free()
570 : 0 : myMatches.reset();
571 : 0 : myTup.reset();
572 : 0 : myCD.reset();
573 : 0 : }
574 : :
575 : : // Simple quick sort to real
576 : 0 : void ParallelMergeMesh::SortTuplesByReal( TupleList& tup, double eps )
577 : : {
578 [ # # ]: 0 : bool canWrite = tup.get_writeEnabled();
579 [ # # ][ # # ]: 0 : if( !canWrite ) tup.enableWriteAccess();
580 : :
581 : : uint mi, ml, mul, mr;
582 [ # # ]: 0 : tup.getTupleSize( mi, ml, mul, mr );
583 [ # # ][ # # ]: 0 : PerformRealSort( tup, 0, tup.get_n(), eps, mr );
584 : :
585 [ # # ][ # # ]: 0 : if( !canWrite ) tup.disableWriteAccess();
586 : 0 : }
587 : :
588 : : // Swap around tuples
589 : 0 : void ParallelMergeMesh::SwapTuples( TupleList& tup, unsigned long a, unsigned long b )
590 : : {
591 [ # # ]: 0 : if( a == b ) return;
592 : :
593 : : uint mi, ml, mul, mr;
594 [ # # ]: 0 : tup.getTupleSize( mi, ml, mul, mr );
595 : :
596 : : // Swap mi
597 : 0 : unsigned long a_val = a * mi, b_val = b * mi;
598 [ # # ]: 0 : for( unsigned long i = 0; i < mi; i++ )
599 : : {
600 : 0 : sint t = tup.vi_rd[a_val];
601 : 0 : tup.vi_wr[a_val] = tup.vi_rd[b_val];
602 : 0 : tup.vi_wr[b_val] = t;
603 : 0 : a_val++;
604 : 0 : b_val++;
605 : : }
606 : : // Swap ml
607 : 0 : a_val = a * ml;
608 : 0 : b_val = b * ml;
609 [ # # ]: 0 : for( unsigned long i = 0; i < ml; i++ )
610 : : {
611 : 0 : slong t = tup.vl_rd[a_val];
612 : 0 : tup.vl_wr[a_val] = tup.vl_rd[b_val];
613 : 0 : tup.vl_wr[b_val] = t;
614 : 0 : a_val++;
615 : 0 : b_val++;
616 : : }
617 : : // Swap mul
618 : 0 : a_val = a * mul;
619 : 0 : b_val = b * mul;
620 [ # # ]: 0 : for( unsigned long i = 0; i < mul; i++ )
621 : : {
622 : 0 : Ulong t = tup.vul_rd[a_val];
623 : 0 : tup.vul_wr[a_val] = tup.vul_rd[b_val];
624 : 0 : tup.vul_wr[b_val] = t;
625 : 0 : a_val++;
626 : 0 : b_val++;
627 : : }
628 : : // Swap mr
629 : 0 : a_val = a * mr;
630 : 0 : b_val = b * mr;
631 [ # # ]: 0 : for( unsigned long i = 0; i < mr; i++ )
632 : : {
633 : 0 : realType t = tup.vr_rd[a_val];
634 : 0 : tup.vr_wr[a_val] = tup.vr_rd[b_val];
635 : 0 : tup.vr_wr[b_val] = t;
636 : 0 : a_val++;
637 : 0 : b_val++;
638 : : }
639 : : }
640 : :
641 : : // Perform the sorting of a tuple by real
642 : : // To sort an entire tuple_list, call (tup,0,tup.n,epsilon)
643 : 0 : void ParallelMergeMesh::PerformRealSort( TupleList& tup, unsigned long left, unsigned long right, double eps,
644 : : uint tup_mr )
645 : : {
646 : : // If list size is only 1 or 0 return
647 [ # # ]: 0 : if( left + 1 >= right ) { return; }
648 : 0 : unsigned long swap = left, tup_l = left * tup_mr, tup_t = tup_l + tup_mr;
649 : :
650 : : // Swap the median with the left position for a (hopefully) better split
651 : 0 : SwapTuples( tup, left, ( left + right ) / 2 );
652 : :
653 : : // Partition the data
654 [ # # ]: 0 : for( unsigned long t = left + 1; t < right; t++ )
655 : : {
656 : : // If the left value(pivot) is greater than t_val, swap it into swap
657 [ # # ]: 0 : if( TupleGreaterThan( tup, tup_l, tup_t, eps, tup_mr ) )
658 : : {
659 : 0 : swap++;
660 : 0 : SwapTuples( tup, swap, t );
661 : : }
662 : 0 : tup_t += tup_mr;
663 : : }
664 : :
665 : : // Swap so that position swap is in the correct position
666 : 0 : SwapTuples( tup, left, swap );
667 : :
668 : : // Sort left and right of swap
669 : 0 : PerformRealSort( tup, left, swap, eps, tup_mr );
670 : 0 : PerformRealSort( tup, swap + 1, right, eps, tup_mr );
671 : : }
672 : :
673 : : // Note, this takes the actual tup.vr[] index (aka i*tup.mr)
674 : 0 : bool ParallelMergeMesh::TupleGreaterThan( TupleList& tup, unsigned long vrI, unsigned long vrJ, double eps,
675 : : uint tup_mr )
676 : : {
677 : 0 : unsigned check = 0;
678 [ # # ]: 0 : while( check < tup_mr )
679 : : {
680 : : // If the values are the same
681 [ # # ]: 0 : if( fabs( tup.vr_rd[vrI + check] - tup.vr_rd[vrJ + check] ) <= eps )
682 : : {
683 : 0 : check++;
684 : 0 : continue;
685 : : }
686 : : // If I greater than J
687 [ # # ]: 0 : else if( tup.vr_rd[vrI + check] > tup.vr_rd[vrJ + check] )
688 : : {
689 : 0 : return true;
690 : : }
691 : : // If J greater than I
692 : : else
693 : : {
694 : 0 : return false;
695 : : }
696 : : }
697 : : // All Values are the same
698 : 0 : return false;
699 : : }
700 : :
701 [ + - ][ + - ]: 228 : } // End namespace moab
|
