From 83c0a8f0e744c984539b8aea06b4db91359dec52 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Sat, 11 Jul 2015 13:06:57 -0400
Subject: [PATCH] Working on component labeling

---
 tests/ComponentLabel.cpp | 429 +--------------------------------------
 1 file changed, 1 insertion(+), 428 deletions(-)

diff --git a/tests/ComponentLabel.cpp b/tests/ComponentLabel.cpp
index 07b3e9c9..b356a83b 100644
--- a/tests/ComponentLabel.cpp
+++ b/tests/ComponentLabel.cpp
@@ -493,7 +493,7 @@ int main(int argc, char **argv)
     int number_NWP_components = ComputeLocalPhaseComponent(PhaseLabel,1,NWP,false);
     int number_WP_components = ComputeLocalPhaseComponent(PhaseLabel,2,WP,false);
 
-	DoubleArray BlobAverages(NUM_AVERAGES,nblobs);
+	DoubleArray BlobAverages(NUM_AVERAGES,number_NWP_components);
 	
 	// Map the signed distance for the analysis
 	for (i=0; i<Nx*Ny*Nz; i++)	SignDist(i) -= (1.0); 
@@ -513,433 +513,6 @@ int main(int argc, char **argv)
 
 	printf("Media porosity is %f \n",porosity);
 
