diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index 79ef5c93..f17c7b83 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -20,6 +20,7 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
 	Pressure.resize(Nx,Ny,Nz);      Pressure.fill(0);
 	Phi.resize(Nx,Ny,Nz);         	Phi.fill(0);
 	DelPhi.resize(Nx,Ny,Nz);        DelPhi.fill(0);
+	Laplacian.resize(Nx,Ny,Nz);     Laplacian.fill(0);
 	Vel_x.resize(Nx,Ny,Nz);         Vel_x.fill(0);	    // Gradient of the phase indicator field
 	Vel_y.resize(Nx,Ny,Nz);         Vel_y.fill(0);
 	Vel_z.resize(Nx,Ny,Nz);         Vel_z.fill(0);
@@ -172,8 +173,21 @@ void SubPhase::Basic(){
 	double count_w = 0.0;
 	double count_n = 0.0;
 	
-	total_wetting_interaction = count_wetting_interaction = 0.0;
-	total_wetting_interaction_global = count_wetting_interaction_global=0.0;
+	/* compute the laplacian */
+ 	Dm->CommunicateMeshHalo(Phi);
+	for (int k=1; k<Nz-1; k++){
+		for (int j=1; j<Ny-1; j++){
+			for (int i=1; i<Nx-1; i++){
+				// Compute all of the derivatives using finite differences
+				double fx = (Phi(i+1,j,k) - 2.0*Phi(i,j,k) + Phi(i-1,j,k));
+				double fy = (Phi(i,j+1,k) - 2.0*Phi(i,j,k) + Phi(i,j-1,k));
+				double fz = (Phi(i,j,k+1) - 2.0*Phi(i,j,k) + Phi(i,j,k-1));
+				Laplacian(i,j,k) = sqrt(fx*fx+fy*fy+fz*fz);
+			}
+		}
+	}
+	Dm->Comm.barrier();
+ 	Dm->CommunicateMeshHalo(Laplacian);
 	
 	for (k=0; k<Nz; k++){
 		for (j=0; j<Ny; j++){
@@ -234,132 +248,20 @@ void SubPhase::Basic(){
 						count_w += 1.0;
 					}
 				}
-				else {
-					// solid wetting assessment
-					double wetval = Phi(i,j,k);
-					//if (wetval != wetval) printf("%f at %i %i %i \n",wetval,i,j,k);
-					double local_wetting_interaction = 0.0;
-					double local_wetting_weight=0.0;
-					nq = (k)*Nx*Ny+(j)*Nx+(i+1);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 1\n",i,j,k);
-						 local_wetting_interaction += (Phi(nq)-wetval);
-						 local_wetting_weight += 1.0;
-					}
-					nq = (k)*Nx*Ny+(j)*Nx+(i-1);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 2\n",i,j,k);
-						 local_wetting_interaction += (Phi(nq)-wetval);
-						 local_wetting_weight += 1.0;
-					}
-					nq = (k)*Nx*Ny+(j+1)*Nx+(i);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 3\n",i,j,k);
-						 local_wetting_interaction += (Phi(nq)-wetval);
-						 local_wetting_weight += 1.0;
-					}
-					nq = (k)*Nx*Ny+(j-1)*Nx+(i);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 4\n",i,j,k);
-						 local_wetting_interaction += (Phi(nq)-wetval);
-						 local_wetting_weight += 1.0;
-					}
-					nq = (k+1)*Nx*Ny+(j)*Nx+(i);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 5\n",i,j,k);
-						 local_wetting_interaction += (Phi(nq)-wetval);
-						 local_wetting_weight += 1.0;
-					}
-					nq = (k-1)*Nx*Ny+(j)*Nx+(i);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 6\n",i,j,k);
-						 local_wetting_interaction += (Phi(nq)-wetval);
-						 local_wetting_weight += 1.0;
-					}
-					// x, y interactions
-					nq = (k)*Nx*Ny+(j+1)*Nx+(i+1);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 7\n",i,j,k);
-						 local_wetting_interaction += 0.5*(Phi(nq)-wetval);
-						 local_wetting_weight += 0.5;
-					}
-					nq = (k)*Nx*Ny+(j-1)*Nx+(i-1);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 8\n",i,j,k);
-						 local_wetting_interaction += 0.5*(Phi(nq)-wetval);
-						 local_wetting_weight += 0.5;
-					}
-					nq = (k)*Nx*Ny+(j-1)*Nx+(i+1);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 9\n",i,j,k);
-						 local_wetting_interaction += 0.5*(Phi(nq)-wetval);
-						 local_wetting_weight += 0.5;
-					}
-					nq = (k)*Nx*Ny+(j+1)*Nx+(i-1);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 10\n",i,j,k);
-						 local_wetting_interaction += 0.5*(Phi(nq)-wetval);
-						 local_wetting_weight += 0.5;
-					}
-					// xz interactions
-					nq = (k+1)*Nx*Ny+(j)*Nx+(i+1);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 11\n",i,j,k);
-						 local_wetting_interaction += 0.5*(Phi(nq)-wetval);
-						 local_wetting_weight += 0.5;
-					}
-					nq = (k-1)*Nx*Ny+(j)*Nx+(i-1);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 12\n",i,j,k);
-						 local_wetting_interaction += 0.5*(Phi(nq)-wetval);
-						 local_wetting_weight += 0.5;
-					}
-					nq = (k+1)*Nx*Ny+(j)*Nx+(i-1);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 13\n",i,j,k);
-						 local_wetting_interaction += 0.5*(Phi(nq)-wetval);
-						 local_wetting_weight += 0.5;
-					}
-					nq = (k-1)*Nx*Ny+(j)*Nx+(i+1);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 14\n",i,j,k);
-						 local_wetting_interaction += 0.5*(Phi(nq)-wetval);
-						 local_wetting_weight += 0.5;
-					}
-					// yz interactions
-					nq = (k+1)*Nx*Ny+(j+1)*Nx+(i);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 15\n",i,j,k);
-						 local_wetting_interaction += 0.5*(Phi(nq)-wetval);
-						 local_wetting_weight += 0.5;
-					}
-					nq = (k-1)*Nx*Ny+(j-1)*Nx+(i);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 16\n",i,j,k);
-						 local_wetting_interaction += 0.5*(Phi(nq)-wetval);
-						 local_wetting_weight += 0.5;
-					}
-					nq = (k+1)*Nx*Ny+(j-1)*Nx+(i);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 17\n",i,j,k);
-						 local_wetting_interaction += 0.5*(Phi(nq)-wetval);
-						 local_wetting_weight += 0.5;
-					}
-					nq = (k-1)*Nx*Ny+(j+1)*Nx+(i);
-					if ( Dm->id[nq] > 0 ) {
-						 //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 18\n",i,j,k);
-						 local_wetting_interaction += 0.5*(Phi(nq)-wetval);
-						 local_wetting_weight += 0.5;
-					}
-					/*  interaction due to this solid site*/
-					if (local_wetting_interaction == local_wetting_interaction){
-						total_wetting_interaction += 0.5*local_wetting_interaction;
-						if (local_wetting_weight > 0.0)
-							count_wetting_interaction += local_wetting_weight;
-					}
-					else{
-						//printf("Check interaction at %i %i %i \n",i,j,k);
-					}
+			}
+		}
+	}
+	
+	total_wetting_interaction = count_wetting_interaction = 0.0;
+	total_wetting_interaction_global = count_wetting_interaction_global=0.0;
+	for (k=kmin; k<kmax; k++){
+		for (j=jmin; j<Ny-1; j++){
+			for (i=imin; i<Nx-1; i++){
+				n = k*Nx*Ny + j*Nx + i;
+				// compute contribution of wetting terms (within two voxels of solid)
+				if ( Dm->id[n] > 0  && SDs(i,j,k) < 2.0 ){
+					count_wetting_interaction += 1.0;
+					total_wetting_interaction += Laplacian(i,j,k);
 				}
 			}
 		}
@@ -367,11 +269,12 @@ void SubPhase::Basic(){
 	//printf("wetting interaction = %f, count = %f\n",total_wetting_interaction,count_wetting_interaction);
 	total_wetting_interaction_global=Dm->Comm.sumReduce(  total_wetting_interaction);
 	count_wetting_interaction_global=Dm->Comm.sumReduce(  count_wetting_interaction);
-	/* normalize wetting interactions */
+	/* normalize wetting interactions  <-- Don't do this if normalizing laplacian (use solid surface area)
 	if (count_wetting_interaction > 0.0)
 		total_wetting_interaction /= count_wetting_interaction;
 	if (count_wetting_interaction_global > 0.0)
 		total_wetting_interaction_global /= count_wetting_interaction_global;
+	*/
 	
 	gwb.V=Dm->Comm.sumReduce(  wb.V);
 	gnb.V=Dm->Comm.sumReduce(  nb.V);
diff --git a/analysis/SubPhase.h b/analysis/SubPhase.h
index d2ce6b44..f3462133 100644
--- a/analysis/SubPhase.h
+++ b/analysis/SubPhase.h
@@ -68,11 +68,11 @@ public:
 	 * 		b - bulk (total)
 	 */
 	// local entities
-	phase wc,wd,wb,nc,nd,nb;
+	phase wc,wd,wb,nc,nd,nb,solid;
 	interface iwn,iwnc;
 	
 	// global entities
-	phase gwc,gwd,gwb,gnc,gnd,gnb;
+	phase gwc,gwd,gwb,gnc,gnd,gnb,gsolid;
 	interface giwn,giwnc;
 	/* fluid-solid wetting interaction */
 	double total_wetting_interaction, count_wetting_interaction;
@@ -88,6 +88,7 @@ public:
 	DoubleArray Rho_w;	// density field 
 	DoubleArray Phi;		// phase indicator field
 	DoubleArray DelPhi;		// Magnitude of Gradient of the phase indicator field
+	DoubleArray Laplacian;  // laplacian of phase indicator field
 	DoubleArray Pressure; 	// pressure field
 	DoubleArray Vel_x;		// velocity field
 	DoubleArray Vel_y;
diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index e2fca48d..0f9f1c74 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -940,7 +940,10 @@ void runAnalysis::run(int timestep, std::shared_ptr<Database> input_db, TwoPhase
  ******************************************************************/
 void runAnalysis::basic(int timestep, std::shared_ptr<Database> input_db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den)
 {
-    int N = d_N[0]*d_N[1]*d_N[2];
+    int Nx = d_N[0];
+    int Ny = d_N[1];
+    int Nz = d_N[2];
+    int N = Nx*Ny*Nz;
     NULL_USE( N );
     // Check which analysis steps we need to perform
 	auto color_db =  input_db->getDatabase( "Color" );
@@ -1036,7 +1039,7 @@ void runAnalysis::basic(int timestep, std::shared_ptr<Database> input_db, SubPha
     if (timestep%d_visualization_interval==0){
         // Write the vis files
         commWrapper comm = getComm();
-        fillHalo<double> fillData( comm.comm, d_rank_info, d_n, {1,1,1}, 0, 1 );
+        fillHalo<double> fillData( comm.comm, d_rank_info, {Nx-2,Ny-2,Nz-2}, {1,1,1}, 0, 1 );
         auto work = new IOWorkItem( timestep, input_db, d_meshData, Averages, fillData, std::move( comm ) );
         work->add_dependency(d_wait_analysis);
         work->add_dependency(d_wait_subphase);
diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 5f0e6f0e..ea1ec97f 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -238,9 +238,14 @@ void ScaLBL_ColorModel::ReadInput(){
 		}
 	}
 //	MeanFilter(Averages->SDs);
+	Minkowski Solid(Dm);
 	if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
 	CalcDist(Averages->SDs,id_solid,*Mask);
-	
+	Solid.ComputeScalar(Averages->SDs,0.0);
+	if (rank == 0) {
+		printf("Vs As Hs Xs\n");
+		printf("%.8g %.8g %.8g %.8g\n",Solid.Vi_global,Solid.Ai_global,Solid.Ji_global,Solid.Xi_global);
+	}
 	if (rank == 0) cout << "Domain set." << endl;
 	
 	Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);