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