From 26c9c6bf01742e488bf43c3afb9ce6a492e6848c Mon Sep 17 00:00:00 2001
From: Olivier Coulaud <Olivier.Coulaud@inria.fr>
Date: Wed, 9 Apr 2014 15:46:22 +0200
Subject: [PATCH] Add vtk, vtp and cosmos format for visualization

---
 Examples/generateDistributions.cpp  | 162 ++++++++++++++++------------
 Src/Utils/FGenerateDistribution.hpp |  87 ++++++++++++++-
 2 files changed, 177 insertions(+), 72 deletions(-)

diff --git a/Examples/generateDistributions.cpp b/Examples/generateDistributions.cpp
index f5789d19..b6e267ce 100644
--- a/Examples/generateDistributions.cpp
+++ b/Examples/generateDistributions.cpp
@@ -33,6 +33,7 @@
 //!     \param   -filename name: generic name for files (without extension) and save data
 //!                 with following format in name.xxx or name.bin in -bin (not yet implemented) is set
 //!       \param   -visu save output in filename.txt
+//!       \param  -visufmt format for the visu file (vtk, vtp, cvs or cosmo). vtp is the default
 //!
 //!  <b> Geometry arguments:</b>
 //!      \param  -unitSphere uniform distribution on unit sphere
@@ -58,32 +59,33 @@
 
 //
 //