-	/* ****************************************************************
-			RUN TCAT AVERAGING ON EACH BLOB
-	****************************************************************** */
-	int n_nw_tris_beg, n_ns_tris_beg, n_ws_tris_beg, n_nws_seg_beg;
-	int start=0,finish;
-	int a,c;
-	
-	printf("-----------------------------------------------\n");
-	printf("Computing TCAT averages based on connectivity \n");
-	printf("The number of non-wetting phase features is %i \n",nblobs-1);
-	printf("-----------------------------------------------\n");
-
-	// Wetting phase averages assume global connectivity
-	As = 0.0;
-	vol_w = 0.0;
-	paw = 0.0;
-	aws = 0.0;
-	vaw(0) = vaw(1) = vaw(2) = 0.0;
-	Gws(0) = Gws(1) = Gws(2) = 0.0;
-	Gws(3) = Gws(4) = Gws(5) = 0.0;	
-
-	// Don't compute the last blob unless specified
-	// the last blob is the entire wetting phase
-//	nblobs -=1;
-#ifdef WP
-	nblobs+=1;
-#endif
-	
-	for (a=0;a<nblobs;a++){
-		
-		finish = start+b(a);
-
-		/* ****************************************************************
-		 RUN PMMC ON EACH BLOB
-		 ****************************************************************** */
-
-		// Store beginning points for surfaces for blob p
-		n_nw_tris_beg = n_nw_tris;
-		n_ns_tris_beg = n_ns_tris;
-		n_ws_tris_beg = n_ws_tris;
-		n_nws_seg_beg = n_nws_seg;
-		// Loop over all cubes
-		nwp_volume = 0.0;
-		
-		// Compute phase averages
-		vol_n =0.0;
-		pan = 0.0;
-		awn  = ans = lwns = 0.0;
-		van(0) = van(1) = van(2) = 0.0;
-		vawn(0) = vawn(1) = vawn(2) = 0.0;
-		Gwn(0) = Gwn(1) = Gwn(2) = 0.0;
-		Gwn(3) = Gwn(4) = Gwn(5) = 0.0;
-		Gns(0) = Gns(1) = Gns(2) = 0.0;
-		Gns(3) = Gns(4) = Gns(5) = 0.0;
-		Jwn = Kwn = efawns = 0.0;
-		trJwn = trawn = trRwn = 0.0;
-				
-		for (c=start;c<finish;c++){
-			// Get cube from the list
-			i = blobs(0,c);
-			j = blobs(1,c);
-			k = blobs(2,c);
-	
-			// Use the cube to compute volume averages
-			for (p=0;p<8;p++){
-				if ( SignDist(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
-					
-					n = i+cube[p][0] + Nx*(j+cube[p][1]) + Nx*Ny*(k+cube[p][2]);
-
-					// Compute the non-wetting phase volume contribution
-					if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 )
-						nwp_volume += 0.125;
-
-					// volume averages over the non-wetting phase
-					if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0.99 ){
-						// volume the excludes the interfacial region
-						vol_n += 0.125;
-						// pressure
-						pan += 0.125*Press(n);
-						// velocity
-						van(0) += 0.125*Vel_x(n);
-						van(1) += 0.125*Vel_y(n);
-						van(2) += 0.125*Vel_z(n);
-					}
-
-					// volume averages over the wetting phase
-					if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) < -0.99 ){
-						// volume the excludes the interfacial region
-						vol_w += 0.125;
-						// pressure
-						paw += 0.125*Press(n);
-						// velocity
-						vaw(0) += 0.125*Vel_x(n);
-						vaw(1) += 0.125*Vel_y(n);
-						vaw(2) += 0.125*Vel_z(n);
-					}
-				}
-			}
-
-			// Interface and common curve averages
-			n_local_sol_tris = 0;
-			n_local_sol_pts = 0;
-			n_local_nws_pts = 0;
-	
-			//...........................................................................
-			// Construct the interfaces and common curve
-			pmmc_ConstructLocalCube(SignDist, Phase, vS, vF,
-					nw_pts, nw_tris, values, ns_pts, ns_tris, ws_pts, ws_tris,
-					local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris,
-					n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
-					n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg,
-					i, j, k, Nx, Ny, Nz);
-
-			// Integrate the contact angle
-			efawns += pmmc_CubeContactAngle(CubeValues,Values,Phase_x,Phase_y,Phase_z,SignDist_x,SignDist_y,SignDist_z,
-											local_nws_pts,i,j,k,n_local_nws_pts);
-
-			// Integrate the mean curvature
-			Jwn    += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
-			Kwn    += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
-
-			// Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels)
-			pmmc_CubeTrimSurfaceInterpValues(CubeValues,MeanCurvature,SignDist,nw_pts,nw_tris,Values,DistValues,
-						i,j,k,n_nw_pts,n_nw_tris,trimdist,trawn,trJwn);
-			
-			pmmc_CubeTrimSurfaceInterpInverseValues(CubeValues,MeanCurvature,SignDist,nw_pts,nw_tris,Values,DistValues,
-						i,j,k,n_nw_pts,n_nw_tris,trimdist,dummy,trRwn);
-			
-			// Compute the normal speed of the interface
-			pmmc_InterfaceSpeed(dPdt, Phase_x, Phase_y, Phase_z, CubeValues, nw_pts, nw_tris,
-								NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
-			
-			As  += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris);
-
-			// Compute the surface orientation and the interfacial area
-			awn += pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris);
-			ans += pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris);
-			aws += pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris);
-			lwns +=  pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts);
-			//...........................................................................
-	
-			//*******************************************************************
-			// Reset the triangle counts to zero
-			n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0;
-			n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
-			
-			n_nw_tris_beg = n_nw_tris;
-			n_ns_tris_beg = n_ns_tris;
-			n_ws_tris_beg = n_ws_tris;
-			n_nws_seg_beg = n_nws_seg;
-			//*******************************************************************
-			
-		}
-		start = finish;
-
-		if (a < nblobs-1){
-		if (vol_n > 0.0){
-			pan /= vol_n;
-			for (i=0;i<3;i++)	van(i) /= vol_n;
-		}
-		if (awn > 0.0){
-			Jwn /= awn;
-			Kwn /= awn;
-			for (i=0;i<3;i++)	vawn(i) /= awn;
-			for (i=0;i<6;i++)	Gwn(i) /= awn;
-		}
-		if (lwns > 0.0){
-			efawns /= lwns;
-		}
-		if (ans > 0.0){
-			for (i=0;i<6;i++)	Gns(i) /= ans;
-		}
-		if (trawn > 0.0){
-			trJwn /= trawn;
-		}
-		
-		BlobAverages(0,a) = nwp_volume;
-		BlobAverages(1,a) = pan;
-		BlobAverages(2,a) = awn;
-		BlobAverages(3,a) = ans;
-		BlobAverages(4,a) = Jwn;
-		BlobAverages(5,a) = Kwn;
-		BlobAverages(6,a) = lwns;
-		BlobAverages(7,a) = efawns;
-		BlobAverages(8,a) = van(0);
-		BlobAverages(9,a) = van(1);
-		BlobAverages(10,a) = van(2);
-		BlobAverages(11,a) = vawn(0);
-		BlobAverages(12,a) = vawn(1);
-		BlobAverages(13,a) = vawn(2);
-		BlobAverages(14,a) = Gwn(0);
-		BlobAverages(15,a) = Gwn(1);
-		BlobAverages(16,a) = Gwn(2);
-		BlobAverages(17,a) = Gwn(3);
-		BlobAverages(18,a) = Gwn(4);
-		BlobAverages(19,a) = Gwn(5);
-		BlobAverages(20,a) = Gns(0);
-		BlobAverages(21,a) = Gns(1);
-		BlobAverages(22,a) = Gns(2);
-		BlobAverages(23,a) = Gns(3);
-		BlobAverages(24,a) = Gns(4);
-		BlobAverages(25,a) = Gns(5);
-		BlobAverages(26,a) = trawn;
-		BlobAverages(27,a) = trJwn;
-		BlobAverages(28,a) = vol_n;
-		BlobAverages(29,a) = trRwn;
-
-		printf("Computed TCAT averages for feature = %i \n", a);
-		}
-		
-	}  // End of the blob loop
-	NULL_USE(n_nw_tris_beg);
-    NULL_USE(n_ns_tris_beg);
-    NULL_USE(n_ws_tris_beg);
-    NULL_USE(n_nws_seg_beg);
-		
-	nblobs -= 1;
-	printf("-----------------------------------------------\n");
-	printf("Sorting the blobs based on volume \n");
-	printf("-----------------------------------------------\n");
-	int TempLabel,aa,bb;
-	double TempValue;
-	IntArray OldLabel(nblobs);
-	for (a=0; a<nblobs; a++)	OldLabel(a) = a;
-	// Sort the blob averages based on volume
-	for (aa=0; aa<nblobs-1; aa++){
-		for ( bb=aa+1; bb<nblobs; bb++){
-			if (BlobAverages(0,aa) < BlobAverages(0,bb)){
-				// Exchange location of blobs aa and bb
-				//printf("Switch blob %i with %i \n", OldLabel(aa),OldLabel(bb));
-				// switch the label
-				TempLabel = OldLabel(bb);
-				OldLabel(bb) = OldLabel(aa);
-				OldLabel(aa) = TempLabel;
-				// switch the averages
-				for (idx=0; idx<NUM_AVERAGES; idx++){
-					TempValue = BlobAverages(idx,bb);
-					BlobAverages(idx,bb) = BlobAverages(idx,aa);
-					BlobAverages(idx,aa) = TempValue;
-				}
-			}
-		}		
-	}
-	
-	IntArray NewLabel(nblobs);
-	for (aa=0; aa<nblobs; aa++){
-		// Match the new label for original blob aa
-		bb=0;
-		while (OldLabel(bb) != aa)	bb++;
-		NewLabel(aa) = bb;
-	}
-	
-	// Re-label the blob ID
-	printf("Re-labeling the blobs, now indexed by volume \n");
-	for (k=0; k<Nz; k++){
-		for (j=0; j<Ny; j++){
-			for (i=0; i<Nx; i++){
-				if (LocalBlobID(i,j,k) > -1){
-					TempLabel = NewLabel(LocalBlobID(i,j,k));
-					LocalBlobID(i,j,k) = TempLabel;
-				}
-			}
-		}
-	}
-		
-	FILE *BLOBLOG= fopen("blobs.tcat","a");
-	for (a=0; a<nblobs; a++){
-		//printf("Blob id =%i \n",a);
-		//printf("Original Blob id = %i \n",OldLabel(a));
-		//printf("Blob volume (voxels) = %f \n", BlobAverages(0,a));
-		for (idx=0; idx<28; idx++){
-			fprintf(BLOBLOG,"%.8g ",BlobAverages(idx,a));
-		}
-		fprintf(BLOBLOG,"\n");
-	}
-	fclose(BLOBLOG);
-	
-	double iVol = 1.0/Nx/Ny/Nz;
-	sw = 1.0;
-	// Compute the Sauter mean grain diamter
-	double D = 6.0*Nx*Ny*Nz*(1.0-porosity) / As;
-	double pw,pn,pc,awnD,ansD,awsD,JwnD,trJwnD,lwnsDD,cwns;
-	pw = paw/vol_w;
-	printf("paw = %f \n", paw/vol_w);
-	printf("vol_w = %f \n", vol_w);
-	
-	printf("-----------------------------------------------\n");
-	double pwn=0.0;
-	vol_n = nwp_volume = 0.0;
-	pan = 0.0;
-	awn  = ans = lwns = 0.0;
-	van(0) = van(1) = van(2) = 0.0;
-	vawn(0) = vawn(1) = vawn(2) = 0.0;
-	Gwn(0) = Gwn(1) = Gwn(2) = 0.0;
-	Gwn(3) = Gwn(4) = Gwn(5) = 0.0;
-	Gns(0) = Gns(1) = Gns(2) = 0.0;
-	Gns(3) = Gns(4) = Gns(5) = 0.0;
-	Jwn = Kwn = efawns = 0.0;
-	trJwn = trawn = trRwn = 0.0;	
-	
-	// Write out the "equilibrium" state with a 0.5 % change in saturation"
-	// Always write the largest blob 
-	// sw, pw, pn, pc*D/gamma, awn*D, aws*D, ans*D,lwns*D*D, Jwn*D, trJwn*D, cwns, nblobs
-	printf("Computing equilibria from blobs \n");
-	printf("Sauter mean diamter = %f \n",D);
-	printf("WARNING: lazy coder hard-coded the surface tension as 0.058 \n");
-	FILE *BLOBSTATES= fopen("blobstates.tcat","a");
-	fprintf(BLOBSTATES,"%.5g %.5g %.5g\n",vol_w,pw,aws);
-	// Compute the averages over the entire non-wetting phsae
-	for (a=0; a<nblobs; a++){
-		nwp_volume += BlobAverages(0,a);
-		pwn += (BlobAverages(1,a)-pw)*BlobAverages(2,a);
-		pan += BlobAverages(1,a)*BlobAverages(28,a);
-		awn += BlobAverages(2,a);
-		ans += BlobAverages(3,a);
-		Jwn += BlobAverages(4,a)*BlobAverages(2,a);
-		Kwn += BlobAverages(5,a)*BlobAverages(2,a);
-		lwns += BlobAverages(6,a);
-		efawns += BlobAverages(7,a)*BlobAverages(6,a);
-		van(0) += BlobAverages(8,a)*BlobAverages(28,a);
-		van(1) += BlobAverages(9,a)*BlobAverages(28,a);
-		van(2) += BlobAverages(10,a)*BlobAverages(28,a);
-		vawn(0) += BlobAverages(11,a)*BlobAverages(2,a);
-		vawn(1) += BlobAverages(12,a)*BlobAverages(2,a);
-		vawn(2) += BlobAverages(13,a)*BlobAverages(2,a);
-		Gwn(0) += BlobAverages(14,a)*BlobAverages(2,a);
-		Gwn(1) += BlobAverages(15,a)*BlobAverages(2,a);
-		Gwn(2) += BlobAverages(16,a)*BlobAverages(2,a);
-		Gwn(3) += BlobAverages(17,a)*BlobAverages(2,a);
-		Gwn(4) += BlobAverages(18,a)*BlobAverages(2,a);
-		Gwn(5) += BlobAverages(19,a)*BlobAverages(2,a);
-		Gns(0) += BlobAverages(20,a)*BlobAverages(3,a);
-		Gns(1) += BlobAverages(21,a)*BlobAverages(3,a);
-		Gns(2) += BlobAverages(22,a)*BlobAverages(3,a);
-		Gns(3) += BlobAverages(23,a)*BlobAverages(3,a);
-		Gns(4) += BlobAverages(24,a)*BlobAverages(3,a);
-		Gns(5) += BlobAverages(25,a)*BlobAverages(3,a);
-		trawn += BlobAverages(26,a);
-		trJwn += BlobAverages(27,a)*BlobAverages(26,a);
-		vol_n += BlobAverages(28,a);
-	}	
-	
-	// Subtract off portions of non-wetting phase in order of size
-	for (a=nblobs-1; a>0; a--){
-		// Subtract the features one-by-one
-		nwp_volume -= BlobAverages(0,a);
-		pan -= BlobAverages(1,a)*BlobAverages(28,a);
-		pwn -= (BlobAverages(1,a)-pw)*BlobAverages(2,a);
-		awn -= BlobAverages(2,a);
-		ans -= BlobAverages(3,a);
-		Jwn -= BlobAverages(4,a)*BlobAverages(2,a);
-		Kwn -= BlobAverages(5,a)*BlobAverages(2,a);
-		lwns -= BlobAverages(6,a);
-		efawns -= BlobAverages(7,a)*BlobAverages(6,a);
-		van(0) -= BlobAverages(8,a)*BlobAverages(28,a);
-		van(1) -= BlobAverages(9,a)*BlobAverages(28,a);
-		van(2) -= BlobAverages(10,a)*BlobAverages(28,a);
-		vawn(0) -= BlobAverages(11,a)*BlobAverages(2,a);
-		vawn(1) -= BlobAverages(12,a)*BlobAverages(2,a);
-		vawn(2) -= BlobAverages(13,a)*BlobAverages(2,a);
-		Gwn(0) -= BlobAverages(14,a)*BlobAverages(2,a);
-		Gwn(1) -= BlobAverages(15,a)*BlobAverages(2,a);
-		Gwn(2) -= BlobAverages(16,a)*BlobAverages(2,a);
-		Gwn(3) -= BlobAverages(17,a)*BlobAverages(2,a);
-		Gwn(4) -= BlobAverages(18,a)*BlobAverages(2,a);
-		Gwn(5) -= BlobAverages(19,a)*BlobAverages(2,a);
-		Gns(0) -= BlobAverages(20,a)*BlobAverages(3,a);
-		Gns(1) -= BlobAverages(21,a)*BlobAverages(3,a);
-		Gns(2) -= BlobAverages(22,a)*BlobAverages(3,a);
-		Gns(3) -= BlobAverages(23,a)*BlobAverages(3,a);
-		Gns(4) -= BlobAverages(24,a)*BlobAverages(3,a);
-		Gns(5) -= BlobAverages(25,a)*BlobAverages(3,a);
-		trawn -= BlobAverages(26,a);
-		trJwn -= BlobAverages(27,a)*BlobAverages(26,a);
-		vol_n -= BlobAverages(28,a);
-		
-		// Update wetting phase averages
-		aws += BlobAverages(3,a);
-		Gws(0) += BlobAverages(20,a)*BlobAverages(3,a);
-		Gws(1) += BlobAverages(21,a)*BlobAverages(3,a);
-		Gws(2) += BlobAverages(22,a)*BlobAverages(3,a);
-		Gws(3) += BlobAverages(23,a)*BlobAverages(3,a);
-		Gws(4) += BlobAverages(24,a)*BlobAverages(3,a);
-		Gws(5) += BlobAverages(25,a)*BlobAverages(3,a);
-		
-		if (fabs(1.0 - nwp_volume*iVol/porosity - sw) > 0.0025 || a == 1){
-			sw = 1.0 - nwp_volume*iVol/porosity;
-			
-			JwnD = -Jwn*D/awn;
-			pn = pan/vol_n;
-			trJwnD = -trJwn*D/trawn;
-			cwns = -efawns / lwns;
-			awnD = awn*D*iVol;
-			awsD = aws*D*iVol;
-			ansD = ans*D*iVol;
-			lwnsDD = lwns*D*D*iVol;
-			pc = pwn*D/0.058/awn; 	// hard-coded surface tension due to being lazy
-			
-			fprintf(BLOBSTATES,"%.5g %.5g %.5g ",sw,pn,pw);
-			fprintf(BLOBSTATES,"%.5g %.5g %.5g %.5g ",awnD,awsD,ansD,lwnsDD);
-			fprintf(BLOBSTATES,"%.5g %.5g %.5g %.5g %i\n",pc,JwnD,trJwnD,cwns,a);
-		}
-	}
-	fclose(BLOBSTATES);
-
-	start = 0;
-	for (a=0;a<nblobs;a++){
-
-		finish = start+b(a);
-
-		for (c=start;c<finish;c++){
-			// Get cube from the list
-			i = blobs(0,c);
-			j = blobs(1,c);
-			k = blobs(2,c);
-			// Label the entire cube so that interfaces can be re-labled easily
-			LocalBlobID(i,j,k) = NewLabel(a);
-			LocalBlobID(i+1,j,k) = NewLabel(a);
-			LocalBlobID(i,j+1,k) = NewLabel(a);
-			LocalBlobID(i+1,j+1,k) = NewLabel(a);
-			LocalBlobID(i,j,k+1) = NewLabel(a);
-			LocalBlobID(i+1,j,k+1) = NewLabel(a);
-			LocalBlobID(i,j+1,k+1) = NewLabel(a);
-			LocalBlobID(i+1,j+1,k+1) = NewLabel(a);
-		}
-		start=finish;
-	}
 
 	FILE *BLOBS;
 	BLOBS = fopen("Blobs.dat","wb");