dolfin team mailing list archive

["Garth N. Wells" <g.n.wells@xxxxxxxxxxxxxxx>]

Attached are files for a Paraview/VTK interface (src/kernel/io/VtkFile.cpp and
src/kernel/io/dolfin/VtkFile.h). Also attached are patches for the necessary
modifications of files in src/kernel/io and src/kernel/io/dolfin. 

Results are output in the VTK XML format (file extension .vtu) and can be read
in by Paraview. Examples of postprocessed output from the demos "convdiff" and
"elasticity" are attached. The interface cannot handle time dependent data in
one file. A series of files can be made file001.vtu, file002.vtu, etc, as is
done in the elasticity demo. Paraview can easily add these files automtically to
produce animations.

Paraview produces bitmap Postscript files. The small program MayaVi
(http://mayavi.sourceforge.net/) can read the .vtu files and produces nice
vector Postscript.

Separate issue: the last backslash in /src/kernel/form/dolfin/Makefile.am should
be removed. It gives me problems when running automake.

Garth

// Copyright (C) 2005 Garth N. Wells.
// Licensed under the GNU GPL Version 2.
//

#include <dolfin/Mesh.h>
#include <dolfin/Function.h>
#include <dolfin/FiniteElement.h>
#include <dolfin/VtkFile.h>

using namespace dolfin;

//-­---------------------------------------------------------------------------
VtkFile::VtkFile(const std::string filename) : GenericFile(filename)
{
  type = "VTK";
}
//-­---------------------------------------------------------------------------
VtkFile::~VtkFile()
{
  // Do nothing
}
//-­---------------------------------------------------------------------------
void VtkFile::operator<<(Mesh& mesh)
{

  dolfin_info("Saving mesh to Vtk file.");
  
	// Write headers
  VtkHeaderOpen(mesh);

	MeshWrite(mesh);

  // Write headers
  VtkHeaderClose();

}
//-­---------------------------------------------------------------------------
void VtkFile::operator<<(Function& u)
{

  dolfin_info("Writing Function to Vtk file");
  

	const Mesh& mesh = u.mesh(); 

  // Write headers
  VtkHeaderOpen(mesh);

	// Write Mesh
  MeshWrite(mesh);

	// Write results
   ResultsWrite(u);

	// Close headers
	VtkHeaderClose();
	
  
  // Increase the number of times we have saved the function
  ++u;

  cout << "Saved function " << u.name() << " (" << u.label()
       << ") to file " << filename << " in VTK format." << endl;

}
//-­---------------------------------------------------------------------------
void VtkFile::MeshWrite(const Mesh& mesh) const
{

  // Open file
  FILE* fp = fopen(filename.c_str(), "a");

  // Write node positions
  fprintf(fp, "<Points>  \n");
  fprintf(fp, "<DataArray  type=\"Float32\"  NumberOfComponents=\"3\"  format=\"ascii\">  \n");
  for (NodeIterator n(mesh); !n.end(); ++n)
  {
    Point   p = n->coord();
    fprintf(fp," %f %f %f \n", p.x, p.y, p.z);
  }
  fprintf(fp, "</DataArray>  \n");
  fprintf(fp, "</Points>  \n");
  
  // Write cell connectivity
  fprintf(fp, "<Cells>  \n");
  fprintf(fp, "<DataArray  type=\"Int32\"  Name=\"connectivity\"  format=\"ascii\">  \n");
  for (CellIterator c(mesh); !c.end(); ++c)
  {
    for (NodeIterator n(c); !n.end(); ++n) fprintf(fp," %8d ",n->id());
    fprintf(fp," \n");
  }  
  fprintf(fp, "</DataArray> \n");

  // Write offset into connectivity array for the end of each cell
  fprintf(fp, "<DataArray  type=\"Int32\"  Name=\"offsets\"  format=\"ascii\">  \n");
  for (int offsets = 1; offsets <= mesh.noCells(); offsets++)
  {
    if (mesh.type() == Mesh::tetrahedra ) fprintf(fp, " %8d \n",  offsets*4);
    if (mesh.type() == Mesh::triangles )    fprintf(fp, " %8d \n", offsets*3);
  }
  fprintf(fp, "</DataArray> \n");
  
  //Write cell type
  fprintf(fp, "<DataArray  type=\"UInt8\"  Name=\"types\"  format=\"ascii\">  \n");
  for (int types = 1; types <= mesh.noCells(); types++)
  {
    if (mesh.type() == Mesh::tetrahedra ) fprintf(fp, " 10 \n");
    if (mesh.type() == Mesh::triangles )    fprintf(fp, " 5 \n");
  }
  fprintf(fp, "</DataArray> \n");
  fprintf(fp, "</Cells> \n"); 
  
  // Close file
  fclose(fp);

}
//-­---------------------------------------------------------------------------
void VtkFile::ResultsWrite(Function& u) const
{

  uint no_components = 0;
  const FiniteElement& element = u.element();

  if ( element.rank() == 0 )
  {
    no_components = 1;
  }
  else if ( element.rank() == 1 )
  {
    no_components = element.tensordim(0);
  }
  else
    dolfin_error("Cannot handle tensor valued functions.");

	// Open file
  FILE *fp = fopen(filename.c_str(), "a");

	  
	//Write PointData displacement	
	if(no_components == 1)
	{
		fprintf(fp, "<PointData  Scalars=\"U\"> \n");
		fprintf(fp, "<DataArray  type=\"Float32\"  Name=\"U\"  format=\"ascii\">	 \n");
  }
	else
	{
	fprintf(fp, "<PointData  Vectors=\"U\"> \n");
	fprintf(fp, "<DataArray  type=\"Float32\"  Name=\"U\"  NumberOfComponents=\"3\" format=\"ascii\">	 \n");	
	}
	
	for (NodeIterator n(u.mesh()); !n.end(); ++n)
  {    
		for(uint i =0; i < no_components; ++i)
		{
			fprintf(fp," %e ",u(*n, i));
  	}
			fprintf(fp,"\n");
	}	 
	fprintf(fp, "</DataArray> \n");
	fprintf(fp, "</PointData> \n");
	
    
  // Close file
  fclose(fp);

}
//-­---------------------------------------------------------------------------
void VtkFile::VtkHeaderOpen(const Mesh& mesh) const
{

	// Open file
  FILE *fp = fopen(filename.c_str(), "a");

  // Write headers
  fprintf(fp, "<VTKFile type=\"UnstructuredGrid\"  version=\"0.1\"  byte_order=\"LittleEndian\"  compressor=\"vtkZLibDataCompressor\">  \n");
  fprintf(fp, "<UnstructuredGrid>  \n");
  fprintf(fp, "<Piece  NumberOfPoints=\" %8d\"  NumberOfCells=\" %8d\">  \n", mesh.noNodes(), mesh.noCells());

  // Close file
  fclose(fp);
  
}
//-­---------------------------------------------------------------------------
void VtkFile::VtkHeaderClose() const
{

	// Open file
  FILE *fp = fopen(filename.c_str(), "a");

	// Close headers
  fprintf(fp, "</Piece> \n </UnstructuredGrid> \n </VTKFile>"); 	

  // Close file
  fclose(fp);
  

}
//-­---------------------------------------------------------------------------

// Copyright (C) 2005 Garth N. Wells.
// Licensed under the GNU GPL Version 2.
//

#ifndef __VTK_FILE_H
#define __VTK_FILE_H

#include <dolfin/GenericFile.h>

namespace dolfin {

  class VtkFile : public GenericFile {
  public:

    VtkFile(const std::string filename);
    ~VtkFile();

    void operator<< (Mesh& mesh);
    void operator<< (Function& u);

	private:
		void MeshWrite(const Mesh& mesh) const;
		void ResultsWrite(Function& u) const;
		void VtkHeaderOpen(const Mesh& mesh) const;
		void VtkHeaderClose() const;

  };
  
}

#endif

Attachment: patch_io_1
Description: Binary data

Attachment: patch_io_2
Description: Binary data

Attachment: ConvectionDiffusion.jpg
Description: JPEG image

Attachment: elasticity.avi
Description: MS Video