diff --git a/tests/GenerateSphereTest.cpp b/tests/GenerateSphereTest.cpp
index 48b4840f..fb17464d 100644
--- a/tests/GenerateSphereTest.cpp
+++ b/tests/GenerateSphereTest.cpp
@@ -328,13 +328,6 @@ int main(int argc, char **argv)
 		int nprocx,nprocy,nprocz;
 		int iproc,jproc,kproc;
 		//*****************************************
-		// MPI ranks for all 18 neighbors
-		//**********************************
-		int rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z;
-		int rank_xy,rank_XY,rank_xY,rank_Xy;
-		int rank_xz,rank_XZ,rank_xZ,rank_Xz;
-		int rank_yz,rank_YZ,rank_yZ,rank_Yz;
-		//**********************************
 
 		if (rank == 0){
 			printf("********************************************************\n");
@@ -343,72 +336,36 @@ int main(int argc, char **argv)
 		}
 
 		// Variables that specify the computational domain  
-		string FILENAME;
+		string filename;
 		int Nx,Ny,Nz;		// local sub-domain size
 		int nspheres;		// number of spheres in the packing
 		double Lx,Ly,Lz;	// Domain length
 		double D = 1.0;		// reference length for non-dimensionalization
 
 		int i,j,k,n;
+		filename = argv[1];
+		auto db = std::make_shared<Database>( filename );
+		auto domain_db = db->getDatabase( "Domain" );
 		
-		if (argc > 1)
-			nspheres=atoi(argv[1]);
-		else nspheres=0;
-		
-		if (rank==0){
-			//.......................................................................
-			// Reading the domain information file
-			//.......................................................................
-			ifstream domain("Domain.in");
-			domain >> nprocx;
-			domain >> nprocy;
-			domain >> nprocz;
-			domain >> Nx;
-			domain >> Ny;
-			domain >> Nz;
-			domain >> Lx;
-			domain >> Ly;
-			domain >> Lz;
-			//.......................................................................
-		}
-		// **************************************************************
-		// Broadcast simulation parameters from rank 0 to all other procs
-		MPI_Barrier(comm);
-		//.................................................
-		// Computational domain
-		MPI_Bcast(&Nx,1,MPI_INT,0,comm);
-		MPI_Bcast(&Ny,1,MPI_INT,0,comm);
-		MPI_Bcast(&Nz,1,MPI_INT,0,comm);
-		MPI_Bcast(&nprocx,1,MPI_INT,0,comm);
-		MPI_Bcast(&nprocy,1,MPI_INT,0,comm);
-		MPI_Bcast(&nprocz,1,MPI_INT,0,comm);
-		MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm);
-		MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm);
-		MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm);
-		//.................................................
-		MPI_Barrier(comm);
-
-		// **************************************************************
-
-		if (nprocs != nprocx*nprocy*nprocz){
-			printf("nprocx =  %i \n",nprocx);
-			printf("nprocy =  %i \n",nprocy);
-			printf("nprocz =  %i \n",nprocz);
-			INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!");
-		}
-
-		InitializeRanks( rank, nprocx, nprocy, nprocz, iproc, jproc, kproc, 
-				rank_x, rank_y, rank_z, rank_X, rank_Y, rank_Z,
-				rank_xy, rank_XY, rank_xY, rank_Xy, rank_xz, rank_XZ, rank_xZ, rank_Xz,
-				rank_yz, rank_YZ, rank_yZ, rank_Yz );
-
-		MPI_Barrier(comm);
+		auto Dm  = std::shared_ptr<Domain>(new Domain(domain_db,comm));      // full domain for analysis
+		Nx = Dm->Nx;
+		Ny = Dm->Ny;
+		Nz = Dm->Nz;
+		Lx = Dm->Lx;
+		Ly = Dm->Ly;
+		Lz = Dm->Lz;
+		iproc = Dm->iproc();
+		jproc = Dm->jproc();
+		kproc = Dm->kproc();
+		nprocx = Dm->nprocx();
+		nprocy = Dm->nprocy();
+		nprocz = Dm->nprocz();
+	    nspheres = domain_db->getScalar<int>( "nspheres");
 
+		printf("Set domain \n");
 		int BoundaryCondition=1;
-		Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
-
-		Nz += 2;
-		Nx = Ny = Nz;	// Cubic domain
+		//Nz += 2;
+		//Nx = Ny = Nz;	// Cubic domain
 		int N = Nx*Ny*Nz;
 
 		// Define Dm.Communication sub-domain -- everywhere
@@ -416,11 +373,11 @@ int main(int argc, char **argv)
 			for (int j=0; j<Ny; j++){
 				for (int i=0; i<Nx; i++){
 					n = k*Nx*Ny+j*Nx+i;
-					Dm.id[n] = 1;
+					Dm->id[n] = 1;
 				}
 			}
 		}
-		Dm.CommInit();
+		Dm->CommInit();
 
 		if (rank==0) printf("Number of nodes per side = %i \n", Nx);
 		if (rank==0) printf("Total Number of nodes = %i \n", N);
@@ -490,8 +447,6 @@ int main(int argc, char **argv)
 		SignedDistance(SignDist.data(),nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,
 				iproc,jproc,kproc,nprocx,nprocy,nprocz);
 		//.......................................................................
-
-
 		// Assign the phase ID field based on the signed distance
 		//.......................................................................
 		for (k=0;k<Nz;k++){
@@ -526,7 +481,7 @@ int main(int argc, char **argv)
 		// Run Morphological opening to initialize 50% saturation
 		double SW=0.50;
 		if (rank==0) printf("MorphOpen: Initializing with saturation %f \n",SW);
-		MorphOpen(SignDist, id, Dm, Nx, Ny, Nz, rank, SW);
+		MorphOpen(SignDist, id, *Dm, Nx, Ny, Nz, rank, SW);
 
 		//.........................................................
 		// don't perform computations at the eight corners