Merge branch 'master' of https://github.com/JamesEMcClure/LBPM-WIA

common/Domain.h

@@ -34,6 +34,7 @@ struct Domain{
 		BlobLabel.resize(Nx,Ny,Nz);
 		BlobGraph.resize(18,MAX_BLOB_COUNT,MAX_BLOB_COUNT);
 		BoundaryCondition = BC;
+		rank_info=RankInfoStruct(rank,nprocx,nprocy,nprocz);
 	}
 	~Domain();
 
@@ -44,6 +45,7 @@ struct Domain{
     double Lx,Ly,Lz,Volume;
 	int rank;
 	int BoundaryCondition;
+	RankInfoStruct rank_info;
 	MPI_Group Group;	// Group of processors associated with this domain
 	MPI_Comm Comm;		// MPI Communicator for this domain
 
@@ -775,6 +777,142 @@ void Domain::BlobComm(MPI_Comm Communicator){
 	UnpackBlobData(recvList_YZ, recvCount_YZ ,recvBuf_YZ, BlobLabelData);
 	//......................................................................................
 }
+inline void SSO(DoubleArray &Distance, char *ID, Domain &Dm, int timesteps){
+	/*
+	 * This routine converts the data in the Distance array to a signed distance
+	 * by solving the equation df/dt = sign(1-|grad f|), where Distance provides
+	 * the values of f on the mesh associated with domain Dm
+	 * It has been tested with segmented data initialized to values [-1,1]
+	 * and will converge toward the signed distance to the surface bounding the associated phases
+	 */
+
+	int Q=26;
+    int q,i,j,k,n;
+    double dt=0.1;
+    int in,jn,kn,nn;
+    double Dqx,Dqy,Dqz,Dx,Dy,Dz,W;
+    double nx,ny,nz,Cqx,Cqy,Cqz,sign,norm;
+
+    const static int D3Q27[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}};
+
+    double weights[26];
+    // Compute the weights from the finite differences
+    for (q=0; q<Q; q++){
+        weights[q] = sqrt(1.0*(D3Q27[q][0]*D3Q27[q][0]) + 1.0*(D3Q27[q][1]*D3Q27[q][1]) + 1.0*(D3Q27[q][2]*D3Q27[q][2]));
+    }
+
+    int xdim,ydim,zdim;
+    xdim=Dm.Nx-2;
+    ydim=Dm.Ny-2;
+    zdim=Dm.Nz-2;
+    fillHalo<double> fillData(Dm.rank_info,xdim,ydim,zdim,1,1,1,0,1);
+
+    int count = 0;
+    while (count < timesteps){
+
+        // Communicate the halo of values
+        fillData.fill(Distance);
+
+        // Execute the next timestep
+        for (k=1;k<Dm.Nz-1;k++){
+            for (j=1;j<Dm.Ny-1;j++){
+                for (i=1;i<Dm.Nx-1;i++){
+
+                	n = k*Dm.Nx*Dm.Ny+j*Dm.Nx+i;
+                    sign = Distance(i,j,k) / fabs(Distance(i,j,k));
+
+              /*
+                    if (!(i+1<Nx)) 	nx=0.5*Distance(i,j,k);
+                    else 			nx=0.5*Distance(i+1,j,k);;
+                    if (!(j+1<Ny)) 	ny=0.5*Distance(i,j,k);
+                    else 			ny=0.5*Distance(i,j+1,k);
+                    if (!(k+1<Nz)) 	nz=0.5*Distance(i,j,k);
+                    else 			nz=0.5*Distance(i,j,k+1);
+                    if (i<1)  		nx-=0.5*Distance(i,j,k);
+                    else 			nx-=0.5*Distance(i-1,j,k);
+                    if (j<1) 		ny-=0.5*Distance(i,j,k);
+                    else 			ny-=0.5*Distance(i,j-1,k);
+                    if (k<1)  		nz-=0.5*Distance(i,j,k);
+                    else 			nz-=0.5*Distance(i,j,k-1);
+                    */
+
+                    //............Compute the Gradient...................................
+                    nx = 0.5*(Distance(i+1,j,k) - Distance(i-1,j,k));
+                    ny = 0.5*(Distance(i,j+1,k) - Distance(i,j-1,k));
+                    nz = 0.5*(Distance(i,j,k+1) - Distance(i,j,k-1));
+
+                    W = 0.0;	Dx = Dy = Dz = 0.0;
+                    // Ignore any values that have distances less than a lattice unit
+                    // since sometimes the positions may be guessed more accurately from
+                    // another source (such as a simulation)
+                    // also ignore places where the gradient is zero since this will not
+                    // result in any local change to Distance
+                    if (nx*nx+ny*ny+nz*nz > 0.0 && !(Distance(i,j,k)*Distance(i,j,k) < 1.0) ){
+                        for (q=0; q<26; q++){
+                            Cqx = 1.0*D3Q27[q][0];
+                            Cqy = 1.0*D3Q27[q][1];
+                            Cqz = 1.0*D3Q27[q][2];
+                            // get the associated neighbor
+                            in = i + D3Q27[q][0];
+                            jn = j + D3Q27[q][1];
+                            kn = k + D3Q27[q][2];
+
+                            // make sure the neighbor is in the domain (periodic BC)
+                            /*				if (in < 0 ) in +=Nx;
+                             * 	don't need this in parallel since MPI handles the halos
+                             if (jn < 0 ) jn +=Ny;
+                             if (kn < 0 ) kn +=Nz;
+                             if (!(in < Nx) ) in -=Nx;
+                             if (!(jn < Ny) ) jn -=Ny;
+                             if (!(kn < Nz) ) kn -=Nz;
+                             			// symmetric boundary
+                            if (in < 0 ) in = i;
+                            if (jn < 0 ) jn = j;
+                            if (kn < 0 ) kn = k;
+                            if (!(in < Nx) ) in = i;
+                            if (!(jn < Ny) ) jn = k;
+                            if (!(kn < Nz) ) kn = k;
+                            */
+
+                            // Compute the gradient using upwind finite differences
+                            Dqx = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqx;
+                            Dqy = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqy;
+                            Dqz = weights[q]*(Distance(i,j,k) - Distance(in,jn,kn))*Cqz;
+
+                            // Only include upwind derivatives
+                            if (sign*(nx*Cqx + ny*Cqy + nz*Cqz) < 0.0 ){
+
+                                Dx += Dqx;
+                                Dy += Dqy;
+                                Dz += Dqz;
+                                W += weights[q];
+                            }
+                        }
+                        // Normalize by the weight to get the approximation to the gradient
+                        Dx /= W;
+                        Dy /= W;
+                        Dz /= W;
+
+                        norm = sqrt(Dx*Dx+Dy*Dy+Dz*Dz);
+                    }
+                    else{
+                        norm = 0.0;
+                    }
+                    Distance(i,j,k) += dt*sign*(1.0 - norm);
+
+                    // Disallow any change in phase
+                    if (Distance(i,j,k)*2.0*(ID[n]-1.0) < 0) Distance(i,j,k) = -Distance(i,j,k);
+                }
+            }
+        }
+
+        count++;
+    }
+}
+
 
 inline void ReadSpherePacking(int nspheres, double *List_cx, double *List_cy, double *List_cz, double *List_rad)
 {

common/TwoPhase.h

@@ -550,7 +550,7 @@ void TwoPhase::ComputeLocalBlob(){
     int i,j,k,n,label;
 	double vF,vS;
 	vF = 0.0; vS= -1.0;
-    const RankInfoStruct rank_info(Dm.rank,Dm.nprocx,Dm.nprocy,Dm.nprocz);
+//    const RankInfoStruct rank_info(Dm.rank,Dm.nprocx,Dm.nprocy,Dm.nprocz);
 
 	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}};
         // get the maximum label locally -- then compute number of global blobs
