gmshgridfactory.hh

 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_MMESH_GRID_GMSHGRIDFACTORY_HH
#define DUNE_MMESH_GRID_GMSHGRIDFACTORY_HH
 
// System includes
#include <memory>
#include <utility>
#include <vector>
 
// Dune includes
#include <dune/grid/common/intersection.hh>
 
// MMesh includes
#include <dune/mmesh/grid/explicitgridfactory.hh>
#include <dune/mmesh/grid/implicitgridfactory.hh>
 
namespace Dune
{
 
  // External Forward Declarations
  // -----------------------------
  template< class GridImp, class IntersectionImp >
  class Intersection;
 
  // GmshGridFactory for MMesh
  // ------------------------
  template< class Grid, bool useImplicitGridFactory = false >
  struct GmshGridFactory
  {
    const static int dimension = Grid::dimension;
    typedef MPIHelper::MPICommunicator MPICommunicatorType;
    typedef typename Grid::template Codim<0>::Entity Element;
    typedef typename Grid::template Codim<dimension>::Entity Vertex;
    typedef typename Grid::FieldType FieldType;
    typedef FieldVector< FieldType, dimension > GlobalPosition;
 
    // Choose explicit or implicit grid factory
    typedef std::conditional_t< useImplicitGridFactory,
      MMeshImplicitGridFactory<Grid>,
      MMeshExplicitGridFactory<Grid>
    > GridFactory;
 
    explicit GmshGridFactory ( std::istream &input )
    {
      input.clear();
      input.seekg( 0 );
      if( !input )
        DUNE_THROW(DGFException, "Error resetting input stream." );
      generate( input );
    }
 
    explicit GmshGridFactory ( const std::string &filename )
    {
      factory_ = new GridFactory();
      std::ifstream input( filename.c_str() );
      if( !input )
        DUNE_THROW( DGFException, "Macrofile " << filename << " not found." );
      if( !generate( input ) )
        DUNE_THROW( DGFException, "Could not generate MMesh from macrofile " << filename );
      input.close();
    }
 
    explicit GmshGridFactory ( Dune::GridFactory< Grid >& gridFactory, const std::string &filename )
    {
      factory_ = &gridFactory;
      std::ifstream input( filename.c_str() );
      if( !input )
        DUNE_THROW( DGFException, "Macrofile " << filename << " not found." );
      if( !generate( input ) )
        DUNE_THROW( DGFException, "Could not generate MMesh from macrofile " << filename );
      input.close();
    }
 
    std::shared_ptr<Grid> grid() const
    {
      if ( !grid_ )
        grid_ = factory_->createGrid();
 
      return grid_;
    }
 
    template< class Intersection >
    bool wasInserted ( const Intersection &intersection ) const
    {
      return factory_->wasInserted( intersection );
    }
 
  private:
    bool generate( std::istream &gridFile )
    {
      // read file until we get to the list of nodes
      std::string line;
      std::getline(gridFile, line);
      while (line.find("$Nodes") == std::string::npos)
      {
        if (!std::getline(gridFile, line))
          DUNE_THROW(IOError, "Gmsh file not valid!");
      }
 
      // read all vertices
      std::getline(gridFile, line);
      const auto numVertices = convertString<std::size_t>(line);
 
      std::getline(gridFile, line);
      std::size_t vertexCount = 0;
      while (line.find("$EndNodes")  == std::string::npos)
      {
          // drop first entry in line (vertex index) and read coordinates
          std::istringstream stream(line);
          std::string buf; stream >> buf;
          GlobalPosition v;
          for (auto& coord : v)
          {
              stream >> coord;
              if (stream.fail()) DUNE_THROW(Dune::IOError, "Could not read vertex coordinate");
          }
 
          // insert vertex and move to next line
          factory_->insertVertex( v );
          if (!std::getline(gridFile, line))
            DUNE_THROW(IOError, "Gmsh file not valid!");
          vertexCount++;
      }
 
      // we should always find as many vertices as the mesh file states
      if (vertexCount != numVertices)
          DUNE_THROW(Dune::InvalidStateException, "Couldn't find as many vertices as stated in the .msh file");
 
      // read file until we get to the list of elements
      while(line.find("$Elements") == std::string::npos)
        if (!std::getline(gridFile, line))
          DUNE_THROW(IOError, "Gmsh file not valid!");
 
      // read elements
      std::getline(gridFile, line);
      const auto numElements = convertString<std::size_t>(line);
 
      std::size_t elemCount = 0;
      std::getline(gridFile, line);
      while (line.find("$EndElements") == std::string::npos)
      {
          // pass all indices into vector
          std::istringstream stream(line);
          std::string buf;
          std::vector<std::size_t> lineData;
          while (stream >> buf) lineData.push_back(convertString<std::size_t>(buf));
          assert(lineData.size() >= 4 && "Grid format erroneous or unsupported");
 
          // obtain geometry type
          const auto gt = obtainGeometryType( lineData[1] );
 
          std::vector< unsigned int > cornersIndices;
          auto it = lineData.begin()+2+lineData[2]+1;
          for (; it != lineData.end(); ++it)
              cornersIndices.push_back(*it-1); // gmsh indices start from 1
 
          // insert element
          if ( gt.dim() == dimension )
              factory_->insertElement( gt, cornersIndices, lineData[2+lineData[2]] );
 
          // insert interface/boundary segments
          if ( gt.dim() == dimension-1 )
          {
              factory_->insertInterface( cornersIndices, lineData[2+lineData[2]] );
              factory_->insertBoundarySegment( cornersIndices );
          }
 
          // insert interface's boundary segments
          if ( gt.dim() == dimension-2 )
              factory_->insertInterfaceBoundarySegment( cornersIndices );
 
          // get next line
          if (!std::getline(gridFile, line))
            DUNE_THROW(IOError, "Gmsh file not valid!");
          elemCount++;
      }
 
      // make sure we read all elements
      if (elemCount != numElements)
          DUNE_THROW(Dune::InvalidStateException, "Didn't read as many elements as stated in the .msh file");
 
      return true;
    }
 
  private:
    template<class T>
    T convertString(const std::string& string) const
    {
        T value;
        std::istringstream stream(string);
        stream >> value;
        if (stream.fail())
            DUNE_THROW(Dune::InvalidStateException, "Couldn't convert string: " << string << "to type: " << typeid(T).name());
        return value;
    }
 
    Dune::GeometryType obtainGeometryType(std::size_t gmshElemType) const
    {
        switch (gmshElemType)
        {
            case 15: return Dune::GeometryTypes::vertex;        // points
            case 1:  return Dune::GeometryTypes::line;          // lines
            case 2:  return Dune::GeometryTypes::triangle;      // triangle
            case 3:  return Dune::GeometryTypes::quadrilateral; // quadrilateral
            case 4:  return Dune::GeometryTypes::tetrahedron;   // tetrahedron
            case 5:  return Dune::GeometryTypes::hexahedron;    // hexahedron
            default:
                DUNE_THROW(Dune::NotImplemented, "FacetCoupling gmsh reader for gmsh element type " << gmshElemType);
        }
    }
 
    mutable std::shared_ptr<Grid> grid_;
    GridFactory* factory_;
  };
 
} // end namespace Dune
 
#endif