From c9cdcbf53971c2df32d27beaece523d1706d165a Mon Sep 17 00:00:00 2001
From: James McClure <mcclurej@vt.edu>
Date: Thu, 22 May 2014 14:43:49 -0400
Subject: [PATCH] Added two functions to assign blob IDs in parallel

---
 pmmc/AssignBlobID.cpp | 207 +++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 206 insertions(+), 1 deletion(-)

diff --git a/pmmc/AssignBlobID.cpp b/pmmc/AssignBlobID.cpp
index 1227aaff..a53b6658 100644
--- a/pmmc/AssignBlobID.cpp
+++ b/pmmc/AssignBlobID.cpp
@@ -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