gmshwriter.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_GRID_IO_FILE_GMSHWRITER_HH
#define DUNE_GRID_IO_FILE_GMSHWRITER_HH
 
 
#include <fstream>
#include <iostream>
 
#include <string>
 
#include <dune/common/exceptions.hh>
 
#include <dune/geometry/type.hh>
 
#include <dune/grid/common/grid.hh>
#include <dune/grid/common/mcmgmapper.hh>
 
 
namespace Dune {
 
  template <class GridView>
  class GmshWriter
  {
  private:
    const GridView gv;
 
    static const int dim = GridView::dimension;
    static const int dimWorld = GridView::dimensionworld;
    static_assert( (dimWorld <= 3), "GmshWriter requires dimWorld <= 3." );
 
    typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
 
 
    size_t nodeIndexFromIterator(const ElementIterator& eIt, int i) const {
      return gv.indexSet().subIndex(*eIt, i, dim)+1;
    }
 
    static size_t translateDuneToGmshType(const GeometryType& type) {
      // Probably the non-clever, but hopefully readable way of translating the GeometryTypes to the gmsh types
      size_t element_type;
 
      if (type.isLine())
        element_type = 1;
      else if (type.isTriangle())
        element_type = 2;
      else if (type.isQuadrilateral())
        element_type = 3;
      else if (type.isTetrahedron())
        element_type = 4;
      else if (type.isHexahedron())
        element_type = 5;
      else if (type.isPrism())
        element_type = 6;
      else if (type.isPyramid())
        element_type = 7;
      else if (type.isVertex())
        element_type = 15;
      else
        DUNE_THROW(Dune::IOError, "GeometryType " << type << " is not supported by gmsh.");
 
      return element_type;
    }
 
    void outputElements(std::ofstream& file,
                        const std::vector<int>* physicalEntities = nullptr) const {
      MultipleCodimMultipleGeomTypeMapper<GridView, MCMGElementLayout> elementMapper(gv);
      ElementIterator eIt    = gv.template begin<0>();
      ElementIterator eEndIt = gv.template end<0>();
 
      for (size_t i = 1; eIt != eEndIt; ++eIt, ++i) {
        // Check whether the type is compatible. If not, close file and rethrow exception.
        try {
          size_t element_type = translateDuneToGmshType(eIt->type());
 
          file << i << " " << element_type;
          // If present, set the first tag to the physical entity
          if (physicalEntities) {
            file << " " << 1 << " " << (*physicalEntities)[elementMapper.index(*eIt)];
          } else {
            file << " " << 0; // "0" for "I do not use any tags."
          }
 
          // Output list of nodes.
          // 3, 5 and 7 got different vertex numbering compared to Dune
          if (3 == element_type)
            file << " "
                 << nodeIndexFromIterator(eIt, 0) << " " << nodeIndexFromIterator(eIt, 1) << " "
                 << nodeIndexFromIterator(eIt, 3) << " " << nodeIndexFromIterator(eIt, 2);
          else if (5 == element_type)
            file << " "
                 << nodeIndexFromIterator(eIt, 0) << " " << nodeIndexFromIterator(eIt, 1) << " "
                 << nodeIndexFromIterator(eIt, 3) << " " << nodeIndexFromIterator(eIt, 2) << " "
                 << nodeIndexFromIterator(eIt, 4) << " " << nodeIndexFromIterator(eIt, 5) << " "
                 << nodeIndexFromIterator(eIt, 7) << " " << nodeIndexFromIterator(eIt, 6);
          else if (7 == element_type)
            file << " "
                 << nodeIndexFromIterator(eIt, 0) << " " << nodeIndexFromIterator(eIt, 1) << " "
                 << nodeIndexFromIterator(eIt, 3) << " " << nodeIndexFromIterator(eIt, 2) << " "
                 << nodeIndexFromIterator(eIt, 4);
          else {
            for (int k = 0; k < eIt->geometry().corners(); ++k)
              file << " " << nodeIndexFromIterator(eIt, k);
          }
 
          file << std::endl;
 
        } catch(Exception& e) {
          file.close();
          throw;
        }
      }
    }
 
 
    void outputNodes(std::ofstream& file) const {
      VertexIterator vIt    = gv.template begin<dim>();
      VertexIterator vEndIt = gv.template end<dim>();
 
      for (; vIt != vEndIt; ++vIt) {
        typename VertexIterator::Entity::Geometry::GlobalCoordinate globalCoord = vIt->geometry().center();
        int nodeIndex = gv.indexSet().index(*vIt)+1; // Start counting indices by "1".
 
        if (1 == dimWorld)
          file << nodeIndex << " " << globalCoord[0] << " " << 0 << " " << 0 << std::endl;
        else if (2 == dimWorld)
          file << nodeIndex << " " << globalCoord[0] << " " << globalCoord[1] << " " << 0 << std::endl;
        else // (3 == dimWorld)
          file << nodeIndex << " " << globalCoord[0] << " " << globalCoord[1] << " " << globalCoord[2] << std::endl;
      }
    }
 
    void writeImpl(const std::string& fileName, const std::vector<int>* physicalEntities = nullptr) const {
      std::ofstream file(fileName.c_str());
 
      if (!file.is_open())
        DUNE_THROW(Dune::IOError, "Could not open " << fileName << " with write access.");
 
      // Output Header
      file << "$MeshFormat" << std::endl
           << "2.0 0 " << sizeof(double) << std::endl // "2.0" for "version 2.0", "0" for ASCII
           << "$EndMeshFormat" << std::endl;
 
 
      // Output Nodes
      const size_t number_of_nodes = gv.size(dim);
      file << "$Nodes" << std::endl
           << number_of_nodes << std::endl;
 
      outputNodes(file);
 
      file << "$EndNodes" << std::endl;
 
 
      // Output Elements
      const size_t number_of_elements = gv.size(0);
      file << "$Elements" << std::endl
           << number_of_elements << std::endl;
 
      outputElements(file, physicalEntities);
 
      file << "$EndElements" << std::endl; // Add an additional new line for good measure.
 
 
      file.close();
    }
 
 
  public:
    GmshWriter(const GridView& gridView) : gv(gridView) {}
 
    void write(const std::string& fileName) const {
      writeImpl(fileName);
    }
 
    void write(const std::string& fileName, const std::vector<int>& physicalEntities) const {
      writeImpl(fileName, &physicalEntities);
    }
  };
 
} // namespace Dune
 
#endif // DUNE_GRID_IO_FILE_GMSHWRITER_HH