From 84690e28407945bdc8b76127ad5eb007dca822a7 Mon Sep 17 00:00:00 2001
From: James McClure <mcclurej@vt.edu>
Date: Sun, 8 Dec 2013 21:35:48 -0500
Subject: [PATCH] Contact angle support validated

---
 include/pmmc.h             | 17 +++++++
 tests/TestContactAngle.cpp | 94 +++++++++++---------------------------
 2 files changed, 43 insertions(+), 68 deletions(-)

diff --git a/include/pmmc.h b/include/pmmc.h
index 7e3cc973..f68fe3cb 100644
--- a/include/pmmc.h
+++ b/include/pmmc.h
@@ -3460,6 +3460,23 @@ inline double pmmc_CubeSurfaceArea(DTMutableList<Point> &Points, IntArray &Trian
 	}
 	return area;
 }
+inline double pmmc_CubeCurveLength(DTMutableList<Point> &Points, int npts)
+{
+	int p;
+	double s,lwns;
+	Point A,B;
+	lwns = 0.0;
+	for (p=0; p < npts-1; p++){
+		// Extract the line segment
+		A = Points(p);
+		B = Points(p+1);
+		// Compute the length of the segment
+		s = 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));
+		// Add the length to the common line 
+		lwns += s;
+	}
+	return lwns;
+}
 //--------------------------------------------------------------------------------------------------------
 inline double pmmc_CubeSurfaceInterpValue(DoubleArray &CubeValues, DTMutableList<Point> &Points, IntArray &Triangles,
 									  DoubleArray &SurfaceValues, int i, int j, int k, int npts, int ntris)
diff --git a/tests/TestContactAngle.cpp b/tests/TestContactAngle.cpp
index b04d8c42..3d78ffb0 100644
--- a/tests/TestContactAngle.cpp
+++ b/tests/TestContactAngle.cpp
@@ -5,6 +5,7 @@
 //#include "Array.h"
 
 #define RADIUS 15
+#define CAPRAD 20
 #define HEIGHT 15.5
 #define N 60
 #define PI 3.14159
@@ -24,6 +25,12 @@ int main (int argc, char *argv[])
 //	Phase = new double [SIZE];
 	DoubleArray SignDist(Nx,Ny,Nz);
 	DoubleArray Phase(Nx,Ny,Nz);
+	DoubleArray Fx(Nx,Ny,Nz);
+	DoubleArray Fy(Nx,Ny,Nz);
+	DoubleArray Fz(Nx,Ny,Nz);
+	DoubleArray Sx(Nx,Ny,Nz);
+	DoubleArray Sy(Nx,Ny,Nz);
+	DoubleArray Sz(Nx,Ny,Nz);
 	double fluid_isovalue = 0.0;
 	double solid_isovalue = 0.0;
 
@@ -35,6 +42,7 @@ int main (int argc, char *argv[])
 	// Averaging variables
 	//...........................................................................
 	double awn,ans,aws,lwns,nwp_volume;
+	double efawns;
 	double As;
 	double dEs,dAwn,dAns;			 // Global surface energy (calculated by rank=0)
 	double awn_global,ans_global,aws_global,lwns_global,nwp_volume_global;	
@@ -74,6 +82,9 @@ int main (int argc, char *argv[])
 	DTMutableList<Point> local_nws_pts(20);
 	int n_local_nws_pts;
 	
+	DoubleArray CubeValues(2,2,2);
+	DoubleArray ContactAngle(20);
+	
 	int c;
 	//...........................................................................
 	int ncubes = (Nx-2)*(Ny-2)*(Nz-2);	// Exclude the "upper" halo
@@ -86,17 +97,19 @@ int main (int argc, char *argv[])
 	for (k=0; k<N; k++){
 		for (j=0; j<N; j++){
 			for (i=0; i<N; i++){
+				
 				dist1 = sqrt((i-Cx)*(i-Cx)+(j-Cy)*(j-Cy)) - RADIUS;
-				dist2 = fabs(Cz-k)-HEIGHT;
-				//	printf("distances = %f, %f \n",dist1,dist2);
-				//Solid.data[k*Nx*Ny+j*Nx+i] = dist1;
-				//Phase[k*Nx*Ny+j*Nx+i] = dist2;
+				dist2 = sqrt((i-Cx)*(i-Cx)+(j-Cy)*(j-Cy)+(k-Cz)*(k-Cz)) - CAPRAD;
+
 				SignDist(i,j,k) = -dist1;
 				Phase(i,j,k) = dist2;
 			}
 		}   
 	}
 
