MOAB: Mesh Oriented datABase
(version 5.4.1)
|
00001 #include "moab/ParallelMergeMesh.hpp" 00002 #include "moab/Core.hpp" 00003 #include "moab/CartVect.hpp" 00004 #include "moab/BoundBox.hpp" 00005 #include "moab/Skinner.hpp" 00006 #include "moab/MergeMesh.hpp" 00007 #include "moab/CN.hpp" 00008 #include <cfloat> 00009 #include <algorithm> 00010 00011 #ifdef MOAB_HAVE_MPI 00012 #include "moab_mpi.h" 00013 #endif 00014 00015 namespace moab 00016 { 00017 00018 // Constructor 00019 /*Get Merge Data and tolerance*/ 00020 ParallelMergeMesh::ParallelMergeMesh( ParallelComm* pc, const double epsilon ) : myPcomm( pc ), myEps( epsilon ) 00021 { 00022 myMB = pc->get_moab(); 00023 mySkinEnts.resize( 4 ); 00024 } 00025 00026 // Have a wrapper function on the actual merge to avoid memory leaks 00027 // Merges elements within a proximity of epsilon 00028 ErrorCode ParallelMergeMesh::merge( EntityHandle levelset, bool skip_local_merge, int dim ) 00029 { 00030 ErrorCode rval = PerformMerge( levelset, skip_local_merge, dim );MB_CHK_ERR( rval ); 00031 CleanUp(); 00032 return rval; 00033 } 00034 00035 // Perform the merge 00036 ErrorCode ParallelMergeMesh::PerformMerge( EntityHandle levelset, bool skip_local_merge, int dim ) 00037 { 00038 // Get the mesh dimension 00039 ErrorCode rval; 00040 if( dim < 0 ) 00041 { 00042 rval = myMB->get_dimension( dim );MB_CHK_ERR( rval ); 00043 } 00044 00045 // Get the local skin elements 00046 rval = PopulateMySkinEnts( levelset, dim, skip_local_merge ); 00047 // If there is only 1 proc, we can return now 00048 if( rval != MB_SUCCESS || myPcomm->size() == 1 ) 00049 { 00050 return rval; 00051 } 00052 00053 // Determine the global bounding box 00054 double gbox[6]; 00055 rval = GetGlobalBox( gbox );MB_CHK_ERR( rval ); 00056 00057 /* Assemble The Destination Tuples */ 00058 // Get a list of tuples which contain (toProc, handle, x,y,z) 00059 myTup.initialize( 1, 0, 1, 3, mySkinEnts[0].size() ); 00060 rval = PopulateMyTup( gbox );MB_CHK_ERR( rval ); 00061 00062 /* Gather-Scatter Tuple 00063 -tup comes out as (remoteProc,handle,x,y,z) */ 00064 myCD.initialize( myPcomm->comm() ); 00065 00066 // 1 represents dynamic tuple, 0 represents index of the processor to send to 00067 myCD.gs_transfer( 1, myTup, 0 ); 00068 00069 /* Sort By X,Y,Z 00070 -Utilizes a custom quick sort incorporating epsilon*/ 00071 SortTuplesByReal( myTup, myEps ); 00072 00073 // Initialize another tuple list for matches 00074 myMatches.initialize( 2, 0, 2, 0, mySkinEnts[0].size() ); 00075 00076 // ID the matching tuples 00077 rval = PopulateMyMatches();MB_CHK_ERR( rval ); 00078 00079 // We can free up the tuple myTup now 00080 myTup.reset(); 00081 00082 /*Gather-Scatter Again*/ 00083 // 1 represents dynamic list, 0 represents proc index to send tuple to 00084 myCD.gs_transfer( 1, myMatches, 0 ); 00085 // We can free up the crystal router now 00086 myCD.reset(); 00087 00088 // Sort the matches tuple list 00089 SortMyMatches(); 00090 00091 // Tag the shared elements 00092 rval = TagSharedElements( dim );MB_CHK_ERR( rval ); 00093 00094 // Free up the matches tuples 00095 myMatches.reset(); 00096 return rval; 00097 } 00098 00099 // Sets mySkinEnts with all of the skin entities on the processor 00100 ErrorCode ParallelMergeMesh::PopulateMySkinEnts( const EntityHandle meshset, int dim, bool skip_local_merge ) 00101 { 00102 /*Merge Mesh Locally*/ 00103 // Get all dim dimensional entities 00104 Range ents; 00105 ErrorCode rval = myMB->get_entities_by_dimension( meshset, dim, ents );MB_CHK_ERR( rval ); 00106 00107 if( ents.empty() && dim == 3 ) 00108 { 00109 dim--; 00110 rval = myMB->get_entities_by_dimension( meshset, dim, ents );MB_CHK_ERR( rval ); // maybe dimension 2 00111 } 00112 00113 // Merge Mesh Locally 00114 if( !skip_local_merge ) 00115 { 00116 MergeMesh merger( myMB, false ); 00117 merger.merge_entities( ents, myEps ); 00118 // We can return if there is only 1 proc 00119 if( rval != MB_SUCCESS || myPcomm->size() == 1 ) 00120 { 00121 return rval; 00122 } 00123 00124 // Rebuild the ents range 00125 ents.clear(); 00126 rval = myMB->get_entities_by_dimension( meshset, dim, ents );MB_CHK_ERR( rval ); 00127 } 00128 00129 /*Get Skin 00130 -Get Range of all dimensional entities 00131 -skinEnts[i] is the skin entities of dimension i*/ 00132 Skinner skinner( myMB ); 00133 for( int skin_dim = dim; skin_dim >= 0; skin_dim-- ) 00134 { 00135 rval = skinner.find_skin( meshset, ents, skin_dim, mySkinEnts[skin_dim] );MB_CHK_ERR( rval ); 00136 } 00137 return MB_SUCCESS; 00138 } 00139 00140 // Determine the global assembly box 00141 ErrorCode ParallelMergeMesh::GetGlobalBox( double* gbox ) 00142 { 00143 ErrorCode rval; 00144 00145 /*Get Bounding Box*/ 00146 BoundBox box; 00147 if( mySkinEnts[0].size() != 0 ) 00148 { 00149 rval = box.update( *myMB, mySkinEnts[0] );MB_CHK_ERR( rval ); 00150 } 00151 00152 // Invert the max 00153 box.bMax *= -1; 00154 00155 /*Communicate to all processors*/ 00156 MPI_Allreduce( (void*)&box, gbox, 6, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD ); 00157 00158 /*Assemble Global Bounding Box*/ 00159 // Flip the max back 00160 for( int i = 3; i < 6; i++ ) 00161 { 00162 gbox[i] *= -1; 00163 } 00164 return MB_SUCCESS; 00165 } 00166 00167 // Assemble the tuples with their processor destination 00168 ErrorCode ParallelMergeMesh::PopulateMyTup( double* gbox ) 00169 { 00170 /*Figure out how do partition the global box*/ 00171 double lengths[3]; 00172 int parts[3]; 00173 ErrorCode rval = PartitionGlobalBox( gbox, lengths, parts );MB_CHK_ERR( rval ); 00174 00175 /* Get Skin Coordinates, Vertices */ 00176 double* x = new double[mySkinEnts[0].size()]; 00177 double* y = new double[mySkinEnts[0].size()]; 00178 double* z = new double[mySkinEnts[0].size()]; 00179 rval = myMB->get_coords( mySkinEnts[0], x, y, z ); 00180 if( rval != MB_SUCCESS ) 00181 { 00182 // Prevent Memory Leak 00183 delete[] x; 00184 delete[] y; 00185 delete[] z; 00186 return rval; 00187 } 00188 00189 // Initialize variable to be used in the loops 00190 std::vector< int > toProcs; 00191 int xPart, yPart, zPart, xEps, yEps, zEps, baseProc; 00192 unsigned long long tup_i = 0, tup_ul = 0, tup_r = 0, count = 0; 00193 // These are boolean to determine if the vertex is on close enough to a given border 00194 bool xDup, yDup, zDup; 00195 bool canWrite = myTup.get_writeEnabled(); 00196 if( !canWrite ) myTup.enableWriteAccess(); 00197 // Go through each vertex 00198 for( Range::iterator it = mySkinEnts[0].begin(); it != mySkinEnts[0].end(); ++it ) 00199 { 00200 xDup = false; 00201 yDup = false; 00202 zDup = false; 00203 // Figure out which x,y,z partition the element is in. 00204 xPart = static_cast< int >( floor( ( x[count] - gbox[0] ) / lengths[0] ) ); 00205 xPart = ( xPart < parts[0] ? xPart : parts[0] - 1 ); // Make sure it stays within the bounds 00206 00207 yPart = static_cast< int >( floor( ( y[count] - gbox[1] ) / lengths[1] ) ); 00208 yPart = ( yPart < parts[1] ? yPart : parts[1] - 1 ); // Make sure it stays within the bounds 00209 00210 zPart = static_cast< int >( floor( ( z[count] - gbox[2] ) / lengths[2] ) ); 00211 zPart = ( zPart < parts[2] ? zPart : parts[2] - 1 ); // Make sure it stays within the bounds 00212 00213 // Figure out the partition with the addition of Epsilon 00214 xEps = static_cast< int >( floor( ( x[count] - gbox[0] + myEps ) / lengths[0] ) ); 00215 yEps = static_cast< int >( floor( ( y[count] - gbox[1] + myEps ) / lengths[1] ) ); 00216 zEps = static_cast< int >( floor( ( z[count] - gbox[2] + myEps ) / lengths[2] ) ); 00217 00218 // Figure out if the vertex needs to be sent to multiple procs 00219 xDup = ( xPart != xEps && xEps < parts[0] ); 00220 yDup = ( yPart != yEps && yEps < parts[1] ); 00221 zDup = ( zPart != zEps && zEps < parts[2] ); 00222 00223 // Add appropriate processors to the vector 00224 baseProc = xPart + yPart * parts[0] + zPart * parts[0] * parts[1]; 00225 toProcs.push_back( baseProc ); 00226 if( xDup ) 00227 { 00228 toProcs.push_back( baseProc + 1 ); // Get partition to the right 00229 } 00230 if( yDup ) 00231 { 00232 // Partition up 1 00233 toProcs.push_back( baseProc + parts[0] ); 00234 } 00235 if( zDup ) 00236 { 00237 // Partition above 1 00238 toProcs.push_back( baseProc + parts[0] * parts[1] ); 00239 } 00240 if( xDup && yDup ) 00241 { 00242 // Partition up 1 and right 1 00243 toProcs.push_back( baseProc + parts[0] + 1 ); 00244 } 00245 if( xDup && zDup ) 00246 { 00247 // Partition right 1 and above 1 00248 toProcs.push_back( baseProc + parts[0] * parts[1] + 1 ); 00249 } 00250 if( yDup && zDup ) 00251 { 00252 // Partition up 1 and above 1 00253 toProcs.push_back( baseProc + parts[0] * parts[1] + parts[0] ); 00254 } 00255 if( xDup && yDup && zDup ) 00256 { 00257 // Partition right 1, up 1, and above 1 00258 toProcs.push_back( baseProc + parts[0] * parts[1] + parts[0] + 1 ); 00259 } 00260 // Grow the tuple list if necessary 00261 while( myTup.get_n() + toProcs.size() >= myTup.get_max() ) 00262 { 00263 myTup.resize( myTup.get_max() ? myTup.get_max() + myTup.get_max() / 2 + 1 : 2 ); 00264 } 00265 00266 // Add each proc as a tuple 00267 for( std::vector< int >::iterator proc = toProcs.begin(); proc != toProcs.end(); ++proc ) 00268 { 00269 myTup.vi_wr[tup_i++] = *proc; 00270 myTup.vul_wr[tup_ul++] = *it; 00271 myTup.vr_wr[tup_r++] = x[count]; 00272 myTup.vr_wr[tup_r++] = y[count]; 00273 myTup.vr_wr[tup_r++] = z[count]; 00274 myTup.inc_n(); 00275 } 00276 count++; 00277 toProcs.clear(); 00278 } 00279 delete[] x; 00280 delete[] y; 00281 delete[] z; 00282 if( !canWrite ) myTup.disableWriteAccess(); 00283 return MB_SUCCESS; 00284 } 00285 00286 // Partition the global box by the number of procs 00287 ErrorCode ParallelMergeMesh::PartitionGlobalBox( double* gbox, double* lengths, int* parts ) 00288 { 00289 // Determine the length of each side 00290 double xLen = gbox[3] - gbox[0]; 00291 double yLen = gbox[4] - gbox[1]; 00292 double zLen = gbox[5] - gbox[2]; 00293 unsigned numProcs = myPcomm->size(); 00294 00295 // Partition sides from the longest to shortest lengths 00296 // If x is the longest side 00297 if( xLen >= yLen && xLen >= zLen ) 00298 { 00299 parts[0] = PartitionSide( xLen, yLen * zLen, numProcs, true ); 00300 numProcs /= parts[0]; 00301 // If y is second longest 00302 if( yLen >= zLen ) 00303 { 00304 parts[1] = PartitionSide( yLen, zLen, numProcs, false ); 00305 parts[2] = numProcs / parts[1]; 00306 } 00307 // If z is the longer 00308 else 00309 { 00310 parts[2] = PartitionSide( zLen, yLen, numProcs, false ); 00311 parts[1] = numProcs / parts[2]; 00312 } 00313 } 00314 // If y is the longest side 00315 else if( yLen >= zLen ) 00316 { 00317 parts[1] = PartitionSide( yLen, xLen * zLen, numProcs, true ); 00318 numProcs /= parts[1]; 00319 // If x is the second longest 00320 if( xLen >= zLen ) 00321 { 00322 parts[0] = PartitionSide( xLen, zLen, numProcs, false ); 00323 parts[2] = numProcs / parts[0]; 00324 } 00325 // If z is the second longest 00326 else 00327 { 00328 parts[2] = PartitionSide( zLen, xLen, numProcs, false ); 00329 parts[0] = numProcs / parts[2]; 00330 } 00331 } 00332 // If z is the longest side 00333 else 00334 { 00335 parts[2] = PartitionSide( zLen, xLen * yLen, numProcs, true ); 00336 numProcs /= parts[2]; 00337 // If x is the second longest 00338 if( xLen >= yLen ) 00339 { 00340 parts[0] = PartitionSide( xLen, yLen, numProcs, false ); 00341 parts[1] = numProcs / parts[0]; 00342 } 00343 // If y is the second longest 00344 else 00345 { 00346 parts[1] = PartitionSide( yLen, xLen, numProcs, false ); 00347 parts[0] = numProcs / parts[1]; 00348 } 00349 } 00350 00351 // Divide up each side to give the lengths 00352 lengths[0] = xLen / (double)parts[0]; 00353 lengths[1] = yLen / (double)parts[1]; 00354 lengths[2] = zLen / (double)parts[2]; 00355 return MB_SUCCESS; 00356 } 00357 00358 // Partition a side based on the length ratios 00359 int ParallelMergeMesh::PartitionSide( double sideLen, double restLen, unsigned numProcs, bool altRatio ) 00360 { 00361 // If theres only 1 processor, then just return 1 00362 if( numProcs == 1 ) 00363 { 00364 return 1; 00365 } 00366 // Initialize with the ratio of 1 proc 00367 double ratio = -DBL_MAX; 00368 unsigned factor = 1; 00369 // We need to be able to save the last ratio and factor (for comparison) 00370 double oldRatio = ratio; 00371 double oldFactor = 1; 00372 00373 // This is the ratio were shooting for 00374 double goalRatio = sideLen / restLen; 00375 00376 // Calculate the divisor and numerator power 00377 // This avoid if statements in the loop and is useful since both calculations are similar 00378 double divisor, p; 00379 if( altRatio ) 00380 { 00381 divisor = (double)numProcs * sideLen; 00382 p = 3; 00383 } 00384 else 00385 { 00386 divisor = (double)numProcs; 00387 p = 2; 00388 } 00389 00390 // Find each possible factor 00391 for( unsigned i = 2; i <= numProcs / 2; i++ ) 00392 { 00393 // If it is a factor... 00394 if( numProcs % i == 0 ) 00395 { 00396 // We need to save the past factor 00397 oldRatio = ratio; 00398 oldFactor = factor; 00399 // There are 2 different ways to calculate the ratio: 00400 // Comparing 1 side to 2 sides: (i*i*i)/(numProcs*x) 00401 // Justification: We have a ratio x:y:z (side Lengths) == a:b:c (procs). So a=kx, 00402 // b=ky, c=kz. Also, abc=n (numProcs) => bc = n/a. Also, a=kx => k=a/x => 1/k=x/a And so 00403 // x/(yz) == (kx)/(kyz) == (kx)/(kykz(1/k)) == a/(bc(x/a)) == a/((n/a)(x/a)) == a^3/(nx). 00404 // Comparing 1 side to 1 side: (i*i)/numprocs 00405 // Justification: i/(n/i) == i^2/n 00406 ratio = pow( (double)i, p ) / divisor; 00407 factor = i; 00408 // Once we have passed the goal ratio, we can break since we'll only move away from the 00409 // goal ratio 00410 if( ratio >= goalRatio ) 00411 { 00412 break; 00413 } 00414 } 00415 } 00416 // If we haven't reached the goal ratio yet, check out factor = numProcs 00417 if( ratio < goalRatio ) 00418 { 00419 oldRatio = ratio; 00420 oldFactor = factor; 00421 factor = numProcs; 00422 ratio = pow( (double)numProcs, p ) / divisor; 00423 } 00424 00425 // Figure out if our oldRatio is better than ratio 00426 if( fabs( ratio - goalRatio ) > fabs( oldRatio - goalRatio ) ) 00427 { 00428 factor = oldFactor; 00429 } 00430 // Return our findings 00431 return factor; 00432 } 00433 00434 // Id the tuples that are matching 00435 ErrorCode ParallelMergeMesh::PopulateMyMatches() 00436 { 00437 // Counters for accessing tuples more efficiently 00438 unsigned long i = 0, mat_i = 0, mat_ul = 0, j = 0, tup_r = 0; 00439 double eps2 = myEps * myEps; 00440 00441 uint tup_mi, tup_ml, tup_mul, tup_mr; 00442 myTup.getTupleSize( tup_mi, tup_ml, tup_mul, tup_mr ); 00443 00444 bool canWrite = myMatches.get_writeEnabled(); 00445 if( !canWrite ) myMatches.enableWriteAccess(); 00446 00447 while( ( i + 1 ) < myTup.get_n() ) 00448 { 00449 // Proximity Comparison 00450 double xi = myTup.vr_rd[tup_r], yi = myTup.vr_rd[tup_r + 1], zi = myTup.vr_rd[tup_r + 2]; 00451 00452 bool done = false; 00453 while( !done ) 00454 { 00455 j++; 00456 tup_r += tup_mr; 00457 if( j >= myTup.get_n() ) 00458 { 00459 break; 00460 } 00461 CartVect cv( myTup.vr_rd[tup_r] - xi, myTup.vr_rd[tup_r + 1] - yi, myTup.vr_rd[tup_r + 2] - zi ); 00462 if( cv.length_squared() > eps2 ) 00463 { 00464 done = true; 00465 } 00466 } 00467 // Allocate the tuple list before adding matches 00468 while( myMatches.get_n() + ( j - i ) * ( j - i - 1 ) >= myMatches.get_max() ) 00469 { 00470 myMatches.resize( myMatches.get_max() ? myMatches.get_max() + myMatches.get_max() / 2 + 1 : 2 ); 00471 } 00472 00473 // We now know that tuples [i to j) exclusive match. 00474 // If n tuples match, n*(n-1) match tuples will be made 00475 // tuples are of the form (proc1,proc2,handle1,handle2) 00476 if( i + 1 < j ) 00477 { 00478 int kproc = i * tup_mi; 00479 unsigned long khand = i * tup_mul; 00480 for( unsigned long k = i; k < j; k++ ) 00481 { 00482 int lproc = kproc + tup_mi; 00483 unsigned long lhand = khand + tup_mul; 00484 for( unsigned long l = k + 1; l < j; l++ ) 00485 { 00486 myMatches.vi_wr[mat_i++] = myTup.vi_rd[kproc]; // proc1 00487 myMatches.vi_wr[mat_i++] = myTup.vi_rd[lproc]; // proc2 00488 myMatches.vul_wr[mat_ul++] = myTup.vul_rd[khand]; // handle1 00489 myMatches.vul_wr[mat_ul++] = myTup.vul_rd[lhand]; // handle2 00490 myMatches.inc_n(); 00491 00492 myMatches.vi_wr[mat_i++] = myTup.vi_rd[lproc]; // proc1 00493 myMatches.vi_wr[mat_i++] = myTup.vi_rd[kproc]; // proc2 00494 myMatches.vul_wr[mat_ul++] = myTup.vul_rd[lhand]; // handle1 00495 myMatches.vul_wr[mat_ul++] = myTup.vul_rd[khand]; // handle2 00496 myMatches.inc_n(); 00497 lproc += tup_mi; 00498 lhand += tup_mul; 00499 } 00500 kproc += tup_mi; 00501 khand += tup_mul; 00502 } // End for(int k... 00503 } 00504 i = j; 00505 } // End while(i+1<tup.n) 00506 00507 if( !canWrite ) myMatches.disableWriteAccess(); 00508 return MB_SUCCESS; 00509 } 00510 00511 // Sort the matching tuples so that vertices can be tagged accurately 00512 ErrorCode ParallelMergeMesh::SortMyMatches() 00513 { 00514 TupleList::buffer buf( mySkinEnts[0].size() ); 00515 // Sorts are necessary to check for doubles 00516 // Sort by remote handle 00517 myMatches.sort( 3, &buf ); 00518 // Sort by matching proc 00519 myMatches.sort( 1, &buf ); 00520 // Sort by local handle 00521 myMatches.sort( 2, &buf ); 00522 buf.reset(); 00523 return MB_SUCCESS; 00524 } 00525 00526 // Tag the shared elements using existing PComm functionality 00527 ErrorCode ParallelMergeMesh::TagSharedElements( int dim ) 00528 { 00529 // Manipulate the matches list to tag vertices and entities 00530 // Set up proc ents 00531 Range proc_ents; 00532 ErrorCode rval; 00533 00534 // get the entities in the partition sets 00535 for( Range::iterator rit = myPcomm->partitionSets.begin(); rit != myPcomm->partitionSets.end(); ++rit ) 00536 { 00537 Range tmp_ents; 00538 rval = myMB->get_entities_by_handle( *rit, tmp_ents, true ); 00539 if( MB_SUCCESS != rval ) 00540 { 00541 return rval; 00542 } 00543 proc_ents.merge( tmp_ents ); 00544 } 00545 if( myMB->dimension_from_handle( *proc_ents.rbegin() ) != myMB->dimension_from_handle( *proc_ents.begin() ) ) 00546 { 00547 Range::iterator lower = proc_ents.lower_bound( CN::TypeDimensionMap[0].first ), 00548 upper = proc_ents.upper_bound( CN::TypeDimensionMap[dim - 1].second ); 00549 proc_ents.erase( lower, upper ); 00550 } 00551 00552 // This vector doesn't appear to be used but its in resolve_shared_ents 00553 int maxp = -1; 00554 std::vector< int > sharing_procs( MAX_SHARING_PROCS ); 00555 std::fill( sharing_procs.begin(), sharing_procs.end(), maxp ); 00556 00557 // get ents shared by 1 or n procs 00558 std::map< std::vector< int >, std::vector< EntityHandle > > proc_nranges; 00559 Range proc_verts; 00560 rval = myMB->get_adjacencies( proc_ents, 0, false, proc_verts, Interface::UNION ); 00561 if( rval != MB_SUCCESS ) 00562 { 00563 return rval; 00564 } 00565 00566 rval = myPcomm->tag_shared_verts( myMatches, proc_nranges, proc_verts ); 00567 if( rval != MB_SUCCESS ) 00568 { 00569 return rval; 00570 } 00571 00572 // get entities shared by 1 or n procs 00573 rval = myPcomm->get_proc_nvecs( dim, dim - 1, &mySkinEnts[0], proc_nranges ); 00574 if( rval != MB_SUCCESS ) 00575 { 00576 return rval; 00577 } 00578 00579 // create the sets for each interface; store them as tags on 00580 // the interface instance 00581 Range iface_sets; 00582 rval = myPcomm->create_interface_sets( proc_nranges ); 00583 if( rval != MB_SUCCESS ) 00584 { 00585 return rval; 00586 } 00587 // establish comm procs and buffers for them 00588 std::set< unsigned int > procs; 00589 rval = myPcomm->get_interface_procs( procs, true ); 00590 if( rval != MB_SUCCESS ) 00591 { 00592 return rval; 00593 } 00594 00595 // resolve shared entity remote handles; implemented in ghost cell exchange 00596 // code because it's so similar 00597 rval = myPcomm->exchange_ghost_cells( -1, -1, 0, true, true ); 00598 if( rval != MB_SUCCESS ) 00599 { 00600 return rval; 00601 } 00602 // now build parent/child links for interface sets 00603 rval = myPcomm->create_iface_pc_links(); 00604 return rval; 00605 } 00606 00607 // Make sure to free up any allocated data 00608 // Need to avoid a double free 00609 void ParallelMergeMesh::CleanUp() 00610 { 00611 // The reset operation is now safe and avoids a double free() 00612 myMatches.reset(); 00613 myTup.reset(); 00614 myCD.reset(); 00615 } 00616 00617 // Simple quick sort to real 00618 void ParallelMergeMesh::SortTuplesByReal( TupleList& tup, double eps ) 00619 { 00620 bool canWrite = tup.get_writeEnabled(); 00621 if( !canWrite ) tup.enableWriteAccess(); 00622 00623 uint mi, ml, mul, mr; 00624 tup.getTupleSize( mi, ml, mul, mr ); 00625 PerformRealSort( tup, 0, tup.get_n(), eps, mr ); 00626 00627 if( !canWrite ) tup.disableWriteAccess(); 00628 } 00629 00630 // Swap around tuples 00631 void ParallelMergeMesh::SwapTuples( TupleList& tup, unsigned long a, unsigned long b ) 00632 { 00633 if( a == b ) return; 00634 00635 uint mi, ml, mul, mr; 00636 tup.getTupleSize( mi, ml, mul, mr ); 00637 00638 // Swap mi 00639 unsigned long a_val = a * mi, b_val = b * mi; 00640 for( unsigned long i = 0; i < mi; i++ ) 00641 { 00642 sint t = tup.vi_rd[a_val]; 00643 tup.vi_wr[a_val] = tup.vi_rd[b_val]; 00644 tup.vi_wr[b_val] = t; 00645 a_val++; 00646 b_val++; 00647 } 00648 // Swap ml 00649 a_val = a * ml; 00650 b_val = b * ml; 00651 for( unsigned long i = 0; i < ml; i++ ) 00652 { 00653 slong t = tup.vl_rd[a_val]; 00654 tup.vl_wr[a_val] = tup.vl_rd[b_val]; 00655 tup.vl_wr[b_val] = t; 00656 a_val++; 00657 b_val++; 00658 } 00659 // Swap mul 00660 a_val = a * mul; 00661 b_val = b * mul; 00662 for( unsigned long i = 0; i < mul; i++ ) 00663 { 00664 Ulong t = tup.vul_rd[a_val]; 00665 tup.vul_wr[a_val] = tup.vul_rd[b_val]; 00666 tup.vul_wr[b_val] = t; 00667 a_val++; 00668 b_val++; 00669 } 00670 // Swap mr 00671 a_val = a * mr; 00672 b_val = b * mr; 00673 for( unsigned long i = 0; i < mr; i++ ) 00674 { 00675 realType t = tup.vr_rd[a_val]; 00676 tup.vr_wr[a_val] = tup.vr_rd[b_val]; 00677 tup.vr_wr[b_val] = t; 00678 a_val++; 00679 b_val++; 00680 } 00681 } 00682 00683 // Perform the sorting of a tuple by real 00684 // To sort an entire tuple_list, call (tup,0,tup.n,epsilon) 00685 void ParallelMergeMesh::PerformRealSort( TupleList& tup, 00686 unsigned long left, 00687 unsigned long right, 00688 double eps, 00689 uint tup_mr ) 00690 { 00691 // If list size is only 1 or 0 return 00692 if( left + 1 >= right ) 00693 { 00694 return; 00695 } 00696 unsigned long swap = left, tup_l = left * tup_mr, tup_t = tup_l + tup_mr; 00697 00698 // Swap the median with the left position for a (hopefully) better split 00699 SwapTuples( tup, left, ( left + right ) / 2 ); 00700 00701 // Partition the data 00702 for( unsigned long t = left + 1; t < right; t++ ) 00703 { 00704 // If the left value(pivot) is greater than t_val, swap it into swap 00705 if( TupleGreaterThan( tup, tup_l, tup_t, eps, tup_mr ) ) 00706 { 00707 swap++; 00708 SwapTuples( tup, swap, t ); 00709 } 00710 tup_t += tup_mr; 00711 } 00712 00713 // Swap so that position swap is in the correct position 00714 SwapTuples( tup, left, swap ); 00715 00716 // Sort left and right of swap 00717 PerformRealSort( tup, left, swap, eps, tup_mr ); 00718 PerformRealSort( tup, swap + 1, right, eps, tup_mr ); 00719 } 00720 00721 // Note, this takes the actual tup.vr[] index (aka i*tup.mr) 00722 bool ParallelMergeMesh::TupleGreaterThan( TupleList& tup, 00723 unsigned long vrI, 00724 unsigned long vrJ, 00725 double eps, 00726 uint tup_mr ) 00727 { 00728 unsigned check = 0; 00729 while( check < tup_mr ) 00730 { 00731 // If the values are the same 00732 if( fabs( tup.vr_rd[vrI + check] - tup.vr_rd[vrJ + check] ) <= eps ) 00733 { 00734 check++; 00735 continue; 00736 } 00737 // If I greater than J 00738 else if( tup.vr_rd[vrI + check] > tup.vr_rd[vrJ + check] ) 00739 { 00740 return true; 00741 } 00742 // If J greater than I 00743 else 00744 { 00745 return false; 00746 } 00747 } 00748 // All Values are the same 00749 return false; 00750 } 00751 00752 } // End namespace moab