starcdreader.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_STARCD_READER_HH
#define DUNE_STARCD_READER_HH
 
#include <dune/common/exceptions.hh>
#include <dune/common/to_unique_ptr.hh>
 
#include <dune/geometry/type.hh>
#include <dune/grid/common/gridfactory.hh>
#include <iostream>
#include <fstream>
 
namespace Dune {
 
  template <class GridType>
  class StarCDReader {
 
  public:
 
    static ToUniquePtr<GridType> read(const std::string& fileName, bool verbose = true)
    {
      // extract the grid dimension
      const int dim = GridType::dimension;
 
      // currently only dim = 3 is implemented
      if (dim != 3)
        DUNE_THROW(Dune::NotImplemented,
                   "Reading Star-CD format is not implemented for dimension " << dim);
 
      // set up the grid factory
      GridFactory<GridType> factory;
 
      // set the name of the vertex file
      std::string vertexFileName = fileName + ".vrt";
 
      // set the vertex input stream
      std::ifstream vertexFile(vertexFileName.c_str());
      if (!vertexFile)
        DUNE_THROW(Dune::IOError, "Could not open " << vertexFileName);
 
      // read the vertices
      int dummyIdx;
      int numberOfVertices = 0;
      while (vertexFile >> dummyIdx) {
        numberOfVertices++;
 
        Dune::FieldVector<double,dim> position;
 
        for (int k = 0; k < dim; k++)
          vertexFile >> position[k];
 
        factory.insertVertex(position);
      }
      if (verbose)
        std::cout << numberOfVertices << " vertices read." << std::endl;
 
      // set the name of the element file
      std::string elementFileName = fileName + ".cel";
 
      // set the element input stream
      std::ifstream elementFile(elementFileName.c_str());
      if (!elementFile)
        DUNE_THROW(Dune::IOError, "Could not open " << elementFileName);
 
      // read the elements
      int numberOfElements = 0;
      int numberOfSimplices = 0;
      int numberOfPyramids = 0;
      int numberOfPrisms = 0;
      int numberOfCubes = 0;;
      int maxNumberOfVertices = (int)pow(2, dim);
      int isVolume = 1;
      while (elementFile >> dummyIdx) {
        std::vector<unsigned int> vertices(maxNumberOfVertices);
        for (int k = 0; k < maxNumberOfVertices; k++)
          elementFile >> vertices[k];
 
        int boundaryId;
        elementFile >> boundaryId;
 
        int volumeOrSurface[2];
        elementFile >> volumeOrSurface[0] >> volumeOrSurface[1];
 
        if (volumeOrSurface[0] == isVolume) {
          numberOfElements++;
 
          if (vertices[2] == vertices[3]) {           // simplex or prism
            if (vertices[4] == vertices[5]) {             // simplex
              numberOfSimplices++;
              std::vector<unsigned int> simplexVertices(4);
              for (int k = 0; k < 3; k++)
                simplexVertices[k] = vertices[k] - 1;
              simplexVertices[3] = vertices[4] - 1;
              factory.insertElement(Dune::GeometryTypes::tetrahedron, simplexVertices);
            }
            else {             // prism
              numberOfPrisms++;
              std::vector<unsigned int> prismVertices(6);
              for (int k = 0; k < 3; k++)
                prismVertices[k] = vertices[k] - 1;
              for (int k = 3; k < 6; k++)
                prismVertices[k] = vertices[k+1] - 1;
              factory.insertElement(Dune::GeometryTypes::prism, prismVertices);
            }
          }
          else {           // cube or pyramid
            if (vertices[4] == vertices[5]) {             // pyramid
              numberOfPyramids++;
              std::vector<unsigned int> pyramidVertices(5);
              for (int k = 0; k < 5; k++)
                pyramidVertices[k] = vertices[k] - 1;
              factory.insertElement(Dune::GeometryTypes::pyramid, pyramidVertices);
            }
            else {             // cube
              numberOfCubes++;
              std::vector<unsigned int> cubeVertices(8);
              for (int k = 0; k < 8; k++)
                cubeVertices[k] = vertices[k] - 1;
              std::swap(cubeVertices[2], cubeVertices[3]);
              std::swap(cubeVertices[6], cubeVertices[7]);
              factory.insertElement(Dune::GeometryTypes::hexahedron, cubeVertices);
            }
          }
        }
      }
      if (verbose)
        std::cout << numberOfElements << " elements read: "
                  << numberOfSimplices << " simplices, " << numberOfPyramids << " pyramids, "
                  << numberOfPrisms << " prisms, " << numberOfCubes << " cubes." << std::endl;
 
      // finish off the construction of the grid object
      if (verbose)
        std::cout << "Starting createGrid() ... " << std::flush;
 
      return factory.createGrid();
 
    }
 
  };
 
}
 
#endif