diff --git a/tests/TestSphereCurvature.cpp b/tests/TestSphereCurvature.cpp
index bd6d732d..0ffff8d6 100644
--- a/tests/TestSphereCurvature.cpp
+++ b/tests/TestSphereCurvature.cpp
@@ -11,86 +11,60 @@ using namespace std;
  *
  */
 
+std::shared_ptr<Database> loadInputs( int nprocs )
+{
+  //auto db = std::make_shared<Database>( "Domain.in" );
+    auto db = std::make_shared<Database>();
+    db->putScalar<int>( "BC", 0 );
+    db->putVector<int>( "nproc", { 1, 1, 1 } );
+    db->putVector<int>( "n", { 32, 32, 32 } );
+    db->putScalar<int>( "nspheres", 1 );
+    db->putVector<double>( "L", { 1, 1, 1 } );
+    return db;
+}
+
 int main(int argc, char **argv)
 {
-	int i,j,k;
-	int Nx,Ny,Nz;
-	double Lx,Ly,Lz;
-	double fluid_isovalue=0.0;
-	double solid_isovalue=0.0;
-
-	Lx = Ly = Lz = 1.0;
-	Nx = Ny = Nz = 64;
-	DoubleArray Phase(Nx,Ny,Nz);
-	DoubleArray CubeValues(2,2,2);
-	
-	printf("Set distance map \n");
-	for (k=0; k<Nz; k++){
-		for (j=0; j<Ny; j++){
-			for (i=0; i<Nx; i++){
-				Phase(i,j,k) = sqrt((1.0*i-0.5*Nx)*(1.0*i-0.5*Nx)+(1.0*j-0.5*Ny)*(1.0*j-0.5*Ny)+(1.0*k-0.5*Nz)*(1.0*k-0.5*Nz))-0.3*Nx;
-			}
-		}
-	}
-	
-	printf("Construct local isosurface \n");
-	DECL sphere;
-	Point P1,P2,P3;
-	unsigned long int e1,e2,e3;
-	double s,s1,s2,s3;
-	double area = 0.f;
-	double Xi = 0.f;
-	double Vx,Vy,Vz,Wx,Wy,Wz,nx,ny,nz,norm;
-	for (k=0; k<Nz-1; k++){
-		for (j=0; j<Ny-1; j++){
-			for (i=0; i<Nx-1; i++){
-				sphere.LocalIsosurface(Phase,0.f,i,j,k);
-				for (unsigned long int idx=0; idx<sphere.TriangleCount; idx++){
-					e1 = sphere.Face(idx); 
-					e2 = sphere.halfedge.next(e1);
-					e3 = sphere.halfedge.next(e2);
-					P1 = sphere.vertex.coords(sphere.halfedge.v1(e1));
-					P2 = sphere.vertex.coords(sphere.halfedge.v1(e2));
-					P3 = sphere.vertex.coords(sphere.halfedge.v1(e3));
-					// compute the area
-					s1 = sqrt((P1.x-P2.x)*(P1.x-P2.x)+(P1.y-P2.y)*(P1.y-P2.y)+(P1.z-P2.z)*(P1.z-P2.z));
-					s2 = sqrt((P1.x-P3.x)*(P1.x-P3.x)+(P1.y-P3.y)*(P1.y-P3.y)+(P1.z-P3.z)*(P1.z-P3.z));
-					s3 = sqrt((P2.x-P3.x)*(P2.x-P3.x)+(P2.y-P3.y)*(P2.y-P3.y)+(P2.z-P3.z)*(P2.z-P3.z));
-					s = 0.5*(s1+s2+s3);
-					area += sqrt(s*(s-s1)*(s-s2)*(s-s3));
-					// compute the normal vector
-					Vx=P2.x-P1.x;
-					Vy=P2.y-P1.y;
-					Vz=P2.z-P1.z;
-					Wx=P3.x-P2.x;
-					Wy=P3.y-P2.y;
-					Wz=P3.z-P2.z;
-					nx = Vy*Wz-Vz*Wy;
-					ny = Vz*Wx-Vx*Wz;
-					nz = Vx*Wy-Vy*Wx;
-					norm = 1.f/sqrt(nx*nx+ny*ny+nz*nz);
-					nx *= norm;
-					ny *= norm;
-					nz *= norm;
-					// Euler characteristic (half edge rule: one face - 0.5*(three edges))
-					Xi -= 0.5;
-				}
-				// Euler characteristic -- each vertex shared by four cubes
-				Xi += 0.25*double(sphere.VertexCount);
-			}
-		}
-	}
-
-	printf("Surface area  = %f (analytical = %f) \n", area,4*3.14159*0.3*0.3*double(Nx*Nx));
-	printf("Euler characteristic  = %f (analytical = 2.0) \n",Xi);
-
-	//	printf("Mean Curvature Average =  %f, Analytical = %f \n", wn_curvature_sum/wn_area_sum, 2.0/rad[0]/101 );
-
+	MPI_Init(&argc,&argv);
+	MPI_Comm comm = MPI_COMM_WORLD;
+	int rank = MPI_WORLD_RANK();
+	int nprocs = MPI_WORLD_SIZE();
 	int toReturn = 0;
-/*	if ( fabs(wn_curvature_sum/wn_area_sum -2.0/rad[0]/101)*rad[0]*101.0*0.5 > 0.01 ){
-		toReturn = 1;
-		printf("Mean curvature test error exceeds relative error tolerance \n ");
+	{
+		int i,j,k;
+
+		// Load inputs
+		auto db = loadInputs( nprocs );
+		int Nx = db->getVector<int>( "n" )[0];
+		int Ny = db->getVector<int>( "n" )[1];
+		int Nz = db->getVector<int>( "n" )[2];
+		int nprocx = db->getVector<int>( "nproc" )[0];
+		int nprocy = db->getVector<int>( "nproc" )[1];
+		int nprocz = db->getVector<int>( "nproc" )[2];
+		std::shared_ptr<Domain> Dm(new Domain(db,comm));
+
+		DoubleArray Phase(Nx,Ny,Nz);
+		
+        Minkowski sphere(new Minkowski(Dm));
+
+		printf("Set distance map \n");
+		for (k=0; k<Nz; k++){
+			for (j=0; j<Ny; j++){
+				for (i=0; i<Nx; i++){
+					Phase(i,j,k) = sqrt((1.0*i-0.5*Nx)*(1.0*i-0.5*Nx)+(1.0*j-0.5*Ny)*(1.0*j-0.5*Ny)+(1.0*k-0.5*Nz)*(1.0*k-0.5*Nz))-0.3*Nx;
+				}
+			}
+		}
+
+		printf("Construct local isosurface \n");
+		sphere.ComputeScalar(Phase,0.f);
+
+		printf("Surface area  = %f (analytical = %f) \n", sphere.Ai,4*3.14159*0.3*0.3*double(Nx*Nx));
+		printf("Euler characteristic  = %f (analytical = 2.0) \n",sphere.Xi);
+
+
 	}
-	*/
+	MPI_Barrier(comm);
+	MPI_Finalize();
 	return toReturn;
 }