diff --git a/common/TwoPhase.cpp b/common/TwoPhase.cpp
index a505612a..794bacc8 100644
--- a/common/TwoPhase.cpp
+++ b/common/TwoPhase.cpp
@@ -507,7 +507,7 @@ void TwoPhase::ComponentAverages()
 			for (i=1; i<Nx-1; i++){
 				LabelWP=GetCubeLabel(i,j,k,Label_WP);
 				LabelNWP=GetCubeLabel(i,j,k,Label_NWP);
-				
+
 				n_nw_pts=n_ns_pts=n_ws_pts=n_nws_pts=n_local_sol_pts=n_local_nws_pts=0;
 				n_nw_tris=n_ns_tris=n_ws_tris=n_nws_seg=n_local_sol_tris=0;
 
@@ -592,12 +592,12 @@ void TwoPhase::ComponentAverages()
 				 */
 				if (n_nw_pts+n_ns_pts > 0){
 
-				/*	double euler;
+					/*	double euler;
 					euler =  geomavg_EulerCharacteristic(nw_pts,nw_tris,n_nw_pts,n_nw_tris,i,j,k);
 					euler += geomavg_EulerCharacteristic(ns_pts,ns_tris,n_ns_pts,n_ns_tris,i,j,k);
 					// adjust for double-counted vertices and edges from the common curve
 					if (n_nws_pts > 0) euler += 1.0*n_nws_pts;
-				 */
+
 					int nvert = n_nw_pts+n_ns_pts-n_nws_pts;
 					int nside = 2*nvert-3;
 					int nface = nvert-2;
@@ -678,6 +678,8 @@ void TwoPhase::ComponentAverages()
 					//ComponentAverages_NWP(NVERT,LabelNWP) += nvert;
 					//ComponentAverages_NWP(NSIDE,LabelNWP) += nside;
 					//ComponentAverages_NWP(NFACE,LabelNWP) += nface;
+					 *
+					 */
 				}
 
 				//...........................................................................
@@ -778,6 +780,54 @@ void TwoPhase::ComponentAverages()
 				}
 				//...........................................................................
 
+				// Compute the Euler characteristic
+				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);
+				if (n_nw_pts > 0 ){
+					int nvert = n_nw_pts;
+					int nside = 2*nvert-3;
+					int nface = nvert-2;
+					//...........................................................
+					// Check that this point is not on a previously computed face
+					// Note direction that the marching cubes algorithm marches
+					// In parallel, other sub-domains fill in the lower boundary
+					for (int p=0; p<n_nw_pts; p++){
+						Point PT = nw_pts(p);
+						if (PT.x - double(i) < 1e-12)		nvert-=1;
+						else if (PT.y - double(j) < 1e-12) 	nvert-=1;
+						else if (PT.z - double(k) < 1e-12) 	nvert-=1;
+					}
+					// Remove previously computed edges
+					for (int p=0; p<n_nw_tris; p++){
+						Point A = nw_pts(nw_tris(0,p));
+						Point B = nw_pts(nw_tris(1,p));
+						Point C = nw_pts(nw_tris(2,p));
+
+						// Check side A-B
+						bool newside = true;
+						if (A.x - double(i) < 1e-12 && B.x - double(i) < 1e-12) newside=false;
+						if (A.y - double(j) < 1e-12 && B.y - double(j) < 1e-12) newside=false;
+						if (A.z - double(k) < 1e-12 && B.z - double(k) < 1e-12) newside=false;
+						if (!newside) nside-=1;
+
+						// Check side A-C
+						newside = true;
+						if (A.x - double(i)< 1e-12 && C.x - double(i) < 1e-12) newside=false;
+						if (A.y - double(j)< 1e-12 && C.y - double(j) < 1e-12) newside=false;
+						if (A.z - double(k)< 1e-12 && C.z - double(k) < 1e-12) newside=false;
+						if (!newside) nside-=1;
+
+						// Check side B-C
+						newside = true;
+						if (B.x - double(i) < 1e-12 && C.x - double(i) < 1e-12) newside=false;
+						if (B.y - double(j) < 1e-12 && C.y - double(j) < 1e-12) newside=false;
+						if (B.z - double(k) < 1e-12 && C.z - double(k) < 1e-12) newside=false;
+						if (!newside) nside-=1;
+
+					}
+					int euler=nvert-nside+nface; // euler characteristic for the cube
+					ComponentAverages_NWP(EULER,LabelNWP) += 1.0*euler;
+				}
 			}
 		}
 	}
