MOAB: Mesh Oriented datABase  (version 5.4.0)
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 <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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines