From 4c26a3303822557c5ea8d920f880e876554c08f9 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Fri, 4 Sep 2015 16:13:32 -0400
Subject: [PATCH] Added other minkowski measures (curvature, surface area)

---
 common/TwoPhase.cpp            | 22 +++++-----------------
 common/pmmc.h                  | 18 ++++++++++++++++--
 tests/lbpm_color_simulator.cpp |  2 +-
 3 files changed, 22 insertions(+), 20 deletions(-)

diff --git a/common/TwoPhase.cpp b/common/TwoPhase.cpp
index 494336dc..c013b53c 100644
--- a/common/TwoPhase.cpp
+++ b/common/TwoPhase.cpp
@@ -367,14 +367,6 @@ void TwoPhase::ComputeLocal()
 	if (Dm.BoundaryCondition > 0 && Dm.kproc == 0) kmin=4;
 	if (Dm.BoundaryCondition > 0 && Dm.kproc == Dm.nprocz-1) kmax=Nz-4;
 
-	// Map solid to erode the fluid so that interfaces can be calculated accurately
-	for (k=0; k<Nz; k++){
-		for (j=0; j<Ny; j++){
-			for (i=0; i<Nx; i++){
-				SDs(i,j,k) += 1.0;
-			}
-		}
-	}
 	for (k=kmin; k<kmax; k++){
 		for (j=1; j<Ny-1; j++){
 			for (i=1; i<Nx-1; i++){
@@ -472,7 +464,11 @@ void TwoPhase::ComputeLocal()
 				//...........................................................................
 				// Compute the integral curvature of the non-wetting phase
 				n_nw_pts=n_nw_tris=0;
-				geomavg_MarchingCubes(SDn,fluid_isovalue,i,j,k,nw_pts,n_nw_pts,nw_tris,n_nw_tris);
+				// Compute the non-wetting phase surface and associated area
+				An += geomavg_MarchingCubes(SDn,fluid_isovalue,i,j,k,nw_pts,n_nw_pts,nw_tris,n_nw_tris);
+				// Compute the integral of mean curvature
+				Jn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,
+										i,j,k,n_nw_pts,n_nw_tris);
 				// Compute Euler characteristic from integral of gaussian curvature
 				euler += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,
 						i,j,k,n_nw_pts,n_nw_tris);
@@ -482,14 +478,6 @@ void TwoPhase::ComputeLocal()
 		}
 	}
 	
-	// Map solid back
-	for (k=0; k<Nz; k++){
-		for (j=0; j<Ny; j++){
-			for (i=0; i<Nx; i++){
-				SDs(i,j,k) -= 1.0;
-			}
-		}
-	}
 }
 
 
diff --git a/common/pmmc.h b/common/pmmc.h
index ff696d91..cfa16035 100644
--- a/common/pmmc.h
+++ b/common/pmmc.h
@@ -1,7 +1,6 @@
 #ifndef pmmc_INC
 #define pmmc_INC
 
-
 #include <stdio.h>
 #include <stdlib.h>
 #include <iostream>
@@ -1836,7 +1835,7 @@ inline void TRIM(DTMutableList<Point> &local_sol_pts, int &n_local_sol_pts, doub
 		}
 	}
 }
-inline void geomavg_MarchingCubes( DoubleArray &A, double &v, int &i, int &j, int &k,
+inline double geomavg_MarchingCubes( DoubleArray &A, double &v, int &i, int &j, int &k,
 				DTMutableList<Point> &nw_pts, int &n_nw_pts, IntArray &nw_tris,
 				int &n_nw_tris)
 {
@@ -1847,6 +1846,7 @@ inline void geomavg_MarchingCubes( DoubleArray &A, double &v, int &i, int &j, in
     int m;
     int o;
     int p;
+    int area;
 
 	// Go over each corner -- check to see if the corners are themselves vertices
 	//1
@@ -2149,6 +2149,20 @@ inline void geomavg_MarchingCubes( DoubleArray &A, double &v, int &i, int &j, in
     		n_nw_tris++;
     	}
     }
+	// Compute the Interfacial Area
+	double s1,s2,s3,s;
+	for (int r=n_nw_tris_beg;r<n_nw_tris;r++){
+		A = nw_pts(nw_tris(0,r));
+		B = nw_pts(nw_tris(1,r));
+		C = nw_pts(nw_tris(2,r));
+		// Compute length of sides (assume dx=dy=dz)
+		s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
+		s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
+		s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
+		s = 0.5*(s1+s2+s3);
+		area+=sqrt(s*(s-s1)*(s-s2)*(s-s3));
+	}
+    return area;
 }
 //-------------------------------------------------------------------------------
 inline void MC( DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, int &k,
diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index dff124d2..26dc1a9c 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -576,7 +576,7 @@ int main(int argc, char **argv)
 	MPI_Barrier(comm);
 	//.......................................................................
 	// Once phase has been initialized, map solid to account for 'smeared' interface
-	//for (i=0; i<N; i++)	Averages->SDs(i) -= (1.0); //
+	for (i=0; i<N; i++)	Averages->SDs(i) -= (1.0); //
 	//.......................................................................
 	// Finalize setup for averaging domain
 	//Averages->SetupCubes(Dm);