MOAB: Mesh Oriented datABase  (version 5.2.1)
ParallelMergeMesh.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines