diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index 3396c9c7..95e8dc02 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -16,8 +16,7 @@
 #include "common/MPI.h"
 
 
-#define USE_NEW_WRITER
-//#define WRITE_SURFACES
+#define WRITE_SURFACES
 
 /*
  * Simulator for two-phase flow in porous media
@@ -2277,27 +2276,22 @@ int main(int argc, char **argv)
 
 #ifdef WRITE_SURFACES
 			
-			sprintf(tmpstr,"vis%03d",logcount);
-			if (rank==0){
-				mkdir(tmpstr,0777);
-			}
-			MPI_Barrier(MPI_COMM_WORLD);
-			
-			FILE *WN_TRIS;
-			sprintf(LocalRankFilename,"%s/%s%s",tmpstr,"wn-tris.",LocalRankString);
-			WN_TRIS = fopen(LocalRankFilename,"wb");
-
-			FILE *NS_TRIS;
-			sprintf(LocalRankFilename,"%s/%s%s",tmpstr,"ns-tris.",LocalRankString);
-			NS_TRIS = fopen(LocalRankFilename,"wb");
-
-			FILE *WS_TRIS;
-			sprintf(LocalRankFilename,"%s/%s%s",tmpstr,"ws-tris.",LocalRankString);
-			WS_TRIS = fopen(LocalRankFilename,"wb");
-
-			FILE *WNS_PTS;
-			sprintf(LocalRankFilename,"%s/%s%s",tmpstr,"wns-crv.",LocalRankString);
-			WNS_PTS = fopen(LocalRankFilename,"wb");
+			std::shared_ptr<TriList> wn_mesh( new TriList() );
+			wn_mesh->A.reserve(8*ncubes);
+			wn_mesh->B.reserve(8*ncubes);
+			wn_mesh->C.reserve(8*ncubes);
+			std::shared_ptr<TriList> ns_mesh( new TriList() );
+			ns_mesh->A.reserve(8*ncubes);
+			ns_mesh->B.reserve(8*ncubes);
+			ns_mesh->C.reserve(8*ncubes);
+			std::shared_ptr<TriList> ws_mesh( new TriList() );
+			ws_mesh->A.reserve(8*ncubes);
+			ws_mesh->B.reserve(8*ncubes);
+			ws_mesh->C.reserve(8*ncubes);
+			std::shared_ptr<TriList> wns_mesh( new TriList() );
+			wns_mesh->A.reserve(8*ncubes);
+			wns_mesh->B.reserve(8*ncubes);
+			wns_mesh->C.reserve(8*ncubes);
 
 			for (c=0;c<ncubes;c++){
 				// Get cube from the list
@@ -2345,16 +2339,9 @@ int main(int argc, char **argv)
 					C.x += 1.0*iproc*(Nx-2);
 					C.y += 1.0*jproc*(Nx-2);
 					C.z += 1.0*kproc*(Nx-2);
-					// write the triangle
-					fwrite(&A.x,sizeof(A.x),1,WN_TRIS);
-					fwrite(&A.y,sizeof(A.y),1,WN_TRIS);
-					fwrite(&A.z,sizeof(A.z),1,WN_TRIS);
-					fwrite(&B.x,sizeof(B.x),1,WN_TRIS);
-					fwrite(&B.y,sizeof(B.y),1,WN_TRIS);
-					fwrite(&B.z,sizeof(B.z),1,WN_TRIS);
-					fwrite(&C.x,sizeof(C.x),1,WN_TRIS);
-					fwrite(&C.y,sizeof(C.y),1,WN_TRIS);
-					fwrite(&C.z,sizeof(C.z),1,WN_TRIS);
+					wn_mesh->A.push_back(A);
+					wn_mesh->B.push_back(B);
+					wn_mesh->C.push_back(C);
 				}		
 				for (int r=0;r<n_ws_tris;r++){
 					A = ws_pts(ws_tris(0,r));
@@ -2370,16 +2357,9 @@ int main(int argc, char **argv)
 					C.x += 1.0*iproc*(Nx-2);
 					C.y += 1.0*jproc*(Nx-2);
 					C.z += 1.0*kproc*(Nx-2);
-					// write the triangle
-					fwrite(&A.x,sizeof(A.x),1,WS_TRIS);
-					fwrite(&A.y,sizeof(A.y),1,WS_TRIS);
-					fwrite(&A.z,sizeof(A.z),1,WS_TRIS);
-					fwrite(&B.x,sizeof(B.x),1,WS_TRIS);
-					fwrite(&B.y,sizeof(B.y),1,WS_TRIS);
-					fwrite(&B.z,sizeof(B.z),1,WS_TRIS);
-					fwrite(&C.x,sizeof(C.x),1,WS_TRIS);
-					fwrite(&C.y,sizeof(C.y),1,WS_TRIS);
-					fwrite(&C.z,sizeof(C.z),1,WS_TRIS);	
+					ws_mesh->A.push_back(A);
+					ws_mesh->B.push_back(B);
+					ws_mesh->C.push_back(C);
 				}
 				for (int r=0;r<n_ns_tris;r++){
 					A = ns_pts(ns_tris(0,r));
@@ -2395,33 +2375,26 @@ int main(int argc, char **argv)
 					C.x += 1.0*iproc*(Nx-2);
 					C.y += 1.0*jproc*(Nx-2);
 					C.z += 1.0*kproc*(Nx-2);
-					// write the triangle
-					fwrite(&A.x,sizeof(A.x),1,NS_TRIS);
-					fwrite(&A.y,sizeof(A.y),1,NS_TRIS);
-					fwrite(&A.z,sizeof(A.z),1,NS_TRIS);
-					fwrite(&B.x,sizeof(B.x),1,NS_TRIS);
-					fwrite(&B.y,sizeof(B.y),1,NS_TRIS);
-					fwrite(&B.z,sizeof(B.z),1,NS_TRIS);
-					fwrite(&C.x,sizeof(C.x),1,NS_TRIS);
-					fwrite(&C.y,sizeof(C.y),1,NS_TRIS);
-					fwrite(&C.z,sizeof(C.z),1,NS_TRIS);
-				}
-				for (int p=0; p < n_nws_pts; p++){
-					P = nws_pts(p);
-				//	fprintf(WNS_PTS,"%f %f %f \n",P.x, P.y, P.z);
-					fwrite(&P.x,sizeof(P.x),1,WNS_PTS);
-					fwrite(&P.y,sizeof(P.y),1,WNS_PTS);
-					fwrite(&P.z,sizeof(P.z),1,WNS_PTS);
+					ns_mesh->A.push_back(A);
+					ns_mesh->B.push_back(B);
+					ns_mesh->C.push_back(C);
 				}
 			}
-			fclose(WN_TRIS);
-			fclose(NS_TRIS);
-			fclose(WS_TRIS);
-			fclose(WNS_PTS);
+
+			std::vector<MeshDataStruct> meshData(4);
+			meshData[0].meshName = "wn-tris";
+			meshData[0].mesh = wn_mesh;
+			meshData[1].meshName = "ws-tris";
+			meshData[1].mesh = ws_mesh;
+			meshData[2].meshName = "ns-tris";
+			meshData[2].mesh = ns_mesh;
+			meshData[3].meshName = "wns-tris";
+			meshData[3].mesh = wns_mesh;
+			writeData( logcount, meshData );
+
 			logcount++;
 #endif 
 		}
-
 	}
 	//************************************************************************/
 	DeviceBarrier();