Added two functions to assign blob IDs in parallel
This commit is contained in:
James McClure
parent
b5cbdfa213
commit
c9cdcbf539
@ -19,6 +19,212 @@ inline double PhaseAverageScalar(int *PhaseID, double *Scalar, int N, int Np)
|
||||
|
||||
}
|
||||
}
|
||||
inline int AssignGlobalBlobID(IntArray &blobs, int &nblobs, int &ncubes, IntArray &indicator, &IntArray cubes_in_blob, int number_of_blobs_below_rank)
|
||||
{
|
||||
// This routine takes the local blob ID and assigns a globally unique tag to each blob
|
||||
// the local tag is already unique, we just need to add it to the number of blobs determined
|
||||
// for all ranks < the rank of the processor calling this routine
|
||||
// Note that the total number of connected blobs is less than number_of_blobs_below_rank
|
||||
// since the tag that results after calling this routine does not reflect connectivity
|
||||
// between different MPI sub-domains.
|
||||
int i,j,k;
|
||||
int start=0,finish;
|
||||
int a,c;
|
||||
|
||||
int minimumBlobID;
|
||||
|
||||
for (a=0;a<nblobs;a++){
|
||||
finish = start+cubes_in_blob(a);
|
||||
//...............
|
||||
for (c=start;c<finish;c++){
|
||||
// Get cube from the list
|
||||
i = blobs(0,c);
|
||||
j = blobs(1,c);
|
||||
k = blobs(2,c);
|
||||
// reassign a unique ID for the whole blob
|
||||
indicator(i,j,k) += number_of_blobs_below_rank;
|
||||
}
|
||||
start = finish;
|
||||
}
|
||||
}
|
||||
|
||||
inline int ReassignLocalBlobID(IntArray &blobs, int &nblobs, int &ncubes, IntArray &indicator, &IntArray cubes_in_blob)
|
||||
{
|
||||
// This routine reassigns all blobs within given MPI process so that a unique tag is assigned to
|
||||
// all cubes within the blob. If a blob is entirely within the sub-domain, the blob will already
|
||||
// be uniquely tagged. If the blob crosses one or more processor boundaries, multiple tags will be
|
||||
// assigned and this routine will assign the minimum tag found within the blob boundaries.
|
||||
// The routine returns the number of blobs within the MPI sub-domain that have had ID reassignments.
|
||||
// To explore the final connectivity, this routine must be run iteratively until no MPI
|
||||
// process reassigns the label for a cube.
|
||||
int i,j,k;
|
||||
int start=0,finish;
|
||||
int a,c;
|
||||
int count_reassigned_blobs=0;
|
||||
|
||||
int minimumBlobID;
|
||||
|
||||
for (a=0;a<nblobs;a++){
|
||||
finish = start+cubes_in_blob(a);
|
||||
//...............
|
||||
i = blobs(0,start);
|
||||
j = blobs(1,start);
|
||||
k = blobs(2,start);
|
||||
//...............
|
||||
minimumBlobID = indicator(i,j,k);
|
||||
//...............
|
||||
// There may be multiple blob IDs for one blob due to parallel assignment
|
||||
// Find the minimum blob ID within this blob
|
||||
for (c=start;c<finish;c++){
|
||||
// Get cube from the list
|
||||
i = blobs(0,c);
|
||||
j = blobs(1,c);
|
||||
k = blobs(2,c);
|
||||
if (IntArray(i,j,k) < minimumBlobID){
|
||||
minimumBlobID = indicator(i,j,k);
|
||||
count_reassigned_blobs++;
|
||||
}
|
||||
}
|
||||
for (c=start;c<finish;c++){
|
||||
// Get cube from the list
|
||||
i = blobs(0,c);
|
||||
j = blobs(1,c);
|
||||
k = blobs(2,c);
|
||||
// reassign a unique ID for the whole blob
|
||||
indicator(i,j,k) = minimumBlobID;
|
||||
}
|
||||
start = finish;
|
||||
}
|
||||
return count_reassigned_blobs;
|
||||
}
|
||||
|
||||
inline int ComputeLocalBlob(IntArray blobs, int &nblobs, int &ncubes, IntArray indicator,
|
||||
DoubleArray F, DoubleArray S, double vf, double vs, int startx, int starty,
|
||||
int startz, IntArray temp)
|
||||
{
|
||||
// Compute the blob (F>vf|S>vs) starting from (i,j,k) - oil blob
|
||||
// F>vf => oil phase S>vs => in porespace
|
||||
// update the list of blobs, indicator mesh
|
||||
int m = F.m; // maxima for the meshes
|
||||
int n = F.n;
|
||||
int o = F.o;
|
||||
|
||||
int cubes_in_blob=0;
|
||||
int nrecent = 1; // number of nodes added at most recent sweep
|
||||
temp(0,0) = startx; // Set the initial point as a "seed" for the sweeps
|
||||
temp(1,0) = starty;
|
||||
temp(2,0) = startz;
|
||||
int ntotal = 1; // total number of nodes in blob
|
||||
indicator(startx,starty,startz) = nblobs;
|
||||
|
||||
int p,s,x,y,z,start,finish,nodx,nody,nodz;
|
||||
int imin=startx,imax=startx,jmin=starty,jmax=starty; // initialize maxima / minima
|
||||
int kmin=startz,kmax=startz;
|
||||
int d[26][3] = {{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1},
|
||||
{1,1,0},{1,-1,0},{-1,1,0},{-1,-1,0},{1,0,1},{-1,0,1},
|
||||
{1,0,-1},{-1,0,-1},{0,1,1},{0,-1,1},{0,1,-1},{0,-1,-1},
|
||||
{1,1,1},{1,1,-1},{1,-1,1},{1,-1,-1},{-1,1,1},{-1,1,-1},
|
||||
{-1,-1,1},{-1,-1,-1}}; // directions to neighbors
|
||||
int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}}; // cube corners
|
||||
bool status = 1; // status == true => continue to look for points
|
||||
while (status == 1){
|
||||
start = ntotal - nrecent;
|
||||
finish = ntotal;
|
||||
nrecent = 0; // set recent points back to zero for next sweep through
|
||||
for (s=start;s<finish;s++){
|
||||
// Loop over recent points; look for new points
|
||||
x = temp(0,s);
|
||||
y = temp(1,s);
|
||||
z = temp(2,s);
|
||||
// Looop over the directions
|
||||
for (p=0;p<26;p++){
|
||||
nodx=x+d[p][0];
|
||||
if (nodx < 0 ){ nodx = m-1; } // Periodic BC for x
|
||||
if (nodx > m-1 ){ nodx = 0; }
|
||||
nody=y+d[p][1];
|
||||
if (nody < 0 ){ nody = n-1; } // Periodic BC for y
|
||||
if (nody > n-1 ){ nody = 0; }
|
||||
nodz=z+d[p][2];
|
||||
if (nodz < 0 ){ nodz = 0; } // No periodic BC for z
|
||||
if (nodz > o-1 ){ nodz = o-1; }
|
||||
if ( F(nodx,nody,nodz) > vf && S(nodx,nody,nodz) > vs
|
||||
&& indicator(nodx,nody,nodz) == -1 ){
|
||||
// Node is a part of the blob - add it to the list
|
||||
temp(0,ntotal) = nodx;
|
||||
temp(1,ntotal) = nody;
|
||||
temp(2,ntotal) = nodz;
|
||||
ntotal++;
|
||||
nrecent++;
|
||||
// Update the indicator map
|
||||
indicator(nodx,nody,nodz) = nblobs;
|
||||
// Update the min / max for the cube loop
|
||||
if ( nodx < imin ){ imin = nodx; }
|
||||
if ( nodx > imax ){ imax = nodx; }
|
||||
if ( nody < jmin ){ jmin = nody; }
|
||||
if ( nody > jmax ){ jmax = nody; }
|
||||
if ( nodz < kmin ){ kmin = nodz; }
|
||||
if ( nodz > kmax ){ kmax = nodz; }
|
||||
}
|
||||
else if (F(nodx,nody,nodz) > vf && S(nodx,nody,nodz) > vs
|
||||
&& indicator(nodx,nody,nodz) > -1 && indicator(nodx,nody,nodz) != nblobs){
|
||||
// Some kind of error in algorithm
|
||||
printf("Error in blob search algorithm!");
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
if ( nrecent == 0){
|
||||
status = 0;
|
||||
}
|
||||
}
|
||||
// Use points in temporary storage array to add cubes to the list of blobs
|
||||
if ( imin > 0) { imin = imin-1; }
|
||||
// if ( imax < m-1) { imax = imax+1; }
|
||||
if ( jmin > 0) { jmin = jmin-1; }
|
||||
// if ( jmax < n-1) { jmax = jmax+1; }
|
||||
if ( kmin > 0) { kmin = kmin-1; }
|
||||
// if ( kmax < o-1) { kmax = kmax+1; }
|
||||
int i,j,k;
|
||||
bool add;
|
||||
for (k=kmin;k<kmax;k++){
|
||||
for (j=jmin;j<jmax;j++){
|
||||
for (i=imin;i<imax;i++){
|
||||
// If cube(i,j,k) has any nodes in blob, add it to the list
|
||||
// Loop over cube edges
|
||||
add = 0;
|
||||
for (p=0;p<8;p++){
|
||||
nodx = i+cube[p][0];
|
||||
nody = j+cube[p][1];
|
||||
nodz = k+cube[p][2];
|
||||
if ( indicator(nodx,nody,nodz) == nblobs ){
|
||||
// Cube corner is in this blob
|
||||
add = 1;
|
||||
}
|
||||
}
|
||||
if (add == 1){
|
||||
// Add cube to the list
|
||||
blobs(0,ncubes) = i;
|
||||
blobs(1,ncubes) = j;
|
||||
blobs(2,ncubes) = k;
|
||||
ncubes++;
|
||||
cubes_in_blob++;
|
||||
// Loop again to check for overlap
|
||||
for (p=0;p<8;p++){
|
||||
nodx = i+cube[p][0];
|
||||
nody = j+cube[p][1];
|
||||
nodz = k+cube[p][2];
|
||||
if (indicator(nodx,nody,nodz) > -1 && indicator(nodx,nody,nodz) != nblobs){
|
||||
printf("Overlapping cube!");
|
||||
cout << i << ", " << j << ", " << k << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return cubes_in_blob;
|
||||
}
|
||||
|
||||
inline int ComputeBlob(IntArray blobs, int &nblobs, int &ncubes, IntArray indicator,
|
||||
DoubleArray F, DoubleArray S, double vf, double vs, int startx, int starty,
|
||||
@ -31,7 +237,6 @@ inline int ComputeBlob(IntArray blobs, int &nblobs, int &ncubes, IntArray indica
|
||||
int n = F.n;
|
||||
int o = F.o;
|
||||
|
||||
|
||||
int cubes_in_blob=0;
|
||||
int nrecent = 1; // number of nodes added at most recent sweep
|
||||
temp(0,0) = startx; // Set the initial point as a "seed" for the sweeps
|
||||
|
||||
Loading…
Reference in New Issue
Block a user