@@ -563,7 +563,7 @@ void TwoPhase::ComputeLocalBlob(){
 	nblobs_global+=1;
 	if (Dm.rank==0) printf("Number of blobs is %i \n",nblobs_global);
 
-    nblobs_global = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,
+    nblobs_global = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,Dm.rank_info,
     		Phase,SDs,vF,vS,BlobLabel);
 
 	if (Dm.rank==0) printf("Number of blobs is %i \n",nblobs_global);

tests/BlobAnalyzeParallel.cpp

@@ -7,7 +7,9 @@
 #include <math.h>
 #include "common/Communication.h"
 #include "analysis/analysis.h"
+#ifdef PROFILE
 	#include "ProfilerApp.h"
+#endif
 #include "TwoPhase.h"
 
 //#include "Domain.h"
@@ -98,10 +100,12 @@ int main(int argc, char **argv)
 	MPI_Init(&argc,&argv);
 	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
 	MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
+#ifdef PROFILE
 	PROFILE_ENABLE(0);
     PROFILE_DISABLE_TRACE();
     PROFILE_SYNCHRONIZE();
     PROFILE_START("main");
+#endif
 
     if ( rank==0 ) {
         printf("-----------------------------------------------------------\n");
@@ -147,7 +151,7 @@ int main(int argc, char **argv)
 	int BC=0;
     // Get the rank info
 	Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
-    const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
+ //   const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
 	TwoPhase Averages(Dm);
 	int N = (nx+2)*(ny+2)*(nz+2);
 
@@ -157,6 +161,7 @@ int main(int argc, char **argv)
     readRankData( rank, nx+2, ny+2, nz+2, Phase, SignDist );
 
     // Communication the halos
+    const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
     fillHalo<double> fillData(rank_info,nx,ny,nz,1,1,1,0,1);
     fillData.fill(Phase);
     fillData.fill(SignDist);
@@ -166,7 +171,7 @@ int main(int argc, char **argv)
     double vF=0.0;
     double vS=0.0;
     IntArray GlobalBlobID;
-    int nblobs = ComputeGlobalBlobIDs(nx,ny,nz,rank_info,
+    int nblobs = ComputeGlobalBlobIDs(nx,ny,nz,Dm.rank_info,
         Phase,SignDist,vF,vS,GlobalBlobID);
     if ( rank==0 ) { printf("Identified %i blobs\n",nblobs); }
 
@@ -191,9 +196,6 @@ int main(int argc, char **argv)
 	//.........................................................................
 	// Populate the arrays needed to perform averaging
 	if (rank==0) printf("Populate arrays \n");
-	// Compute porosity
-	double porosity,sum,sum_global;
-	sum=0.0;
 	for (int k=0; k<nz+2; k++){
 		for (int j=0; j<ny+2; j++){
 			for (int i=0; i<nx+2; i++){
@@ -237,6 +239,20 @@ int main(int argc, char **argv)
 				Averages.Vel_x(i,j,k)=vx;
 				Averages.Vel_y(i,j,k)=vy;
 				Averages.Vel_z(i,j,k)=vz;
+
+			}
+		}
+	}
+    delete [] DistEven;
+    delete [] DistOdd;
+
+	// Compute porosity
+	double porosity,sum,sum_global;
+	sum=0.0;
+    for (int k=1; k<nz+1; k++){
+        for (int j=1; j<ny+1; j++){
+            for (int i=1; i<nx+1; i++){
+				int n = k*(nx+2)*(ny+2)+j*(nx+2)+i;
 				if (Averages.SDs(i,j,k) > 0.0){
 					Dm.id[n]=1;
 					sum += 1.0;
@@ -245,9 +261,6 @@ int main(int argc, char **argv)
             }
         }
     }
-    delete [] DistEven;
-    delete [] DistOdd;
-
     MPI_Allreduce(&sum,&sum_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
     porosity = sum_global/Dm.Volume;
     if (rank==0) printf("Porosity = %f \n",porosity);
@@ -366,9 +379,10 @@ int main(int argc, char **argv)
     /*FILE *BLOBS = fopen("Blobs.dat","wb");
     fwrite(GlobalBlobID.get(),4,Nx*Ny*Nz,BLOBS);
     fclose(BLOBS);*/
-
+#ifdef PROFILE
     PROFILE_STOP("main");
     PROFILE_SAVE("BlobIdentifyParallel",false);
+#endif
     MPI_Barrier(MPI_COMM_WORLD);
     MPI_Finalize();
     return 0;  

tests/BlobIdentifyParallel.cpp

@@ -8,8 +8,9 @@
 #include "common/pmmc.h"
 #include "common/Communication.h"
 #include "analysis/analysis.h"
+#ifdef PROFILE
  #include "ProfilerApp.h"
-
+#endif
 
 //#include "Domain.h"
 
@@ -50,10 +51,12 @@ int main(int argc, char **argv)
 	MPI_Init(&argc,&argv);
 	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
 	MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
+#ifdef PROFILE
 	PROFILE_ENABLE(0);
     PROFILE_DISABLE_TRACE();
     PROFILE_SYNCHRONIZE();
     PROFILE_START("main");
+#endif
 
     if ( rank==0 ) {
         printf("-----------------------------------------------------------\n");
@@ -121,9 +124,10 @@ int main(int argc, char **argv)
     /*FILE *BLOBS = fopen("Blobs.dat","wb");
     fwrite(GlobalBlobID.get(),4,Nx*Ny*Nz,BLOBS);
     fclose(BLOBS);*/
-
+#ifdef PROFILE
     PROFILE_STOP("main");
     PROFILE_SAVE("BlobIdentifyParallel",false);
+#endif
     MPI_Barrier(MPI_COMM_WORLD);
 	MPI_Finalize();
     return 0;  

tests/CMakeLists.txt

@@ -22,6 +22,7 @@ ADD_LBPM_TEST( TestSphereCurvature )
 ADD_LBPM_TEST( TestContactAngle )
 ADD_LBPM_TEST_1_2_4( TestTwoPhase )
 ADD_LBPM_TEST_PARALLEL( TestTwoPhase 8 )
+ADD_LBPM_TEST_PARALLEL( TestSegDist 8 )
 ADD_LBPM_TEST_1_2_4( testCommunication )
 ADD_LBPM_TEST_1_2_4( testUtilities )
 ADD_LBPM_TEST( TestWriter )

tests/TestComputeSignDist.cpp

@@ -11,13 +11,13 @@
 #include <fstream>
 #include <Array.h>
 
-inline void SSO(DoubleArray &Distance, IntArray &ID, int timesteps){
+inline void SSO(DoubleArray &Distance, char *ID, int timesteps){
 
 	int Q=26;
-	int q,i,j,k;
+	int q,i,j,k,n;
-	int Nx = ID.m;
+	int Nx = Distance.m;
-	int Ny = ID.n;
+	int Ny = Distance.n;
-	int Nz = ID.o;
+	int Nz = Distance.o;
 	const static int D3Q27[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},
@@ -30,40 +30,58 @@ inline void SSO(DoubleArray &Distance, IntArray &ID, int timesteps){
 	}
 
 	// Initialize the Distance from ID
-	for (i=0; i<Nx*Ny*Nz; i++) Distance.data[i] = -0.5;
+//	for (i=0; i<Nx*Ny*Nz; i++) Distance.data[i] = -0.5;
-	for (k=1;k<Nz-1;k++){
+	for (k=0;k<Nz;k++){
-		for (j=1;j<Ny-1;j++){
+		for (j=0;j<Ny;j++){
-			for (i=1;i<Nx-1;i++){
+			for (i=0;i<Nx;i++){
+				n=k*Nx*Ny+j*Nx+i;
 				// Initialize distance to +/- 1
-				Distance(i,j,k) = 1.0*ID(i,j,k)-0.5;
+				Distance(i,j,k) = 1.0*ID[n]-0.5;
 			}
 		}
 	}
 
 	int count = 0;
-	double dt=1.0;
+	double dt=0.1;
-	int n,in,jn,kn,nn;
+	int in,jn,kn,nn;
 	double Dqx,Dqy,Dqz,Dx,Dy,Dz,W;
 	double nx,ny,nz,Cqx,Cqy,Cqz,sign,norm;
 	double f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
 
+	printf("Number of timesteps is %i \n",timesteps);
+	printf("Mesh is %i,%i,%i \n",Nx,Ny,Nz);
+
 	while (count < timesteps){
 
-		for (k=1;k<Nz-1;k++){
+		printf("count=%i \n",count);
-			for (j=1;j<Ny-1;j++){
+		for (k=0;k<Nz;k++){
-				for (i=1;i<Nx-1;i++){
+			for (j=0;j<Ny;j++){
+				for (i=0;i<Nx;i++){
 
 					n = k*Nx*Ny + j*Nx + i;
 					sign = Distance.data[n] / fabs(Distance.data[n]);
 
 					//............Compute the Gradient...................................
-					nx = 0.5*(Distance(i+1,j,k) - Distance(i-1,j,k));
+					if (!(i+1<Nx)) 	nx=0.5*Distance(i,j,k);
-					ny = 0.5*(Distance(i,j+1,k) - Distance(i,j-1,k));
+					else 			nx=0.5*Distance(i+1,j,k);;
-					nz = 0.5*(Distance(i,j,k+1) - Distance(i,j,k-1));
+					if (!(j+1<Ny)) 	ny=0.5*Distance(i,j,k);
+					else 			ny=0.5*Distance(i,j+1,k);
+					if (!(k+1<Nz)) 	nz=0.5*Distance(i,j,k);
+					else 			nz=0.5*Distance(i,j,k+1);
+					if (i<1)  		nx-=0.5*Distance(i,j,k);
+					else 			nx-=0.5*Distance(i-1,j,k);
+					if (j<1) 		ny-=0.5*Distance(i,j,k);
+					else 			ny-=0.5*Distance(i,j-1,k);
+					if (k<1)  		nz-=0.5*Distance(i,j,k);
+					else 			nz-=0.5*Distance(i,j,k-1);
+
+//					nx = 0.5*(Distance(i+1,j,k) - Distance(i-1,j,k));
+	//				ny = 0.5*(Distance(i,j+1,k) - Distance(i,j-1,k));
+		//			nz = 0.5*(Distance(i,j,k+1) - Distance(i,j,k-1));
 
 					W = 0.0;	Dx = Dy = Dz = 0.0;
 					if (nx*nx+ny*ny+nz*nz > 0.0){
-						for (q=0; q<27; q++){
+						for (q=0; q<26; q++){
 							Cqx = 1.0*D3Q27[q][0];
 							Cqy = 1.0*D3Q27[q][1];
 							Cqz = 1.0*D3Q27[q][2];
@@ -72,12 +90,19 @@ inline void SSO(DoubleArray &Distance, IntArray &ID, int timesteps){
 							jn = j + D3Q27[q][1];
 							kn = k + D3Q27[q][2];
 							// make sure the neighbor is in the domain (periodic BC)
-							if (in < 0 ) in +=Nx;
+			/*				if (in < 0 ) in +=Nx;
 							if (jn < 0 ) jn +=Ny;
 							if (kn < 0 ) kn +=Nz;
 							if (!(in < Nx) ) in -=Nx;
 							if (!(jn < Ny) ) jn -=Ny;
 							if (!(kn < Nz) ) kn -=Nz;
+			*/				// symmetric boundary
+							if (in < 0 ) in = i;
+							if (jn < 0 ) jn = j;
+							if (kn < 0 ) kn = k;
+							if (!(in < Nx) ) in = i;
+							if (!(jn < Ny) ) jn = k;
+							if (!(kn < Nz) ) kn = k;
 							// 1-D index
 							nn = kn*Nx*Ny + jn*Nx + in;
 
@@ -103,20 +128,18 @@ inline void SSO(DoubleArray &Distance, IntArray &ID, int timesteps){
 						norm = sqrt(Dx*Dx+Dy*Dy+Dz*Dz);
 					}
 					else{
-						norm = 0.7;
+						norm = 0.0;
 					}
 					Distance.data[n] += dt*sign*(1.0 - norm);
 
-					// Disallow any change in phase position
+					// Disallow any change in phase
-					if (Distance.data[n]*ID.data[n] < 0) Distance.data[n] = -Distance.data[n];
+					if (Distance.data[n]*2.0*(ID[n]-1.0) < 0) Distance.data[n] = -Distance.data[n];
 				}
 			}
 		}
 
 		count++;
 	}
-
-
 }
 
 int main(){
@@ -137,27 +160,30 @@ int main(){
 	double fxy,fXy,fxY,fXY,fxz,fXz,fxZ,fXZ,fyz,fYz,fyZ,fYZ;
 	double nx,ny,nz;
 	
-	
 	int ip,im,jp,jm,kp,km;
 	
 	int count = 0;
 	double sign;
-	double dt=0.01;
+	double dt=1.0;
 	printf("Nx=%i, Ny=%i, Nz= %i, \n",Nx,Ny,Nz);
 	
-	short int *id;
+	char *id;
 
 #ifdef READMEDIA
 	Nx = 347;
 	Ny = 347;
 	Nz = 235;
+	Nx = 512;
+	Ny = 512;
+	Nz = 512;
 	N = Nx*Ny*Nz;
-	id = new short int [N];
+	id = new char [N];
-	FILE *INPUT = fopen("Solid.dat",'rb');
+
-	fread(id,4,N*,Ny*Nz,INPUT)
+	FILE *INPUT = fopen("Solid.dat","rb");
+	fread(id,1,Nx*Ny*Nz,INPUT);
 	fclose(INPUT);
 #else
-	id = new short int [N];
+	id = new char [N];
 	double BubbleRadius = 5;
 	// Initialize the bubble
 	for (k=0;k<Nz;k++){
@@ -176,6 +202,21 @@ int main(){
 	}
 #endif
 	
+	DoubleArray Distance(Nx,Ny,Nz);
+	// Initialize the signed distance function
+	for (k=0;k<Nz;k++){
+		for (j=0;j<Ny;j++){
+			for (i=0;i<Nx;i++){
+				n=k*Nx*Ny+j*Nx+i;
+				// Initialize distance to +/- 1
+				Distance(i,j,k) = 2.0*ID[n]-1.0;
+			}
+		}
+	}
+
+	printf("Initialized! Converting to Signed Distance function \n");
+	SSO(Distance,id,20);
+
 /*	double *f,*f_old,*f_new;
 	f = new double[N];
 	f_old = new double[N];
@@ -184,9 +225,11 @@ int main(){
 	for (int n=0; n<N; n++)	Distance.data[n] = 0.5*(id[n]-1);
 	for (int n=0; n<N; n++)	f_old[n] = Distance.data[n];
 	for (int n=0; n<N; n++)	f_new[n] = Distance.data[n];
+	for (int n=0; n<N; n++)	f[n] = Distance.data[n];
 
 	count=0;
-	while (count < 10000 && dt > 1.0e-6){
+	dt=1.0;
+	while (count < 10 && dt > 1.0e-6){
 
 		err = 0.0;
 		for (k=0;k<Nz;k++){
@@ -306,7 +349,7 @@ int main(){
 					f_y = 0.5*(f[k*Nx*Ny + jp*Nx + i] - f[k*Nx*Ny + jm*Nx + i]);
 					f_z = 0.5*(f[kp*Nx*Ny + j*Nx + i] - f[km*Nx*Ny + j*Nx + i]);					
 
-					if (id[n] > 0){
+					if (id[n] < 0){
 						if ( nx > 0.0)	f_x = fxp;
 						else			f_x = fxm;
 						if ( ny > 0.0)	f_y = fyp;
@@ -348,8 +391,8 @@ int main(){
 		
 		count++;			
 	}
-
 */
+
 	FILE *DIST;
 	DIST = fopen("SignDist","wb");
 	fwrite(Distance.data,8,N,DIST);

tests/TestSegDist.cpp

@@ -0,0 +1,106 @@
+// Compute the signed distance from a digitized image 
+// Two phases are present
+// Phase 1 has value -1
+// Phase 2 has value 1
+// this code uses the segmented image to generate the signed distance 
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <iostream>
+#include <fstream>
+#include <Array.h>
+#include <Domain.h>
+
+int main(int argc, char **argv)
+{
+	// Initialize MPI
+	int rank, nprocs;
+	MPI_Init(&argc,&argv);
+	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+	MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
+	
+	int i,j,k,n,nn;
+	int iproc,jproc,kproc;
+	int nx,ny,nz;
+	int Nx, Ny, Nz, N;
+    int nprocx, nprocy, nprocz, nspheres;
+    double Lx, Ly, Lz;
+	Nx = Ny = Nz = 50;
+	nx = ny = nz = 50;
+	N = Nx*Ny*Nz;
+	nprocx=nprocy=nprocz=2;
+	Lx = Ly = Lz = 1.0;
+	int BC=0;
+
+	
+	if (nprocs != 8){
+		ERROR("TestSegDist: Number of MPI processes must be equal to 8");
+	}
+
+    if (nprocx !=2 || nprocz !=2 || nprocy !=2 ){
+		ERROR("TestSegDist: MPI process grid must be 2x2x2");
+	}
+
+    // Get the rank info
+	Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
+	Dm.CommInit(MPI_COMM_WORLD);
+
+	nx+=2; ny+=2; nz+=2;
+	N = nx*ny*nz;
+	int count = 0;
+
+	char *id;
+	id = new char [N];
+	double BubbleRadius = 25;
+	// Initialize the bubble
+	int x,y,z;
+	for (k=1;k<nz-1;k++){
+		for (j=1;j<ny-1;j++){
+			for (i=1;i<nx-1;i++){
+				x = (nx-2)*Dm.iproc+i;
+				y = (ny-2)*Dm.jproc+j;
+				z = (nz-2)*Dm.kproc+k;
+				n = k*nx*ny+j*nx+i;
+
+				// Initialize phase positions
+				if ((x-nx+1)*(x-nx+1)+(y-ny+1)*(y-ny+1)+(z-nz+1)*(z-nz+1) < BubbleRadius*BubbleRadius){
+					id[n] = 0;
+				}
+				else{
+					id[n]=1;
+				}
+			}
+		}
+	}
+	
+	DoubleArray Distance(nx,ny,nz);
+	// Initialize the signed distance function
+	for (k=0;k<nz;k++){
+		for (j=0;j<ny;j++){
+			for (i=0;i<nx;i++){
+				n=k*nx*ny+j*nx+i;
+				// Initialize distance to +/- 1
+				Distance(i,j,k) = 2.0*id[n]-1.0;
+			}
+		}
+	}
+
+	if (rank==0) printf("Nx = %i \n",(int)Distance.size(0));
+	if (rank==0) printf("Ny = %i \n",(int)Distance.size(1));
+	if (rank==0) printf("Nz = %i \n",(int)Distance.size(2));
+
+	printf("Initialized! Converting to Signed Distance function \n");
+	SSO(Distance,id,Dm,10);
+
+	char LocalRankFilename[40];
+    sprintf(LocalRankFilename,"Dist.%05i",rank);
+    FILE *DIST = fopen(LocalRankFilename,"wb");
+    fwrite(Distance.get(),8,Distance.length(),DIST);
+    fclose(DIST);
+
+    MPI_Barrier(MPI_COMM_WORLD);
+	MPI_Finalize();
+    return 0;
+
+}

tests/TestTwoPhase.cpp

@@ -106,7 +106,6 @@ int main(int argc, char **argv)
 	Averages.PrintAll(timestep);
 	//....................................................................
 
-
 	if (rank==0){
 		FILE *PHASE;
 		PHASE = fopen("Phase.00000","wb");