-		void genDistusage() {
-			std::cout << "Driver to generate N points (non)uniformly distributed on a given geometry"
- << std::endl;
-           std::cout <<	 "Options  "<< std::endl
-		                 <<     "   -help       to see the parameters    " << std::endl
-		                 <<     "   -N       The number of points in the distribution    " << std::endl
-			              <<    std::endl
-			             <<     "    Distributions   " << std::endl
-			             <<     "        Uniform on    " << std::endl
-			             <<     "             -unitSphere  uniform distribution on unit sphere" <<std::endl
-			             <<     "             -sphere  uniform distribution on  sphere of radius given by" <<std::endl
-			             <<     "                     -radius  R - default value for R is 2.0" <<std::endl
-			             <<     "             -prolate ellipsoid with aspect ratio a:a:c" <<std::endl
-			             <<     "                     -ar a:a:c   with  c > a > 0" <<std::endl<<std::endl
-			             <<     "        Non Uniform on    " << std::endl
-			             <<     "             -ellipsoid non uniform distribution on  an ellipsoid of aspect ratio given by" <<std::endl
-			             <<     "                     -ar a:b:c   with a, b and c > 0" <<std::endl
-			             <<     "             -plummer (Highly non unuiform) plummer distrinution (astrophysics)"<<std::endl
-			             <<     "                     -radius  R - default value 10.0" <<std::endl
-			             <<     "    Physical values" <<std::endl
-			             <<     "             -charge generate physical values between -1 and 1 otherwise generate between 0 and 1		" <<std::endl<<std::endl
-			             <<     "     Output " << std::endl
-			             <<     "             -filename name: generic name for files (without extension) and save data" <<std::endl
-			             <<     "                     with following format in name.xxx or name.bin in -bin is set" <<std::endl
-			             <<     "             -visu save output in name.txt" <<std::endl;
-		}
+void genDistusage() {
+	std::cout << "Driver to generate N points (non)uniformly distributed on a given geometry"
+			<< std::endl;
+	std::cout <<	 "Options  "<< std::endl
+			<<     "   -help       to see the parameters    " << std::endl
+			<<     "   -N       The number of points in the distribution    " << std::endl
+			<<    std::endl
+			<<     "    Distributions   " << std::endl
+			<<     "        Uniform on    " << std::endl
+			<<     "             -unitSphere  uniform distribution on unit sphere" <<std::endl
+			<<     "             -sphere  uniform distribution on  sphere of radius given by" <<std::endl
+			<<     "                     -radius  R - default value for R is 2.0" <<std::endl
+			<<     "             -prolate ellipsoid with aspect ratio a:a:c" <<std::endl
+			<<     "                     -ar a:a:c   with  c > a > 0" <<std::endl<<std::endl
+			<<     "        Non Uniform on    " << std::endl
+			<<     "             -ellipsoid non uniform distribution on  an ellipsoid of aspect ratio given by" <<std::endl
+			<<     "                     -ar a:b:c   with a, b and c > 0" <<std::endl
+			<<     "             -plummer (Highly non unuiform) plummer distrinution (astrophysics)"<<std::endl
+			<<     "                     -radius  R - default value 10.0" <<std::endl
+			<<     "    Physical values" <<std::endl
+			<<     "             -charge generate physical values between -1 and 1 otherwise generate between 0 and 1		" <<std::endl<<std::endl
+			<<     "     Output " << std::endl
+			<<     "             -filename name: generic name for files (without extension) and save data" <<std::endl
+			<<     "                     with following format in name.xxx or name.bin in -bin is set" <<std::endl
+			<<     "             -visu save output in name.txt" <<std::endl
+			<<     "             -visufmt  vtk, vtp, cosmo or cvs format " <<std::endl;
+}
 
 int main(int argc, char ** argv){
 	//
@@ -141,7 +143,7 @@ int main(int argc, char ** argv){
 		if(A != B){
 			std::cerr << " A /= B in prolate sllipsoide A =B. Your aspect ratio: "<< aspectRatio<<std::endl;
 		}
-			std::cout << "A: "<<A<<" B "<< B << " C: " << C<<std::endl;
+		std::cout << "A: "<<A<<" B "<< B << " C: " << C<<std::endl;
 		unifRandonPointsOnProlate(NbPoints,A,C,particles);
 		BoxWith =  C;
 	}
@@ -155,52 +157,76 @@ int main(int argc, char ** argv){
 		nonunifRandonPointsOnElipsoid(NbPoints,A,B,C,particles);
 		BoxWith =  FMath::Max( A,FMath::Max( B,C)) ;
 	}
-		else if(FParameters::existParameter(argc, argv, "-plummer")){
-			const FReal Radius  = FParameters::getValue(argc,argv,"-radius",  10.0);
-			unifRandonPlummer(NbPoints, Radius, sum, particles) ;
-			BoxWith = 2.0*Radius ;
-		}
+	else if(FParameters::existParameter(argc, argv, "-plummer")){
+		const FReal Radius  = FParameters::getValue(argc,argv,"-radius",  10.0);
+		unifRandonPlummer(NbPoints, Radius, sum, particles) ;
+		BoxWith = 2.0*Radius ;
+	}
 
-		else {
-			std::cout << "Bad geometry option"<< std::endl;
-			exit(-1) ;
+	else {
+		std::cout << "Bad geometry option"<< std::endl;
+		exit(-1) ;
+	}
+
+	if(FParameters::existParameter(argc, argv, "-visu")){
+		std::string visufile(""), fmt(FParameters::getStr(argc,argv,"-visufmt",   "vtp"));
+		if( fmt == "vtp" ){
+			visufile = genericFileName + ".vtp" ;
 		}
-		if(FParameters::existParameter(argc, argv, "-visu")){
-			std::ofstream file( genericFileName + ".txt", std::ofstream::out);
-			if(!file) {
-				std::cout << "Cannot open file."<< std::endl;
-				exit(-1)	;
-			}	//
-			//
-			// Export data in cvs format
-			//
-			std::cout << "Writes in CVS format  (visualization) in file "<< genericFileName + ".txt" <<std::endl ;
-			exportCVS( file, NbPoints, particles)  ;
-			//
-			// Export data in vtk format
-			//
+		else	if( fmt == "vtk" ){
+			visufile = genericFileName + ".vtk" ;
 		}
-		//
-		//  Generate file for ScalFMM Loader
-		//
-		std::ofstream outfile( genericFileName + ".fma", std::ofstream::out);
-		if(!outfile) {
-			std::cout << "Cannot open file."<< std::endl;
-			exit(-1)	 ;
+		else if( fmt == "cosmo" ){
+			visufile = genericFileName + ".cosmo" ;
 		}
-		BoxWith += 2*extraRadius ;
-		std::cout << "Writes in FMA format  in file "<< genericFileName + ".fma" <<std::endl ;
-		std::cout << " Points are in a cube of size  "<< BoxWith << "  Centered in the Origin"<<std::endl;
-		//
-		outfile << 	NbPoints << "  " << BoxWith << "   0.0   0.0  0.0 " << std::endl;
-		j=0;
-		for(int i = 0 ; i< NbPoints; ++i, j+=4){
-			outfile <<    particles[j]    << "       "    <<   particles[j+1]    << "       "   <<   particles[j+2]    << "       "   <<   particles[j+3]   <<std::endl;
+		else {
+			visufile =   genericFileName + ".csv" ;
 		}
+		std::ofstream file( visufile, std::ofstream::out);
+		if(!file) {
+			std::cout << "Cannot open file."<< std::endl;
+			exit(-1)	;
+		}	//
 		//
-		delete particles ;
-
+		// Export data in cvs format
 		//
-		return 1;
-
+		if( fmt == "vtp" ){
+			std::cout << "Writes in XML VTP format  (visualization) in file "<< visufile <<std::endl ;
+			exportVTKxml( file, NbPoints, particles)  ;
+		}
+		else		if( fmt == "vtk" ){
+			std::cout << "Writes in VTK format  (visualization) in file "<< visufile <<std::endl ;
+			exportVTK( file, NbPoints, particles)  ;
+		}
+		else if( fmt == "cosmo" ){
+			std::cout << "Writes in COSMO format  (visualization) in file "<< visufile <<std::endl ;
+			exportCOSMOS( file, NbPoints, particles)  ;
+		}
+		else {
+			std::cout << "Writes in CVS format  (visualization) in file "<<visufile<<std::endl ;
+			exportCVS( file, NbPoints, particles)  ;
+		}
+	}
+	//
+	//  Generate file for ScalFMM Loader
+	//
+	std::ofstream outfile( genericFileName + ".fma", std::ofstream::out);
+	if(!outfile) {
+		std::cout << "Cannot open file."<< std::endl;
+		exit(-1)	 ;
+	}
+	BoxWith += 2*extraRadius ;
+	std::cout << "Writes in FMA format  in file "<< genericFileName + ".fma" <<std::endl ;
+	std::cout << " Points are in a cube of size  "<< BoxWith << "  Centered in the Origin"<<std::endl;
+	//
+	outfile << 	NbPoints << "  " << BoxWith << "   0.0   0.0  0.0 " << std::endl;
+	j=0;
+	for(int i = 0 ; i< NbPoints; ++i, j+=4){
+		outfile <<    particles[j]    << "       "    <<   particles[j+1]    << "       "   <<   particles[j+2]    << "       "   <<   particles[j+3]   <<std::endl;
 	}
+	//
+	delete particles ;
+
+	//
+	return 1;
+}
diff --git a/Src/Utils/FGenerateDistribution.hpp b/Src/Utils/FGenerateDistribution.hpp
index 311a09c6..cebec439 100644
--- a/Src/Utils/FGenerateDistribution.hpp
+++ b/Src/Utils/FGenerateDistribution.hpp
@@ -181,7 +181,7 @@ void unifRandonPlummer(const int N , const FReal R, const FReal M, FReal * point
 	}
 
 	std::cout << "Total tested points: "<< cpt << " % of rejected points: "
-			      <<100*static_cast<FReal>(cpt-N)/cpt << " %" <<std::endl;
+			<<100*static_cast<FReal>(cpt-N)/cpt << " %" <<std::endl;
 
 } ;
 //! \fn void exportCVS(std::ofstream& file, const int N, const FReal * particles )
@@ -199,14 +199,93 @@ void exportCVS(std::ofstream& file, const int N, const FReal * particles ){
 		file <<    particles[j]    << " , "    <<   particles[j+1]    << " , "   <<   particles[j+2]    << " , "   <<   particles[j+3]   <<std::endl;
 	}
 }
-void exportVTK(std::ofstream& file, const int N, const FReal * particles ){
+//
+void exportCOSMOS(std::ofstream& file, const int N, const FReal * particles ){
 	int j = 0;
 	file << " x ,  y , z, q " <<std::endl;
 	for(int i = 0 ; i< N; ++i, j+=4){
-		file <<    particles[j]    << " , "    <<   particles[j+1]    << " , "   <<   particles[j+2]    << " , "   <<   particles[j+3]   <<std::endl;
+		file <<    particles[j]    << "  "    <<   particles[j+1]    << "  "   <<   particles[j+2]    << "  0.0 0.0 0.0  "   <<   particles[j+3]   <<"  " << i << std::endl;
 	}
-
 }
+//
+void exportVTK(std::ofstream& VTKfile, const int N, const FReal * particles ){
+	int j = 0;
+	//---------------------------
+	// print generic information
+	//---------------------------
+	VTKfile << "# vtk DataFile Version 3.0" << "\n";
+	VTKfile << "#  Generated bt exportVTK" << "\n";
 
+	VTKfile << "ASCII" << "\n";
+	VTKfile << "DATASET POLYDATA" << "\n";
+	//
+	//---------------------------------
+	// print nodes ordered by their TAG
+	//---------------------------------
+	VTKfile << "POINTS " << N << "  float" << "\n";
+	//
+	for(int i = 0 ; i< N; ++i, j+=4){
+		VTKfile <<    particles[j]    << "  "    <<   particles[j+1]    << "   "   <<   particles[j+2]      <<std::endl;
+	}
+	// ------------------------------------------
+	VTKfile << "\n";
+	VTKfile << "VERTICES  " <<  N << " " << 2*N << "\n";
+	for(int i = 0 ; i< N; ++i, j+=4){
+		VTKfile <<    "  1 "    << " "    <<i<<std::endl;
+	}
+	VTKfile << "POINT_DATA  " <<  N << "\n";
+	VTKfile << "SCALARS PhysicalValue  float 1" << "\n"
+			<< "LOOKUP_TABLE default" << "\n" ;
+	j = 0 ;
+	for(int i = 0 ; i< N; ++i, j+=4){
+		VTKfile <<    particles[j+3]    << " "    <<std::endl;
+	}
+	VTKfile << "\n";
+};
+
+void exportVTKxml(std::ofstream& VTKfile, const int N, const FReal * particles ){
+	int j = 0;
+
+	VTKfile << "<?xml version=\"1.0\"?>" <<std::endl
+			<< "<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"LittleEndian\"> "<<std::endl
+			<< "<PolyData>"<<std::endl
+			<< "<Piece NumberOfPoints=\" " << N << " \"  NumberOfVerts=\" "<<N <<" \" NumberOfLines=\" 0\" NumberOfStrips=\"0\" NumberOfPolys=\"0\">"<<std::endl
+			<< "<Points>"<<std::endl
+			<< "<DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\"> "<<std::endl ;
+	j = 0 ;
+	for(int i = 0 ; i< N; ++i, j+=4){
+		VTKfile <<    particles[j]    << "  "    <<   particles[j+1]    << "   "   <<   particles[j+2]      << "   "   ;
+	}
+	VTKfile <<std::endl<< "</DataArray> "<<std::endl
+			<< "</Points> "<<std::endl
+			<< "<PointData Scalars=\"PhysicalValue\" > "<<std::endl
+			<< "<DataArray type=\"Float64\" Name=\"PhysicalValue\"  format=\"ascii\">"<<std::endl ;
+	j = 0 ;
+	for(int i = 0 ; i< N; ++i, j+=4){
+		VTKfile <<    particles[j+3]    << " "   ;
+	}
+	VTKfile <<std::endl << "</DataArray>"<<std::endl
+			<< "	</PointData>"<<std::endl
+			<< "	<CellData>"<<" </CellData>"<<std::endl
+			<< "	<Verts>"<<std::endl
+			<< "	<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">"<<std::endl ;
+	for(int i = 0 ; i< N; ++i){
+		VTKfile <<   i   << " "   ;
+	}
+	VTKfile<<std::endl << "</DataArray>" <<std::endl
+			<< "<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">"<<std::endl ;
+	for(int i = 1 ; i< N+1; ++i){
+		VTKfile <<   i   << " "   ;
+	}
+	VTKfile<<std::endl  << "</DataArray>"<<std::endl
+			<< "	</Verts>"<<std::endl
+			<< "<Lines></Lines>"<<std::endl
+			<< "<Strips></Strips>"<<std::endl
+			<< "<Polys></Polys>"<<std::endl
+			<< "</Piece>"<<std::endl
+			<< "</PolyData>"<<std::endl
+			<< "</VTKFile>"<<std::endl;
+} ;
+//
 
 #endif
-- 
2.24.1