diff --git a/common/pmmc.h b/common/pmmc.h
index 8637ada9..15425a0c 100644
--- a/common/pmmc.h
+++ b/common/pmmc.h
@@ -1837,6 +1837,320 @@ 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,
+				DTMutableList<Point> &nw_pts, int &n_nw_pts, IntArray &nw_tris,
+				int &n_nw_tris)
+{
+	int N = 0;		// n will be the number of vertices in this grid cell only
+    Point P;
+	Point pt;
+    Point PlaceHolder;
+    int m;
+    int o;
+    int p;
+
+	// Go over each corner -- check to see if the corners are themselves vertices
+	//1
+	if (A(i,j,k) == v){
+		P.x = i;
+		P.y = j;
+		P.z = k;
+		nw_pts(n_nw_pts++) = P;
+		N++;
+	}
+	//2
+	if (A(i+1,j,k) == v){
+		P.x = i+1;
+		P.y = j;
+		P.z = k;
+		nw_pts(n_nw_pts++) = P;
+		N++;
+	}
+	//3
+	if (A(i+1,j+1,k) == v){
+		P.x = i+1;
+		P.y = j+1;
+		P.z = k;
+		nw_pts(n_nw_pts++) = P;
+		N++;
+	}
+	//4
+	if (A(i,j+1,k) == v){
+		P.x = i;
+		P.y = j+1;
+		P.z = k;
+		nw_pts(n_nw_pts++) = P;
+		N++;
+	}
+	//5
+	if (A(i,j,k+1) == v){
+		P.x = i;
+		P.y = j;
+		P.z = k+1;
+		nw_pts(n_nw_pts++) = P;
+		N++;
+	}
+	//6
+	if (A(i+1,j,k+1) == v){
+		P.x = i+1;
+		P.y = j;
+		P.z = k+1;
+		nw_pts(n_nw_pts++) = P;
+		N++;
+	}
+	//7
+	if (A(i+1,j+1,k+1) == v){
+		P.x = i+1;
+		P.y = j+1;
+		P.z = k+1;
+		nw_pts(n_nw_pts++) = P;
+		N++;
+	}
+	//8
+	if (A(i,j+1,k+1) == v){
+		P.x = i;
+		P.y = j+1;
+		P.z = k+1;
+
+		nw_pts(n_nw_pts++) = P;
+		N++;
+	}
+
+    // Go through each side, compute P for sides of box spiraling up
+
+
+//	float val;
+	if ((A(i,j,k)-v)*(A(i+1,j,k)-v) < 0)
+	{
+		// If both points are in the fluid region
+		if (A(i,j,k) != 0 && A(i+1,j,k) != 0){
+			P.x = i + (A(i,j,k)-v)/(A(i,j,k)-A(i+1,j,k));
+			P.y = j;
+			P.z = k;
+
+			nw_pts(n_nw_pts++) =  P;
+			N++;
+		}
+	}
+
+	if ((A(i+1,j,k)-v)*(A(i+1,j+1,k)-v) < 0)
+	{
+		if ( A(i+1,j,k) != 0 && A(i+1,j+1,k) != 0 ){
+			P.x = i+1;
+			P.y = j + (A(i+1,j,k)-v)/(A(i+1,j,k)-A(i+1,j+1,k));
+			P.z = k;
+
+			if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){ // P is a new vertex (not counted twice)
+				nw_pts(n_nw_pts++) =  P;
+				N++;
+			}
+		}
+	}
+
+	if ((A(i+1,j+1,k)-v)*(A(i,j+1,k)-v) < 0 )
+	{
+		if ( A(i+1,j+1,k) != 0 && A(i,j+1,k) != 0 ){
+			P.x = i + (A(i,j+1,k)-v) / (A(i,j+1,k)-A(i+1,j+1,k));
+			P.y = j+1;
+			P.z = k;
+			if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){ // P is a new vertex (not counted twice)
+				nw_pts(n_nw_pts++) =  P;
+				N++;
+			}
+		}
+	}
+
+	//4
+	if ((A(i,j+1,k)-v)*(A(i,j,k)-v) < 0 )
+	{
+		if (A(i,j+1,k) != 0 && A(i,j,k) != 0 ){
+			P.x = i;
+			P.y = j + (A(i,j,k)-v) / (A(i,j,k)-A(i,j+1,k));
+			P.z = k;
+			if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){ // P is a new vertex (not counted twice)
+				nw_pts(n_nw_pts++) =  P;
+				N++;
+			}
+		}
+	}
+
+	//5
+	if ((A(i,j,k)-v)*(A(i,j,k+1)-v) < 0 )
+	{
+		if ( A(i,j,k) != 0 && A(i,j,k+1) != 0 ){
+			P.x = i;
+			P.y = j;
+			P.z = k + (A(i,j,k)-v) / (A(i,j,k)-A(i,j,k+1));
+
+			if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){ // P is a new vertex (not counted twice)
+				nw_pts(n_nw_pts++) =  P;
+				N++;
+			}
+		}
+	}
+
+	//6
+	if ((A(i+1,j,k)-v)*(A(i+1,j,k+1)-v) < 0 )
+	{
+		if ( A(i+1,j,k) != 0 && A(i+1,j,k+1) != 0 ){
+			P.x = i+1;
+			P.y = j;
+			P.z = k + (A(i+1,j,k)-v) / (A(i+1,j,k)-A(i+1,j,k+1));
+			if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){
+				nw_pts(n_nw_pts++) =  P;
+				N++;
+			}
+		}
+	}
+
+	//7
+	if ((A(i+1,j+1,k)-v)*(A(i+1,j+1,k+1)-v) < 0 )
+	{
+		if ( A(i+1,j+1,k) != 0 && A(i+1,j+1,k+1) != 0 ){
+			P.x = i+1;
+			P.y = j+1;
+			P.z = k + (A(i+1,j+1,k)-v) / (A(i+1,j+1,k)-A(i+1,j+1,k+1));
+			if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){
+				nw_pts(n_nw_pts++) =  P;
+				N++;
+			}
+		}
+	}
+
+	//8
+	if ((A(i,j+1,k)-v)*(A(i,j+1,k+1)-v) < 0 )
+	{
+		if ( A(i,j+1,k) != 0 && A(i,j+1,k+1) != 0 ){
+			P.x = i;
+			P.y = j+1;
+			P.z = k + (A(i,j+1,k)-v) / (A(i,j+1,k)-A(i,j+1,k+1));
+			if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){
+				nw_pts(n_nw_pts++) =  P;
+				N++;
+			}
+		}
+	}
+
+	//9
+	if ((A(i,j,k+1)-v)*(A(i+1,j,k+1)-v) < 0 )
+	{
+		if ( A(i,j,k+1) != 0 && A(i+1,j,k+1) != 0 ){
+			P.x = i + (A(i,j,k+1)-v) / (A(i,j,k+1)-A(i+1,j,k+1));
+			P.y = j;
+			P.z = k+1;
+			if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){
+				nw_pts(n_nw_pts++) =  P;
+				N++;
+			}
+		}
+	}
+
+	//10
+	if ((A(i+1,j,k+1)-v)*(A(i+1,j+1,k+1)-v) < 0 )
+	{
+		if ( A(i+1,j,k+1) != 0 && A(i+1,j+1,k+1) != 0 ){
+			P.x = i+1;
+			P.y = j + (A(i+1,j,k+1)-v) / (A(i+1,j,k+1)-A(i+1,j+1,k+1));
+			P.z = k+1;
+			if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){
+				nw_pts(n_nw_pts++) =  P;
+				N++;
+			}
+		}
+	}
+
+	//11
+	if ((A(i+1,j+1,k+1)-v)*(A(i,j+1,k+1)-v) < 0 )
+	{
+		if ( A(i+1,j+1,k+1) != 0 && A(i,j+1,k+1) != 0 ){
+			P.x = i+(A(i,j+1,k+1)-v) / (A(i,j+1,k+1)-A(i+1,j+1,k+1));
+			P.y = j+1;
+			P.z = k+1;
+			if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){
+				nw_pts(n_nw_pts++) =  P;
+				N++;
+			}
+		}
+	}
+
+	//12
+	if ((A(i,j+1,k+1)-v)*(A(i,j,k+1)-v) < 0 )
+	{
+		if ( A(i,j+1,k+1) != 0 && A(i,j,k+1) != 0 ){
+			P.x = i;
+			P.y = j + (A(i,j,k+1)-v) / (A(i,j,k+1)-A(i,j+1,k+1));
+			P.z = k+1;
+			if (vertexcheck(P, N, n_nw_pts, nw_pts) == 1){
+				nw_pts(n_nw_pts++) =  P;
+				N++;
+			}
+		}
+	}
+
+
+    // Assemble the triangles as long as points are found
+    if (N > 0){
+    	for (m = n_nw_pts-N; m < n_nw_pts-2; m++) {
+    		for (o = m+2; o < n_nw_pts-1; o++) {
+    			if (ShareSide(nw_pts(m), nw_pts(o)) == 1) {
+    				PlaceHolder = nw_pts(m+1);
+    				nw_pts(m+1) = nw_pts(o);
+    				nw_pts(o) = PlaceHolder;
+    			}
+    		}
+
+    		// make sure other neighbor of vertex 1 is in last spot
+    		if (m == n_nw_pts-N){
+    			for (p = m+2; p < n_nw_pts-1; p++){
+    				if (ShareSide(nw_pts(m), nw_pts(p)) == 1){
+    					PlaceHolder = nw_pts(n_nw_pts-1);
+    					nw_pts(n_nw_pts-1) = nw_pts(p);
+    					nw_pts(p) = PlaceHolder;
+    				}
+    			}
+    		}
+    		if ( ShareSide(nw_pts(n_nw_pts-2), nw_pts(n_nw_pts-3)) != 1 ){
+    			if (ShareSide( nw_pts(n_nw_pts-3), nw_pts(n_nw_pts-1)) == 1 &&
+    					ShareSide( nw_pts(n_nw_pts-N),nw_pts(n_nw_pts-2)) == 1 ){
+    				PlaceHolder = nw_pts(n_nw_pts-2);
+    				nw_pts(n_nw_pts-2) = nw_pts(n_nw_pts-1);
+    				nw_pts(n_nw_pts-1) = PlaceHolder;
+    			}
+    		}
+    		if ( ShareSide(nw_pts(n_nw_pts-1), nw_pts(n_nw_pts-2)) != 1 ){
+    			if (ShareSide( nw_pts(n_nw_pts-3), nw_pts(n_nw_pts-1)) == 1 &&
+    					ShareSide(nw_pts(n_nw_pts-4),nw_pts(n_nw_pts-2)) == 1 ){
+    				PlaceHolder = nw_pts(n_nw_pts-3);
+    				nw_pts(n_nw_pts-3) = nw_pts(n_nw_pts-2);
+    				nw_pts(n_nw_pts-2) = PlaceHolder;
+    			}
+    			if (ShareSide( nw_pts(n_nw_pts-N+1), nw_pts(n_nw_pts-3)) == 1 &&
+    					ShareSide(nw_pts(n_nw_pts-1),nw_pts(n_nw_pts-N+1)) == 1 ){
+    				PlaceHolder = nw_pts(n_nw_pts-2);
+    				nw_pts(n_nw_pts-2) = nw_pts(n_nw_pts-N+1);
+    				nw_pts(n_nw_pts-N+1) = PlaceHolder;
+    			}
+    		}
+    		if ( ShareSide(nw_pts(n_nw_pts-N), nw_pts(n_nw_pts-N+1)) != 1 ){
+    			if (ShareSide( nw_pts(n_nw_pts-N), nw_pts(n_nw_pts-2)) == 1 &&
+    					ShareSide(nw_pts(n_nw_pts-1), nw_pts(n_nw_pts-N+1)) == 1){
+    				PlaceHolder = nw_pts(n_nw_pts-1);
+    				nw_pts(n_nw_pts-1) = nw_pts(n_nw_pts-N);
+    				nw_pts(n_nw_pts-N) = PlaceHolder;
+    			}
+    		}
+    	}
+
+    	// *    *    *   ESTABLISH TRIANGLE CONNECTIONS	   *    *    *
+
+    	for (p=n_nw_pts-N+2; p<n_nw_pts; p++){
+    		nw_tris(0,n_nw_tris) = n_nw_pts-N;
+    		nw_tris(1,n_nw_tris) = p-1;
+    		nw_tris(2,n_nw_tris) = p;
+    		n_nw_tris++;
+    	}
+    }
+}
 //-------------------------------------------------------------------------------
 inline void MC( DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, int &k,
 				DTMutableList<Point> &nw_pts, int &n_nw_pts, IntArray &nw_tris,
@@ -4037,14 +4351,14 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray
 }
 //--------------------------------------------------------------------------------------------------------
 inline double geomavg_EulerCharacteristic(DTMutableList<Point> &Points, IntArray &Triangles,
-		int npts, int ntris, int i, int j, int k){
+		int &npts, int &ntris, int &i, int &j, int &k){
 
 	// Compute the Euler characteristic for triangles in a cube
 	// Exclude edges and vertices shared with between multiple cubes
 	double EulerChar;
 	int nvert=npts;
-	int nside=2*npts-3;
-	int nface=npts-2;
+	int nside=2*vert-3;
+	int nface=nvert-2;
 	//if (ntris != nface){
 	//	nface = ntris;
 	//	nside =
@@ -4087,7 +4401,7 @@ inline double geomavg_EulerCharacteristic(DTMutableList<Point> &Points, IntArray
 		if (!newside) nside-=1;
 
 	}
-	EulerChar = double(nvert - nside + nface);
+	EulerChar = 1.0*(nvert - nside + nface);
 	return EulerChar;
 }