function exportVTK(filename,model,Solution)
% vtk export
% function exportVTK(filename,model,Solution)
% creates a directory with the vtk files for you simulation
% (only work for triangle and wedges based on their number of nodes)
% and additional cell data
%
% input: filename   destination 
%                   (string)
%------------------------------------------------------------------
%        model      this is md 
%------------------------------------------------------------------
%        Solution   Put on the enum of the solution you want 
%	                  to plot         
% Basile de Fleurian:

[path,name,ext]=fileparts(filename);
separator=filesep;
mkdir(filename);

points=[model.mesh.x model.mesh.y model.mesh.z];
[num_of_points,dim]=size(points);
[num_of_elt]=size(model.mesh.elements,1);
[point_per_elt]=size(model.mesh.elements,2);
sol_enum=EnumToString(Solution);

sol_struct=model.results.(sol_enum);

%Select the type of element function of the number of nodes per elements
if point_per_elt==3;
	celltype=5; %triangles
elseif point_per_elt==6;
	celltype=13; %wedges
else
	error('Your Element definition is not taken into account \n');
end

%looking for multiple time steps
num_of_timesteps=size(sol_struct,2);

%getting the number of fields in the solution
fieldnames=fields(sol_struct(1));
num_of_fields=length(fieldnames);

if num_of_timesteps==1; %just one timestep only write one file

	FID = fopen(strcat(path,filesep,name,filesep,name,'.vtk'),'w+');
	fprintf(FID,'# vtk DataFile Version 2.0 \n');
	fprintf(FID,'Data for run %s \n',model.miscellaneous.name);
	fprintf(FID,'ASCII \n');
	fprintf(FID,'DATASET UNSTRUCTURED_GRID \n');

	fprintf(FID,'POINTS %d float\n',num_of_points);
	if(dim==3);
		s='%f %f %f \n';
	elseif(dim==2);
		s='%f %f \n';
  end
	P=[points zeros(num_of_points,3-dim)];
	fprintf(FID,s,P');

	fprintf(FID,'CELLS %d %d\n',num_of_elt,num_of_elt*(point_per_elt+1));
	s='%d';
	for k=1:point_per_elt
    s=horzcat(s,{' %d'});
  end
	s=cell2mat(horzcat(s,{'\n'}));
	fprintf(FID,s,[(point_per_elt)*ones(num_of_elt,1) model.mesh.elements-1]');

	fprintf(FID,'CELL_TYPES %d\n',num_of_elt);
	s='%d\n';
	fprintf(FID,s,celltype*ones(num_of_elt,1));

	%check which field is a real result and print	
	fprintf(FID,'POINT_DATA %s \n',num2str(num_of_points));
	for j=1:num_of_fields

		if (length(sol_struct(1).(fieldnames{j}))==num_of_points);
			fprintf(FID,'SCALARS %s float 1 \n',fieldnames{j});
			fprintf(FID,'LOOKUP_TABLE default\n');
			s='%e\n';
			fprintf(FID,s,sol_struct.(fieldnames{j}));
	  end
  end
		fclose(FID);
else
	for i=1:num_of_timesteps;
		timestep=sol_struct(i).step;
		FID = fopen(strcat(path,filesep,name,filesep,name,'.vtk',int2str(timestep),'.vtk'),'w+');
		fprintf(FID,'# vtk DataFile Version 2.0 \n');
		fprintf(FID,'Data for run %s \n',model.miscellaneous.name);
		fprintf(FID,'ASCII \n');
		fprintf(FID,'DATASET UNSTRUCTURED_GRID \n');

		fprintf(FID,'POINTS %d float\n',num_of_points);
		if(dim==3);
			s='%f %f %f \n';
		elseif(dim==2);
			s='%f %f \n';
   end
		P=[points zeros(num_of_points,3-dim)];
		fprintf(FID,s,P');

		fprintf(FID,'CELLS %d %d\n',num_of_elt,num_of_elt*(point_per_elt+1));
		s='%d';
		for k=1:point_per_elt
			s=horzcat(s,{' %d'});
  	end
		s=cell2mat(horzcat(s,{'\n'}));
		fprintf(FID,s,[(point_per_elt)*ones(num_of_elt,1) model.mesh.elements-1]');

		fprintf(FID,'CELL_TYPES %d\n',num_of_elt);
		s='%d\n';
		fprintf(FID,s,celltype*ones(num_of_elt,1));

		%check which field is a real result and print
		fprintf(FID,'POINT_DATA %s \n',num2str(num_of_points));
		for j=1:num_of_fields

			if (length(sol_struct(1).(fieldnames{j}))==num_of_points);
				fprintf(FID,'SCALARS %s float 1 \n',fieldnames{j});
				fprintf(FID,'LOOKUP_TABLE default\n');
				s='%e\n';
				fprintf(FID,s,sol_struct(i).(fieldnames{j}));
	    end		
    end 
		fclose(FID);
  end
end
