From 2e47e9c306ef361b841c5b31bbfb78083324fc79 Mon Sep 17 00:00:00 2001
From: James McClure <mcclurej@vt.edu>
Date: Sun, 7 Jun 2020 20:54:09 -0400
Subject: [PATCH] fix bug in Euler characteristic at domain boundary

---
 analysis/Minkowski.cpp |  6 ++--
 models/ColorModel.cpp  | 74 +++++++++++++++++++++---------------------
 2 files changed, 40 insertions(+), 40 deletions(-)

diff --git a/analysis/Minkowski.cpp b/analysis/Minkowski.cpp
index faac6142..e98cfdc7 100644
--- a/analysis/Minkowski.cpp
+++ b/analysis/Minkowski.cpp
@@ -58,9 +58,9 @@ void Minkowski::ComputeScalar(const DoubleArray& Field, const double isovalue)
 	//int Nx = Field.size(0);
 	//int Ny = Field.size(1);
 	//int Nz = Field.size(2);
-	for (int k=1; k<Nz-1; k++){
-		for (int j=1; j<Ny-1; j++){
-			for (int i=1; i<Nx-1; i++){
+	for (int k=0; k<Nz-1; k++){
+		for (int j=0; j<Ny-1; j++){
+			for (int i=0; i<Nx-1; i++){
 				object.LocalIsosurface(Field,isovalue,i,j,k);
 				for (int idx=0; idx<object.TriangleCount; idx++){
 					e1 = object.Face(idx); 
diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 3d9d651d..2784e28a 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -1398,48 +1398,48 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 	}
 
 	if (USE_CONNECTED_NWP){
-	if (volume_connected - second_biggest < 2.0*fabs(target_delta_volume) && target_delta_volume < 0.0){
-		// if connected volume is less than 2% just delete the whole thing
-		if (rank==0) printf("Connected region has shrunk! \n");
-		REVERSE_FLOW_DIRECTION = true;
-	}
-	
-/*	else{*/
+		if (volume_connected - second_biggest < 2.0*fabs(target_delta_volume) && target_delta_volume < 0.0){
+			// if connected volume is less than 2% just delete the whole thing
+			if (rank==0) printf("Connected region has shrunk! \n");
+			REVERSE_FLOW_DIRECTION = true;
+		}
+
+		/*	else{*/
 		if (rank==0) printf("Pathway volume / next largest ganglion %f \n",volume_connected/second_biggest );
 	}
-		if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial);
-		double target_delta_volume_incremental = target_delta_volume;
-		if (fabs(target_delta_volume) > 0.01*volume_initial)  
-			target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume);
-		delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental, WallFactor);
+	if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial);
+	double target_delta_volume_incremental = target_delta_volume;
+	if (fabs(target_delta_volume) > 0.01*volume_initial)  
+		target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume);
+	delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental, WallFactor);
 
-		for (int k=0; k<Nz; k++){
-			for (int j=0; j<Ny; j++){
-				for (int i=0; i<Nx; i++){
-					if (phase_distance(i,j,k) < 0.0 ) phase_id(i,j,k) = 0;
-					else 		     				  phase_id(i,j,k) = 1;
-					//if (phase_distance(i,j,k) < 0.0 ) phase(i,j,k) = 1.0;
-				}
-			}
-		}	
-
-		CalcDist(phase_distance,phase_id,*Dm); // re-calculate distance
-
-		// 5. Update phase indicator field based on new distnace
-		for (int k=0; k<Nz; k++){
-			for (int j=0; j<Ny; j++){
-				for (int i=0; i<Nx; i++){
-					double d = phase_distance(i,j,k);
-					if (Averages->SDs(i,j,k) > 0.f){
-						if (d < 3.f){
-							//phase(i,j,k) = -1.0;
-							phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f);	
-						}
-					}
-				} 
+	for (int k=0; k<Nz; k++){
+		for (int j=0; j<Ny; j++){
+			for (int i=0; i<Nx; i++){
+				if (phase_distance(i,j,k) < 0.0 ) phase_id(i,j,k) = 0;
+				else 		     				  phase_id(i,j,k) = 1;
+				//if (phase_distance(i,j,k) < 0.0 ) phase(i,j,k) = 1.0;
 			}
 		}
-		fillDouble.fill(phase);
+	}	
+
+	CalcDist(phase_distance,phase_id,*Dm); // re-calculate distance
+
+	// 5. Update phase indicator field based on new distnace
+	for (int k=0; k<Nz; k++){
+		for (int j=0; j<Ny; j++){
+			for (int i=0; i<Nx; i++){
+				double d = phase_distance(i,j,k);
+				if (Averages->SDs(i,j,k) > 0.f){
+					if (d < 3.f){
+						//phase(i,j,k) = -1.0;
+						phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f);	
+					}
+				}
+			} 
+		}
+	}
+	fillDouble.fill(phase);
 	//}
 
 	count = 0.f;