dolfin team mailing list archive

[Minh Do Quang <minh@xxxxxxxxxxx>]

Hi,

I just tried to wake up the function OpenDXFile. It now can handle the
discrete functions as in VTK. Each time it will write to a separate
files instead of writing all together as before.

/ Minh

// Copyright (C) 2003-2006 Johan Hoffman and Anders Logg.
// Licensed under the GNU LGPL Version 2.1.
//
// Modified by Minh Do-Quang 2006-2008
//
// First added:  2003-07-15
// Last changed: 2008-02-15

#include <stdio.h>
#include <dolfin/log.h>
#include <dolfin/parameters.h>
#include <dolfin/Cell.h>
#include <dolfin/Function.h>
#include <dolfin/Mesh.h>
#include <dolfin/Vertex.h>
#include <dolfin/OpenDXFile.h>

using namespace dolfin;

//-�---------------------------------------------------------------------------
OpenDXFile::OpenDXFile(const std::string filename) : GenericFile(filename)
{
  type = "OpenDX";
}
//-�---------------------------------------------------------------------------
OpenDXFile::~OpenDXFile()
{
  // Do nothing
}
//-�---------------------------------------------------------------------------
void OpenDXFile::operator<<(Mesh& mesh)
{
  // Update dx file name and clear file
  dxNameUpdate(counter);

  // Open file
  FILE* fp = fopen(dx_filename.c_str(), "w");
  
  // Write header first time
  writeHeader(fp);

  // Write mesh
  writeMesh(fp, mesh);

  // Close headers
  writeMeshData(fp, mesh);
  
  // Increase the number of times we have saved the mesh
  counter++;

  message(1, "Saved mesh %s (%s) to file %s in OpenDX format.",
          mesh.name().c_str(), mesh.label().c_str(), filename.c_str());

  // Close file
  fclose(fp);
}
//-�---------------------------------------------------------------------------
void OpenDXFile::operator<<(Function& u)
{
  // Can only save discrete functions
  if ( u.type() != Function::discrete )
    error("Only discrete functions can be saved in OpenDX format.");

  // Update dx file name and clear file
  dxNameUpdate(counter);

  FILE* fp = fopen(dx_filename.c_str(), "w");
  fseek(fp, 0L, SEEK_END);

  // Write header first time
  if ( ftell(fp) == 0 )
    writeHeader(fp);

  // Write mesh
  if ( save_each_mesh )
    writeMesh(fp, u.mesh());

  // Write function
  writeFunction(fp, u);

  // Close headers
  writeHeaderClose(fp);
  
  // Close file
  fclose(fp);

  // Increase the number of times we have saved the function
  counter++;
  
  cout << "Saved function " << u.name() << " (" << u.label()
       << ") to file " << filename << " in OpenDX format." << endl;
}
//-�---------------------------------------------------------------------------
void OpenDXFile::writeHeader(FILE* fp)
{
  fprintf(fp,"# Output from DOLFIN version %s.\n", DOLFIN_VERSION);
  fprintf(fp,"# Format intended for use with OpenDX (Data Explorer).\n");
  fprintf(fp,"#\n");
  fprintf(fp,"\n");
}
//-�---------------------------------------------------------------------------
void OpenDXFile::writeHeaderClose(FILE* fp)
{
  fprintf(fp,"#\n");
  fprintf(fp,"#\n");
  fprintf(fp,"object \"maindata\" class field\n");
  fprintf(fp,"component \"positions\"    \"vertex\"\n");
  fprintf(fp,"component \"connections\"  \"cells\"\n");
  fprintf(fp,"component \"data\"         \"data\"\n");
  fprintf(fp,"end\n");
  
  fprintf(fp,"\n");
}
//-�---------------------------------------------------------------------------
void OpenDXFile::writeMesh(FILE* fp, Mesh& mesh)
{
  int meshdim, celldim;
  if( mesh.type().cellType() == CellType::tetrahedron ) {
    meshdim = 3;
    celldim = 4;
  } else {
    meshdim = 2;
    celldim = 3;
  }

  // Write vertices
  fprintf(fp, "# A list of all vertex positions\n");
  fprintf(fp, "object \"vertex\" class array type float rank 1 shape %d items %u data follows\n", 
	  meshdim, mesh.numVertices());

  for (VertexIterator n(mesh); !n.end(); ++n)
  {
    Point p = n->point();
    
    float x = (float) p.x();
    float y = (float) p.y();
    float z = (float) p.z();
    
    if ( mesh.type().cellType() == CellType::tetrahedron )
      fprintf(fp, "%g %g %g\n", p.x(), p.y(), p.z());
    else if ( mesh.type().cellType() == CellType::triangle )
      fprintf(fp, "%g %g\n", p.x(), p.y());
      
  }
  fprintf(fp,"#\n");

  // Write cells
  fprintf(fp, "# A list of all cells (connections)\n");
  fprintf(fp, "object \"cells\" class array type int rank 1 shape %d items %u data follows\n",
	  celldim, mesh.numCells());

  for (CellIterator c(mesh); !c.end(); ++c)
  {  
    for (VertexIterator n(*c); !n.end(); ++n)
      fprintf(fp, " %8u ", n->index());
    fprintf(fp, "\n");
  }
  fprintf(fp, "#\n");
  if ( mesh.type().cellType() == CellType::tetrahedron )
    fprintf(fp, "attribute \"element type\" string \"tetrahedra\"\n");
  else if ( mesh.type().cellType() == CellType::triangle )
    fprintf(fp, "attribute \"element type\" string \"triangles\"\n");
  fprintf(fp, "attribute \"ref\" string \"positions\"\n");
  fprintf(fp, "\n");
}
//-�---------------------------------------------------------------------------
void OpenDXFile::writeMeshData(FILE* fp, Mesh& mesh)
{
  // Write data for a given mesh to create an object that can be visualized
  // when no data is associated with the mesh. This is necessary when we only
  // want to save the mesh.
  
  // Check that we don't try to write mesh data at the same time as we
  // write a time series

  fprintf(fp,"# Cell diameter\n");
  fprintf(fp,"object \"diameter\" class array type float rank 0 items %u data follows\n",
	  mesh.numCells());
  
  for (CellIterator c(mesh); !c.end(); ++c)
    {
      float value = static_cast<float>(c->diameter());
      fprintf(fp, "%g\n", value);
    }
  fprintf(fp, "\n");
  fprintf(fp, "attribute \"dep\" string \"connections\"\n");
  fprintf(fp, "\n");  
  
  // Write the mesh
  fprintf(fp, "# The mesh\n");
  fprintf(fp, "object \"Mesh\" class field\n");
  fprintf(fp, "component \"positions\" value \"vertex\"\n");
  fprintf(fp, "component \"connections\" value \"cells\"\n");
  fprintf(fp, "component \"data\" value \"diameter\"\n");
}

