moab
|
#include <ParallelMergeMesh.hpp>
Public Member Functions | |
ParallelMergeMesh (ParallelComm *pc, const double epsilon) | |
ErrorCode | merge () |
Private Member Functions | |
ErrorCode | PerformMerge () |
ErrorCode | PopulateMySkinEnts (const EntityHandle meshset, int dim) |
ErrorCode | GetGlobalBox (double *gbox) |
ErrorCode | PopulateMyTup (double *gbox) |
ErrorCode | PopulateMyMatches () |
ErrorCode | SortMyMatches () |
ErrorCode | TagSharedElements (int dim) |
void | CleanUp () |
ErrorCode | PartitionGlobalBox (double *gbox, double *lengths, int *parts) |
Static Private Member Functions | |
static int | PartitionSide (double sideLeng, double restLen, unsigned numProcs, bool altRatio) |
static void | SwapTuples (TupleList &tup, unsigned long a, unsigned long b) |
static void | SortTuplesByReal (TupleList &tup, double eps2=0) |
static void | PerformRealSort (TupleList &tup, unsigned long left, unsigned long right, double eps2, uint tup_mr) |
static bool | TupleGreaterThan (TupleList &tup, unsigned long vrI, unsigned long vrJ, double eps2, uint tup_mr) |
Private Attributes | |
ParallelComm * | myPcomm |
Interface * | myMB |
std::vector< Range > | mySkinEnts |
double | myEps |
TupleList | myTup |
TupleList | myMatches |
gs_data::crystal_data | myCD |
Definition at line 23 of file ParallelMergeMesh.hpp.
moab::ParallelMergeMesh::ParallelMergeMesh | ( | ParallelComm * | pc, |
const double | epsilon | ||
) |
Definition at line 19 of file ParallelMergeMesh.cpp.
: myPcomm(pc), myEps(epsilon) { myMB = pc->get_moab(); mySkinEnts.resize(4); }
void moab::ParallelMergeMesh::CleanUp | ( | ) | [private] |
ErrorCode moab::ParallelMergeMesh::GetGlobalBox | ( | double * | gbox | ) | [private] |
Definition at line 149 of file ParallelMergeMesh.cpp.
{ ErrorCode rval; /*Get Bounding Box*/ BoundBox box; if(mySkinEnts[0].size() != 0){ rval = box.update(*myMB, mySkinEnts[0]); if(rval != MB_SUCCESS) return rval; } //Invert the max box.bMax *= -1; /*Communicate to all processors*/ MPI_Allreduce(&box, gbox, 6, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); /*Assemble Global Bounding Box*/ //Flip the max back for(int i=3; i<6; i++){ gbox[i] *= -1; } return MB_SUCCESS; }
Definition at line 30 of file ParallelMergeMesh.cpp.
{ ErrorCode rval = PerformMerge(); CleanUp(); return rval; }
ErrorCode moab::ParallelMergeMesh::PartitionGlobalBox | ( | double * | gbox, |
double * | lengths, | ||
int * | parts | ||
) | [private] |
Definition at line 284 of file ParallelMergeMesh.cpp.
{ //Determine the length of each side double xLen = gbox[3]-gbox[0]; double yLen = gbox[4]-gbox[1]; double zLen = gbox[5]-gbox[2]; unsigned numProcs = myPcomm->size(); //Partition sides from the longest to shortest lengths //If x is the longest side if(xLen >= yLen && xLen >= zLen){ parts[0] = PartitionSide(xLen, yLen * zLen, numProcs, true); numProcs /= parts[0]; //If y is second longest if(yLen >= zLen){ parts[1] = PartitionSide(yLen, zLen, numProcs, false); parts[2] = numProcs/parts[1]; } //If z is the longer else{ parts[2] = PartitionSide(zLen, yLen, numProcs, false); parts[1] = numProcs/parts[2]; } } //If y is the longest side else if (yLen >= zLen){ parts[1] = PartitionSide(yLen, xLen * zLen, numProcs, true); numProcs /= parts[1]; //If x is the second longest if(xLen >= zLen){ parts[0] = PartitionSide(xLen, zLen, numProcs, false); parts[2] = numProcs/parts[0]; } //If z is the second longest else{ parts[2] = PartitionSide(zLen, xLen, numProcs, false); parts[0] = numProcs/parts[2]; } } //If z is the longest side else{ parts[2] = PartitionSide(zLen, xLen * yLen, numProcs, true); numProcs /= parts[2]; //If x is the second longest if(xLen >= yLen){ parts[0] = PartitionSide(xLen, yLen, numProcs, false); parts[1] = numProcs/parts[0]; } //If y is the second longest else{ parts[1] = PartitionSide(yLen, xLen, numProcs, false); parts[0] = numProcs/parts[1]; } } //Divide up each side to give the lengths lengths[0] = xLen/(double)parts[0]; lengths[1] = yLen/(double)parts[1]; lengths[2] = zLen/(double)parts[2]; return MB_SUCCESS; }
int moab::ParallelMergeMesh::PartitionSide | ( | double | sideLeng, |
double | restLen, | ||
unsigned | numProcs, | ||
bool | altRatio | ||
) | [static, private] |
Definition at line 347 of file ParallelMergeMesh.cpp.
{ //If theres only 1 processor, then just return 1 if(numProcs == 1){ return 1; } //Initialize with the ratio of 1 proc double ratio = -DBL_MAX; unsigned factor = 1; //We need to be able to save the last ratio and factor (for comparison) double oldRatio = ratio; double oldFactor = 1; //This is the ratio were shooting for double goalRatio = sideLen/restLen; //Calculate the divisor and numerator power //This avoid if statements in the loop and is useful since both calculations are similar double divisor, p; if(altRatio){ divisor = (double)numProcs * sideLen; p = 3; } else{ divisor = (double)numProcs; p = 2; } //Find each possible factor for (unsigned i = 2; i <= numProcs/2; i++){ //If it is a factor... if (numProcs % i == 0){ //We need to save the past factor oldRatio = ratio; oldFactor = factor; //There are 2 different ways to calculate the ratio: //Comparing 1 side to 2 sides: (i*i*i)/(numProcs*x) //Justification: We have a ratio x:y:z (side Lengths) == a:b:c (procs). So a=kx, b=ky, c=kz. //Also, abc=n (numProcs) => bc = n/a. Also, a=kx => k=a/x => 1/k=x/a //And so x/(yz) == (kx)/(kyz) == (kx)/(kykz(1/k)) == a/(bc(x/a)) == a/((n/a)(x/a)) == a^3/(nx). //Comparing 1 side to 1 side: (i*i)/numprocs //Justification: i/(n/i) == i^2/n ratio = pow((double)i, p)/divisor; factor = i; //Once we have passed the goal ratio, we can break since we'll only move away from the goal ratio if(ratio >= goalRatio){ break; } } } //If we havent reached the goal ratio yet, check out factor = numProcs if(ratio < goalRatio){ oldRatio = ratio; oldFactor = factor; factor = numProcs; ratio = pow((double)numProcs, p)/divisor; } //Figure out if our oldRatio is better than ratio if(fabs(ratio - goalRatio) > fabs(oldRatio-goalRatio)){ factor = oldFactor; } //Return our findings return factor; }
ErrorCode moab::ParallelMergeMesh::PerformMerge | ( | ) | [private] |
Definition at line 38 of file ParallelMergeMesh.cpp.
{ //Get the mesh dimension int dim; ErrorCode rval = myMB->get_dimension(dim); if(rval != MB_SUCCESS){ return rval; } //Get the local skin elements rval = PopulateMySkinEnts(0,dim); //If there is only 1 proc, we can return now if(rval != MB_SUCCESS || myPcomm->size() == 1){ return rval; } //Determine the global bounding box double gbox[6]; rval = GetGlobalBox(gbox); if(rval != MB_SUCCESS){ return rval; } /* Assemble The Destination Tuples */ //Get a list of tuples which contain (toProc, handle, x,y,z) myTup.initialize(1,0,1,3,mySkinEnts[0].size()); rval = PopulateMyTup(gbox); if(rval != MB_SUCCESS){ return rval; } /* Gather-Scatter Tuple -tup comes out as (remoteProc,handle,x,y,z) */ myCD.initialize(myPcomm->comm()); //1 represents dynamic tuple, 0 represents index of the processor to send to myCD.gs_transfer(1, myTup, 0); /* Sort By X,Y,Z -Utilizes a custom quick sort incoroporating eplison*/ SortTuplesByReal(myTup,myEps); //Initilize another tuple list for matches myMatches.initialize(2,0,2,0,mySkinEnts[0].size()); //ID the matching tuples rval = PopulateMyMatches(); if(rval != MB_SUCCESS){ return rval; } //We can free up the tuple myTup now myTup.reset(); /*Gather-Scatter Again*/ //1 represents dynamic list, 0 represents proc index to send tuple to myCD.gs_transfer(1,myMatches,0); //We can free up the crystal router now myCD.reset(); //Sort the matches tuple list SortMyMatches(); //Tag the shared elements rval = TagSharedElements(dim); //Free up the matches tuples myMatches.reset(); return rval; }
void moab::ParallelMergeMesh::PerformRealSort | ( | TupleList & | tup, |
unsigned long | left, | ||
unsigned long | right, | ||
double | eps2, | ||
uint | tup_mr | ||
) | [static, private] |
Definition at line 653 of file ParallelMergeMesh.cpp.
{ //If list size is only 1 or 0 return if(left+1 >= right){ return; } unsigned long swap = left, tup_l = left*tup_mr, tup_t = tup_l + tup_mr; //Swap the median with the left position for a (hopefully) better split SwapTuples(tup, left, (left+right)/2); //Partition the data for(unsigned long t=left+1;t<right;t++){ //If the left value(pivot) is greater than t_val, swap it into swap if(TupleGreaterThan(tup,tup_l,tup_t,eps, tup_mr)){ swap++; SwapTuples(tup, swap,t); } tup_t+=tup_mr; } //Swap so that position swap is in the correct position SwapTuples(tup,left,swap); //Sort left and right of swap PerformRealSort(tup, left ,swap, eps, tup_mr); PerformRealSort(tup, swap+1,right,eps, tup_mr); }
ErrorCode moab::ParallelMergeMesh::PopulateMyMatches | ( | ) | [private] |
Definition at line 414 of file ParallelMergeMesh.cpp.
{ //Counters for accessing tuples more efficiently unsigned long i = 0, mat_i=0, mat_ul=0, j=0, tup_r=0; double eps2 = myEps*myEps; uint tup_mi, tup_ml, tup_mul, tup_mr; myTup.getTupleSize(tup_mi, tup_ml, tup_mul, tup_mr); bool canWrite = myMatches.get_writeEnabled(); if(!canWrite) myMatches.enableWriteAccess(); while((i+1)<myTup.get_n()){ //Proximity Comparison double xi = myTup.vr_rd[tup_r], yi = myTup.vr_rd[tup_r+1], zi = myTup.vr_rd[tup_r+2]; bool done = false; while(!done){ j++; tup_r += tup_mr; if(j >= myTup.get_n()){ break; } CartVect cv(myTup.vr_rd[tup_r]-xi, myTup.vr_rd[tup_r+1]-yi, myTup.vr_rd[tup_r+2]-zi); if(cv.length_squared() > eps2){ done = true; } } //Allocate the tuple list before adding matches while(myMatches.get_n()+(j-i)*(j-i-1) >= myMatches.get_max()){ myMatches.resize(myMatches.get_max() ? myMatches.get_max() + myMatches.get_max()/2 + 1 : 2); } //We now know that tuples [i to j) exclusive match. //If n tuples match, n*(n-1) match tuples will be made //tuples are of the form (proc1,proc2,handle1,handle2) if(i+1 < j){ int kproc = i*tup_mi; unsigned long khand = i*tup_mul; for(unsigned long k = i; k<j; k++){ int lproc = kproc+tup_mi; unsigned long lhand = khand+tup_mul; for(unsigned long l=k+1; l<j; l++){ myMatches.vi_wr[mat_i++]=myTup.vi_rd[kproc];//proc1 myMatches.vi_wr[mat_i++]=myTup.vi_rd[lproc];//proc2 myMatches.vul_wr[mat_ul++]=myTup.vul_rd[khand];//handle1 myMatches.vul_wr[mat_ul++]=myTup.vul_rd[lhand];//handle2 myMatches.inc_n(); myMatches.vi_wr[mat_i++]=myTup.vi_rd[lproc];//proc1 myMatches.vi_wr[mat_i++]=myTup.vi_rd[kproc];//proc2 myMatches.vul_wr[mat_ul++]=myTup.vul_rd[lhand];//handle1 myMatches.vul_wr[mat_ul++]=myTup.vul_rd[khand];//handle2 myMatches.inc_n(); lproc += tup_mi; lhand += tup_mul; } kproc += tup_mi; khand += tup_mul; }//End for(int k... } i = j; }//End while(i+1<tup.n) if(!canWrite) myMatches.disableWriteAccess(); return MB_SUCCESS; }
ErrorCode moab::ParallelMergeMesh::PopulateMySkinEnts | ( | const EntityHandle | meshset, |
int | dim | ||
) | [private] |
Definition at line 110 of file ParallelMergeMesh.cpp.
{ /*Merge Mesh Locally*/ //Get all dim dimensional entities Range ents; ErrorCode rval = myMB->get_entities_by_dimension(meshset,dim,ents); if(rval != MB_SUCCESS){ return rval; } //Merge Mesh Locally MergeMesh merger(myMB, false); merger.merge_entities(ents,myEps); //We can return if there is only 1 proc if(rval != MB_SUCCESS || myPcomm->size() == 1){ return rval; } //Rebuild the ents range ents.clear(); rval = myMB->get_entities_by_dimension(0,dim,ents); if(rval != MB_SUCCESS){ return rval; } /*Get Skin -Get Range of all dimensional entities -skinEnts[i] is the skin entities of dimension i*/ Skinner skinner(myMB); for(int skin_dim = dim; skin_dim >= 0; skin_dim--){ rval = skinner.find_skin(meshset,ents,skin_dim,mySkinEnts[skin_dim]); if(rval != MB_SUCCESS){ return rval; } } return MB_SUCCESS; }
ErrorCode moab::ParallelMergeMesh::PopulateMyTup | ( | double * | gbox | ) | [private] |
Definition at line 176 of file ParallelMergeMesh.cpp.
{ /*Figure out how do partition the global box*/ double lengths[3]; int parts[3]; ErrorCode rval = PartitionGlobalBox(gbox, lengths, parts); if(rval != MB_SUCCESS){ return rval; } /* Get Skin Coordinates, Vertices */ double *x = new double[mySkinEnts[0].size()]; double *y = new double[mySkinEnts[0].size()]; double *z = new double[mySkinEnts[0].size()]; rval = myMB->get_coords(mySkinEnts[0],x,y,z); if(rval != MB_SUCCESS){ //Prevent Memory Leak delete []x; delete []y; delete []z; return rval; } //Initialize variable to be used in the loops std::vector<int> toProcs; int xPart, yPart, zPart, xEps, yEps, zEps, baseProc; unsigned long long tup_i=0, tup_ul=0, tup_r=0, count=0; //These are boolean to determine if the vertice is on close enought to a given border bool xDup, yDup, zDup; bool canWrite = myTup.get_writeEnabled(); if(!canWrite) myTup.enableWriteAccess(); //Go through each vertice for(Range::iterator it = mySkinEnts[0].begin(); it != mySkinEnts[0].end(); it++){ xDup = false; yDup = false; zDup = false; //Figure out which x,y,z partition the element is in. xPart = static_cast<int>(floor((x[count]-gbox[0])/lengths[0])); xPart = (xPart<parts[0]?xPart:parts[0]-1);//Make sure it stays within the bounds yPart = static_cast<int>(floor((y[count]-gbox[1])/lengths[1])); yPart = (yPart<parts[1]?yPart:parts[1]-1);//Make sure it stays within the bounds zPart = static_cast<int>(floor((z[count]-gbox[2])/lengths[2])); zPart = (zPart<parts[2]?zPart:parts[2]-1);//Make sure it stays within the bounds //Figure out the partition with the addition of Epsilon xEps = static_cast<int>(floor((x[count]-gbox[0]+myEps)/lengths[0])); yEps = static_cast<int>(floor((y[count]-gbox[1]+myEps)/lengths[1])); zEps = static_cast<int>(floor((z[count]-gbox[2]+myEps)/lengths[2])); //Figure out if the vertice needs to be sent to multiple procs xDup = (xPart != xEps && xEps < parts[0]); yDup = (yPart != yEps && yEps < parts[1]); zDup = (zPart != zEps && zEps < parts[2]); //Add appropriate processors to the vector baseProc = xPart+ yPart * parts[0] + zPart * parts[0] * parts[1]; toProcs.push_back(baseProc); if(xDup){ toProcs.push_back(baseProc + 1);//Get partition to the right } if(yDup){ //Partition up 1 toProcs.push_back(baseProc + parts[0]); } if(zDup){ //Partition above 1 toProcs.push_back(baseProc + parts[0] * parts[1]); } if(xDup && yDup){ //Partition up 1 and right 1 toProcs.push_back(baseProc + parts[0] + 1); } if(xDup && zDup){ //Partition right 1 and above 1 toProcs.push_back(baseProc + parts[0] * parts[1] + 1); } if(yDup && zDup){ //Partition up 1 and above 1 toProcs.push_back(baseProc + parts[0] * parts[1] + parts[0]); } if(xDup && yDup && zDup){ //Partition right 1, up 1, and above 1 toProcs.push_back(baseProc + parts[0] * parts[1] + parts[0] + 1); } //Grow the tuple list if necessary while(myTup.get_n() + toProcs.size() >= myTup.get_max()){ myTup.resize(myTup.get_max() ? myTup.get_max() + myTup.get_max()/2 + 1 : 2); } //Add each proc as a tuple for(std::vector<int>::iterator proc = toProcs.begin(); proc != toProcs.end(); proc++){ myTup.vi_wr[tup_i++] = *proc; myTup.vul_wr[tup_ul++] = *it; myTup.vr_wr[tup_r++] = x[count]; myTup.vr_wr[tup_r++] = y[count]; myTup.vr_wr[tup_r++] = z[count]; myTup.inc_n(); } count++; toProcs.clear(); } if(!canWrite) myTup.disableWriteAccess(); return MB_SUCCESS; }
ErrorCode moab::ParallelMergeMesh::SortMyMatches | ( | ) | [private] |
Definition at line 488 of file ParallelMergeMesh.cpp.
{ TupleList::buffer buf(mySkinEnts[0].size()); //Sorts are necessary to check for doubles //Sort by remote handle myMatches.sort(3,&buf); //Sort by matching proc myMatches.sort(1,&buf); //Sort by local handle myMatches.sort(2,&buf); buf.reset(); return MB_SUCCESS; }
void moab::ParallelMergeMesh::SortTuplesByReal | ( | TupleList & | tup, |
double | eps2 = 0 |
||
) | [static, private] |
Definition at line 588 of file ParallelMergeMesh.cpp.
{ bool canWrite = tup.get_writeEnabled(); if(!canWrite) tup.enableWriteAccess(); uint mi, ml, mul, mr; tup.getTupleSize(mi,ml,mul,mr); PerformRealSort(tup, 0, tup.get_n(), eps, mr); if(!canWrite) tup.disableWriteAccess(); }
void moab::ParallelMergeMesh::SwapTuples | ( | TupleList & | tup, |
unsigned long | a, | ||
unsigned long | b | ||
) | [static, private] |
Definition at line 602 of file ParallelMergeMesh.cpp.
{ if(a==b) return; uint mi, ml, mul, mr; tup.getTupleSize(mi, ml, mul, mr); //Swap mi unsigned long a_val = a*mi, b_val=b*mi; for(unsigned long i=0; i< mi;i++){ sint t =tup.vi_rd[a_val]; tup.vi_wr[a_val] = tup.vi_rd[b_val]; tup.vi_wr[b_val] = t; a_val++; b_val++; } //Swap ml a_val = a*ml; b_val = b*ml; for(unsigned long i=0; i< ml;i++){ slong t =tup.vl_rd[a_val]; tup.vl_wr[a_val] = tup.vl_rd[b_val]; tup.vl_wr[b_val] = t; a_val++; b_val++; } //Swap mul a_val = a*mul; b_val = b*mul; for(unsigned long i=0; i< mul;i++){ ulong t =tup.vul_rd[a_val]; tup.vul_wr[a_val] = tup.vul_rd[b_val]; tup.vul_wr[b_val] = t; a_val++; b_val++; } //Swap mr a_val = a*mr; b_val = b*mr; for(unsigned long i=0; i< mr;i++){ realType t =tup.vr_rd[a_val]; tup.vr_wr[a_val] = tup.vr_rd[b_val]; tup.vr_wr[b_val] = t; a_val++; b_val++; } }
ErrorCode moab::ParallelMergeMesh::TagSharedElements | ( | int | dim | ) | [private] |
Definition at line 503 of file ParallelMergeMesh.cpp.
{ //Manipulate the matches list to tag vertices and entities //Set up proc ents Range proc_ents; ErrorCode rval; // get the entities in the partition sets for (Range::iterator rit = myPcomm->partitionSets.begin(); rit != myPcomm->partitionSets.end(); rit++) { Range tmp_ents; rval = myMB->get_entities_by_handle(*rit, tmp_ents, true); if (MB_SUCCESS != rval){ return rval; } proc_ents.merge(tmp_ents); } if (myMB->dimension_from_handle(*proc_ents.rbegin()) != myMB->dimension_from_handle(*proc_ents.begin())) { Range::iterator lower = proc_ents.lower_bound(CN::TypeDimensionMap[0].first), upper = proc_ents.upper_bound(CN::TypeDimensionMap[dim-1].second); proc_ents.erase(lower, upper); } //This vector doesnt appear to be used but its in resolve_shared_ents int maxp = -1; std::vector<int> sharing_procs(MAX_SHARING_PROCS); std::fill(sharing_procs.begin(), sharing_procs.end(), maxp); // get ents shared by 1 or n procs std::map<std::vector<int>, std::vector<EntityHandle> > proc_nranges; Range proc_verts; rval = myMB->get_adjacencies(proc_ents, 0, false, proc_verts, Interface::UNION); if(rval != MB_SUCCESS){ return rval; } rval = myPcomm->tag_shared_verts(myMatches, proc_nranges, proc_verts); if(rval != MB_SUCCESS){ return rval; } // get entities shared by 1 or n procs rval = myPcomm->get_proc_nvecs(dim,dim-1, &mySkinEnts[0],proc_nranges); if(rval != MB_SUCCESS){ return rval; } // create the sets for each interface; store them as tags on // the interface instance Range iface_sets; rval = myPcomm->create_interface_sets(proc_nranges); if(rval != MB_SUCCESS){ return rval; } // establish comm procs and buffers for them std::set<unsigned int> procs; rval = myPcomm->get_interface_procs(procs, true); if(rval != MB_SUCCESS){ return rval; } // resolve shared entity remote handles; implemented in ghost cell exchange // code because it's so similar rval = myPcomm->exchange_ghost_cells(-1, -1, 0, true, true); if(rval != MB_SUCCESS){ return rval; } // now build parent/child links for interface sets rval = myPcomm->create_iface_pc_links(); return rval; }
bool moab::ParallelMergeMesh::TupleGreaterThan | ( | TupleList & | tup, |
unsigned long | vrI, | ||
unsigned long | vrJ, | ||
double | eps2, | ||
uint | tup_mr | ||
) | [static, private] |
Definition at line 688 of file ParallelMergeMesh.cpp.
{ unsigned check=0; while(check < tup_mr){ //If the values are the same if(fabs(tup.vr_rd[vrI+check]-tup.vr_rd[vrJ+check]) <= eps){ check++; continue; } //If I greater than J else if(tup.vr_rd[vrI+check] > tup.vr_rd[vrJ+check]){ return true; } //If J greater than I else{ return false; } } //All Values are the same return false; }
Definition at line 37 of file ParallelMergeMesh.hpp.
double moab::ParallelMergeMesh::myEps [private] |
Definition at line 35 of file ParallelMergeMesh.hpp.
TupleList moab::ParallelMergeMesh::myMatches [private] |
Definition at line 36 of file ParallelMergeMesh.hpp.
Interface* moab::ParallelMergeMesh::myMB [private] |
Definition at line 33 of file ParallelMergeMesh.hpp.
ParallelComm* moab::ParallelMergeMesh::myPcomm [private] |
Definition at line 32 of file ParallelMergeMesh.hpp.
std::vector<Range> moab::ParallelMergeMesh::mySkinEnts [private] |
Definition at line 34 of file ParallelMergeMesh.hpp.
TupleList moab::ParallelMergeMesh::myTup [private] |
Definition at line 36 of file ParallelMergeMesh.hpp.