xmlfile.cpp

利用C

// Copyright (C) 2002-2008 Anders Logg.
// Licensed under the GNU LGPL Version 2.1.
//
// Modified by Erik Svensson 2003.
// Modified by Garth N. Wells 2006.
// Modified by Ola Skavhaug 2006.
// Modified by Magnus Vikstrom 2007.
//
// First added:  2002-12-03
// Last changed: 2008-04-22

#include <stdarg.h>

#include <dolfin/log/log.h>
#include <dolfin/common/constants.h>
#include <dolfin/common/Array.h>
#include <dolfin/la/GenericVector.h>
#include <dolfin/la/GenericMatrix.h>
#include <dolfin/mesh/Mesh.h>
#include <dolfin/mesh/Vertex.h>
#include <dolfin/mesh/Cell.h>
#include <dolfin/graph/Graph.h>
#include <dolfin/mesh/MeshFunction.h>
#include <dolfin/function/Function.h>
#include <dolfin/function/DiscreteFunction.h>
#include <dolfin/fem/DofMap.h>
#include <dolfin/parameter/Parameter.h>
#include <dolfin/parameter/ParameterList.h>

#include "XMLObject.h"
#include "XMLVector.h"
#include "XMLMatrix.h"
#include "XMLMesh.h"
#include "XMLMeshFunction.h"
#include "XMLDofMap.h"
#include "XMLFunction.h"
#include "XMLFiniteElement.h"
#include "XMLParameterList.h"
#include "XMLBLASFormData.h"
#include "XMLGraph.h"
#include "XMLFile.h"

using namespace dolfin;

