Index: /issm/trunk-jpl/src/m/contrib/paraview/writeVTKcell.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/paraview/writeVTKcell.m	(revision 15010)
+++ /issm/trunk-jpl/src/m/contrib/paraview/writeVTKcell.m	(revision 15010)
@@ -0,0 +1,119 @@
+function writeVTKcell(filename,model,Solution)
+% vtk export
+% creates a vtk-file filename.vtk containing simplicial mesh data
+% (only work for triangle now)
+% and additional cell data
+%
+% input: filename   destination 
+%                   (string)
+%------------------------------------------------------------------
+%        model      this is md 
+%------------------------------------------------------------------
+%        u          mesh function assigning a real number to every
+%                   element
+%                   of the triangulation/tetrahedralization
+%                   (Mx1 array)
+% Basile de Fleurian, modified from:
+% (c) Daniel Peterseim, 2009-11-07
+
+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(filename,'.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);
+	s='%f %f %f \n';
+	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(filename,'.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);
+		s='%f %f %f \n';
+		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
Index: /issm/trunk-jpl/src/m/contrib/paraview/writeVTKcell.m.bck
===================================================================
--- /issm/trunk-jpl/src/m/contrib/paraview/writeVTKcell.m.bck	(revision 15010)
+++ /issm/trunk-jpl/src/m/contrib/paraview/writeVTKcell.m.bck	(revision 15010)
@@ -0,0 +1,65 @@
+function writeVTKcell(filename,t,p,u)
+% vtk export
+% creates a vtk-file filename.vtk containing simplicial mesh data
+% (2- or 3d)
+% and additional cell data
+%
+% input: filename   destination 
+%                   (string)
+%        p          array of N points 
+%                   (Nxd matrix where d denotes the dimension)
+%        t          triangulation/tetrahedralization of the points
+%        in p
+%                   (Mxd+1 array, where M denotes the number of
+%                   simplices)
+%        u          mesh function assigning a real number to every
+%        element
+%                   of the triangulation/tetrahedralization
+%                   (Mx1 array)
+%
+% example usage:
+%        2d: p=rand(10,2); 
+%            t=delaunayn(p); 
+%            u=sum(t,2);
+%            writeVTKcell('test2d',t,p,u);
+%        3d: p=rand(10,3); 
+%            t=delaunayn(p); 
+%            u=sum(t,2);
+%            writeVTKcell('test3d',t,p,u);
+% (the result is accessible with paraview!)
+%
+% (c) Daniel Peterseim, 2009-11-07
+
+[np,dim]=size(p);
+[nt]=size(t,1);
+celltype=5; %gives the type of element, 5 for triangles
+
+FID = fopen(strcat(filename,'.vtk'),'w+');
+fprintf(FID,'# vtk DataFile Version 2.0\n');
+fprintf(FID,'Unstructured Grid Example\n');
+fprintf(FID,'ASCII\n');
+fprintf(FID,'DATASET UNSTRUCTURED_GRID\n');
+
+fprintf(FID,'POINTS %d float\n',np);
+s='%f %f %f \n';
+P=[p zeros(np,3-dim)];
+fprintf(FID,s,P');
+
+fprintf(FID,'CELLS %d %d\n',nt,nt*(dim+1));
+s='%d';
+for k=1:dim
+    s=horzcat(s,{' %d'});
+end
+s=cell2mat(horzcat(s,{'\n'}));
+fprintf(FID,s,[(dim)*ones(nt,1) t-1]');
+
+fprintf(FID,'CELL_TYPES %d\n',nt);
+s='%d\n';
+fprintf(FID,s,celltype*ones(nt,1));
+
+fprintf(FID,'POINT_DATA %s\nSCALARS Hwater float 1\n',num2str(np));
+fprintf(FID,'LOOKUP_TABLE default\n');
+s='%f\n';
+fprintf(FID,s,u);
+
+fclose(FID);