+	pmmc_MeshGradient(Phase,Fx,Fy,Fz,Nx,Ny,Nz);
+	pmmc_MeshGradient(SignDist,Sx,Sy,Sz,Nx,Ny,Nz);
+	
 	FILE *STRIS;
 	STRIS = fopen("solid-triangles.out","w");
 	
@@ -146,72 +159,16 @@ int main (int argc, char *argv[])
 				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);
 
+		efawns += pmmc_CubeContactAngle(CubeValues,ContactAngle,Fx,Fy,Fz,Sx,Sy,Sz,local_nws_pts,i,j,k,n_local_nws_pts);
 		
 		//*******************************************************************
-		// Compute the Interfacial Areas, Common Line length for blob p
-		// nw surface
-		double temp;
-		for (r=0;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);
-			temp = s*(s-s1)*(s-s2)*(s-s3);
-			if (temp > 0.0) awn += sqrt(temp);
-			
-		}
-		for (r=0;r<n_ns_tris;r++){
-			A = ns_pts(ns_tris(0,r));
-			B = ns_pts(ns_tris(1,r));
-			C = ns_pts(ns_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);
-			//ans=ans+sqrt(s*(s-s1)*(s-s2)*(s-s3));
-			temp = s*(s-s1)*(s-s2)*(s-s3);
-			if (temp > 0.0) ans += sqrt(temp);
-		}
-		for (r=0;r<n_ws_tris;r++){
-			A = ws_pts(ws_tris(0,r));
-			B = ws_pts(ws_tris(1,r));
-			C = ws_pts(ws_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);
-			//aws=aws+sqrt(s*(s-s1)*(s-s2)*(s-s3));
-			temp = s*(s-s1)*(s-s2)*(s-s3);
-			if (temp > 0.0) aws += sqrt(temp);
-		}
-		for (r=0;r<n_local_sol_tris;r++){
-			A = local_sol_pts(local_sol_tris(0,r));
-			B = local_sol_pts(local_sol_tris(1,r));
-			C = local_sol_pts(local_sol_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);
-			//aws=aws+sqrt(s*(s-s1)*(s-s2)*(s-s3));
-			temp = s*(s-s1)*(s-s2)*(s-s3);
-			if (temp > 0.0) As += sqrt(temp);
-		}
-		for (p=0; p < n_local_nws_pts-1; p++){
-			// Extract the line segment
-			A = local_nws_pts(p);
-			B = local_nws_pts(p+1);
-			// Compute the length of the segment
-			s = 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));
-			// Add the length to the common line 
-			lwns += s;
-		}
+		// Compute the Interfacial Areas, Common Line length
+		awn += pmmc_CubeSurfaceArea(nw_pts,nw_tris,n_nw_tris);
+		ans += pmmc_CubeSurfaceArea(ns_pts,ns_tris,n_ns_tris);
+		aws += pmmc_CubeSurfaceArea(ws_pts,ws_tris,n_ws_tris);
+		As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris);
+		lwns +=  pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts);
+
 		//.......................................................................................
 		// Write the triangle lists to text file
 		for (r=0;r<n_nw_tris;r++){
@@ -267,6 +224,7 @@ int main (int argc, char *argv[])
 	printf("Area ws = %f, Analytical = %f \n", aws, 4*PI*RADIUS*HEIGHT);
 	printf("Area s = %f, Analytical = %f \n", As, 2*PI*RADIUS*(N-2));
 	printf("Length wns = %f, Analytical = %f \n", lwns, 4*PI*RADIUS);
+	printf("Cos(theta_wns) = %f, Analytical = %f \n", efawns/lwns,1.0*RADIUS/CAPRAD);
 	printf("-------------------------------- \n");	
 	//.........................................................................