//-----------------------------------------------------------------------------
XMLFile::XMLFile(const std::string filename) : GenericFile(filename),
					       header_written(false),
					       mark(0)
{
  type = "XML";
  xmlObject = 0;
}
//-----------------------------------------------------------------------------
XMLFile::~XMLFile()
{
  if ( xmlObject )
    delete xmlObject;
}
//-----------------------------------------------------------------------------
void XMLFile::operator>>(GenericVector& x)
{
  message(1, "Reading vector from file %s.", filename.c_str());

  if ( xmlObject )
    delete xmlObject;
  xmlObject = new XMLVector(x);
  parseFile();
}
//-----------------------------------------------------------------------------
void XMLFile::operator>>(GenericMatrix& A)
{
  message(1, "Reading matrix from file %s.", filename.c_str());

  if ( xmlObject )
    delete xmlObject;
  xmlObject = new XMLMatrix(A);
  parseFile();
}
//-----------------------------------------------------------------------------
void XMLFile::operator>>(Mesh& mesh)
{
  message(1, "Reading mesh from file %s.", filename.c_str());

  if ( xmlObject )
    delete xmlObject;
  xmlObject = new XMLMesh(mesh);
  parseFile();
}
//-----------------------------------------------------------------------------
void XMLFile::operator>>(MeshFunction<int>& meshfunction)
{
  message(1, "Reading int-valued mesh function from file %s.", filename.c_str());

  if ( xmlObject )
    delete xmlObject;
  xmlObject = new XMLMeshFunction(meshfunction);
  parseFile();
}
//-----------------------------------------------------------------------------
void XMLFile::operator>>(MeshFunction<unsigned int>& meshfunction)
{
  message(1, "Reading uint-valued mesh function from file %s.", filename.c_str());

  if ( xmlObject )
    delete xmlObject;
  xmlObject = new XMLMeshFunction(meshfunction);
  parseFile();
}
//-----------------------------------------------------------------------------
void XMLFile::operator>>(MeshFunction<double>& meshfunction)
{
  message(1, "Reading real-valued mesh function from file %s.", filename.c_str());

  if ( xmlObject )
    delete xmlObject;
  xmlObject = new XMLMeshFunction(meshfunction);
  parseFile();
}
//-----------------------------------------------------------------------------
void XMLFile::operator>>(MeshFunction<bool>& meshfunction)
{
  message(1, "Reading bool-valued mesh function from file %s.", filename.c_str());

  if ( xmlObject )
    delete xmlObject;
  xmlObject = new XMLMeshFunction(meshfunction);
  parseFile();
}
//-----------------------------------------------------------------------------
void XMLFile::operator>>(Function& f)
{
  // We are cheating here. Instead of actually parsing the XML for
  // Function data nested inside <function></function>, we just ignore
  // the nesting and look for the first occurence of the data which
  // might be outide of <function></function>

  message(1, "Reading function from %s.", filename.c_str());

  // Read the vector
  Vector* x = new Vector();
  *this >> *x;

  // Read the mesh
  Mesh* mesh = new Mesh();
  *this >> *mesh;

  // Read the finite element specification
  std::string finite_element_signature;
  if ( xmlObject )
    delete xmlObject;
  xmlObject = new XMLFiniteElement(finite_element_signature);
  parseFile(); 

  // Read the dof map specification
  std::string dof_map_signature;
  if ( xmlObject )
    delete xmlObject;
  xmlObject = new XMLDofMap(dof_map_signature);
  parseFile(); 

  // Read the function
  if ( xmlObject )
    delete xmlObject;
  xmlObject = new XMLFunction(f);
  parseFile(); 
  
  // Create the function (we're all friends here)
  if ( f.f )
    delete f.f;
  f.f = new DiscreteFunction(*mesh, *x, finite_element_signature, dof_map_signature);
  f._type = Function::discrete;
  f.rename("u", "discrete function from file data");
}
//-----------------------------------------------------------------------------
void XMLFile::operator>>(ParameterList& parameters)
{
  message(1, "Reading parameter list from file %s.", filename.c_str());

  if ( xmlObject )
    delete xmlObject;
  xmlObject = new XMLParameterList(parameters);
  parseFile();
}
//-----------------------------------------------------------------------------
void XMLFile::operator>>(BLASFormData& blas)
{
  if ( xmlObject )
    delete xmlObject;
  xmlObject = new XMLBLASFormData(blas);
  parseFile();
}
//-----------------------------------------------------------------------------
void XMLFile::operator>>(Graph& graph)
{
  message(1, "Reading graph from file %s.", filename.c_str());

  if ( xmlObject )
    delete xmlObject;
  xmlObject = new XMLGraph(graph);
  parseFile();
}
//-----------------------------------------------------------------------------
void XMLFile::operator<<(GenericVector& x)
{
  // Open file
  FILE* fp = openFile();

  // Get vector values
  real* values = new real[x.size()];
  x.get(values);
  
  // Write vector in XML format
  fprintf(fp, "  <vector size=\"%u\"> \n", x.size() );
  for (unsigned int i = 0; i < x.size(); i++) 
  {
    fprintf(fp, "    <entry row=\"%u\" value=\"%.15g\"/>\n", i, values[i]);
    if ( i == (x.size() - 1))
      fprintf(fp, "  </vector>\n");
  }
  
  // Delete vector values
  delete [] values;

  // Close file
  closeFile(fp);
  
//  message(1, "Saved vector %s (%s) to file %s in DOLFIN XML format.", x.name().c_str(), x.label().c_str(), filename.c_str());
  message(1, "Saved vector  to file %s in DOLFIN XML format.", filename.c_str());
}
//-----------------------------------------------------------------------------
void XMLFile::operator<<(GenericMatrix& A)
{
  // Open file
  FILE *fp = openFile();
  
  // Write matrix in XML format
  fprintf(fp, "  <matrix rows=\"%u\" columns=\"%u\">\n", A.size(0), A.size(1));
        
  Array<uint> columns;
  Array<real> values;

  for (unsigned int i = 0; i < A.size(0); i++)
  {
    A.getrow(i, columns, values);
    if (columns.size() > 0)
      fprintf(fp, "    <row row=\"%u\" size=\"%d\">\n", i, (int)columns.size());
    for (uint pos = 0; pos < columns.size(); pos++)
    {
      unsigned int j = columns[pos];
      real aij = values[pos];
      fprintf(fp, "      <entry column=\"%u\" value=\"%.15g\"/>\n", j, aij);
    }
    if (columns.size() > 0 )
      fprintf(fp, "    </row>\n");
  }
  fprintf(fp, "  </matrix>\n");

  // Close file
  closeFile(fp);

  message(1, "Saved matrix file %s in DOLFIN XML format.", filename.c_str());
}
//-----------------------------------------------------------------------------
void XMLFile::operator<<(Mesh& mesh)
{
  // Open file
  FILE *fp = openFile();
  
  // Get cell type
  CellType::Type cell_type = mesh.type().cellType();

  // Write mesh in XML format
  fprintf(fp, "  <mesh celltype=\"%s\" dim=\"%u\">\n",
          CellType::type2string(cell_type).c_str(), mesh.geometry().dim());

  fprintf(fp, "    <vertices size=\"%u\">\n", mesh.numVertices());
  
  for(VertexIterator v(mesh); !v.end(); ++v)
  {
    Point p = v->point();

    switch ( mesh.geometry().dim() ) {
    case 1:
      fprintf(fp, "      <vertex index=\"%u\" x=\"%g\"/>\n",
              v->index(), p.x());
      break;
    case 2:
      fprintf(fp, "      <vertex index=\"%u\" x=\"%g\" y=\"%g\"/>\n",
              v->index(), p.x(), p.y());
      break;
    case 3:
      fprintf(fp, "      <vertex index=\"%u\" x=\"%g\" y=\"%g\" z=\"%g\" />\n",
              v->index(), p.x(), p.y(), p.z());
      break;
    default:
      error("The XML mesh file format only supports 1D, 2D and 3D meshes.");
    }
  }

  fprintf(fp, "    </vertices>\n");
  fprintf(fp, "    <cells size=\"%u\">\n", mesh.numCells());

  for (CellIterator c(mesh); !c.end(); ++c)
  {
    uint* vertices = c->entities(0);
    dolfin_assert(vertices);

    switch ( cell_type )
    {
    case CellType::interval:
      fprintf(fp, "      <interval index=\"%u\" v0=\"%u\" v1=\"%u\"/>\n",
	      c->index(), vertices[0], vertices[1]);
      break;
    case CellType::triangle:
      fprintf(fp, "      <triangle index=\"%u\" v0=\"%u\" v1=\"%u\" v2=\"%u\"/>\n",
	      c->index(), vertices[0], vertices[1], vertices[2]);
      break;
    case CellType::tetrahedron:
      fprintf(fp, "      <tetrahedron index=\"%u\" v0=\"%u\" v1=\"%u\" v2=\"%u\" v3=\"%u\"/>\n",
              c->index(), vertices[0], vertices[1], vertices[2], vertices[3]);
      break;
    default:
      error("Unknown cell type: %u.", cell_type);
    }
  }

  fprintf(fp, "    </cells>\n");
  fprintf(fp, "  </mesh>\n");
 
  // Close file
  closeFile(fp);

  message(1, "Saved mesh to file %s in DOLFIN XML format.", filename.c_str());
}
//-----------------------------------------------------------------------------
void XMLFile::operator<<(MeshFunction<int>& meshfunction)
{
  // Open file
  FILE *fp = openFile();
  
  // Write mesh in XML format
  fprintf(fp, "  <meshfunction type=\"int\" dim=\"%u\" size=\"%u\">\n",
          meshfunction.dim(), meshfunction.size());
  
  Mesh& mesh = meshfunction.mesh();
  for(MeshEntityIterator e(mesh, meshfunction.dim()); !e.end(); ++e)
  {
      fprintf(fp, "    <entity index=\"%u\" value=\"%d\"/>\n",
              e->index(), meshfunction(*e));
  }