//-�---------------------------------------------------------------------------
void OpenDXFile::writeFunction(FILE* fp, Function& u)
{
  Mesh& mesh = u.mesh();

  const uint rank = u.rank();
  if(rank > 1)
    error("Only scalar and vectors functions can be saved in VTK format.");

  // Get number of components
  const uint dim = u.dim(0);

  // Allocate memory for function values at vertices
  uint size = mesh.numVertices();
  for (uint i = 0; i < u.rank(); i++)
    size *= u.dim(i);
  real* values = new real[size];

  // Write header for object
  fprintf(fp,"# Values for [%s] at nodal points\n", u.label().c_str());
  fprintf(fp,"object \"data\" class array type float rank 1 shape %u items %u data follows\n",
	  dim, mesh.numVertices());
  
  // Get function values at vertices
  u.interpolate(values);

  // Write function data at mesh vertices
  if ( rank <=1 )
  {
    for (VertexIterator vertex(mesh); !vertex.end(); ++vertex)
      {
	for(unsigned int i=0; i<dim; i++) 
	  fprintf(fp, " %e ", values[vertex->index()+i*mesh.numVertices()]);
	fprintf(fp, "\n");
      }
  }
  fprintf(fp,"\n");
  fprintf(fp,"attribute \"dep\" string \"positions\"\n\n");
  
}
//----------------------------------------------------------------------------
void OpenDXFile::dxNameUpdate(const int counter) 
{
  std::string filestart, extension;
  std::ostringstream fileid, newfilename;
  
  fileid.fill('0');
  fileid.width(6);
  
  filestart.assign(filename, 0, filename.find("."));
  extension.assign(filename, filename.find("."), filename.size());
  
  fileid << counter;
  newfilename << filestart << fileid.str() << ".dx";
  
  dx_filename = newfilename.str();
  
}

// Copyright (C) 2003-2005 Johan Hoffman and Anders Logg.
// Licensed under the GNU LGPL Version 2.1.
//
// Modified by Minh Do-Quang 2008.
//
// First added:  2003-07-15
// Last changed: 2008

#ifndef __OPEN_DX_FILE_H
#define __OPEN_DX_FILE_H

#include <string>
#include <stdio.h>
#include <dolfin/constants.h>
#include <dolfin/Event.h>
#include <dolfin/Array.h>
#include <dolfin/GenericFile.h>

namespace dolfin
{
    
  class OpenDXFile : public GenericFile
  {
  public:
    
    OpenDXFile(const std::string filename);
    ~OpenDXFile();
    
    void operator<< (Mesh& mesh);
    void operator<< (Function& u);

  private:

    void writeHeader   (FILE* fp);
    void writeHeaderClose(FILE* fp);
    void writeMesh     (FILE* fp, Mesh& mesh);
    void writeMeshData (FILE* fp, Mesh& mesh);
    void writeFunction (FILE* fp, Function& u);
    void dxNameUpdate(const int counter);  

    // dx filename
    std::string dx_filename;

    // Check if we should save each mesh
    bool save_each_mesh;

  };
  
}

#endif