gridfactory.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=2 sw=2 sts=2:
#ifndef DUNE_POLYHEDRALGRID_GRIDFACTORY_HH
#define DUNE_POLYHEDRALGRID_GRIDFACTORY_HH
 
#include <algorithm>
#include <numeric>
 
#include <dune/common/typetraits.hh>
#include <dune/common/version.hh>
 
#include <dune/grid/common/gridfactory.hh>
#include <opm/grid/polyhedralgrid/grid.hh>
 
namespace Dune
{
 
 
  // GridFactory for PolyhedralGrid
  // ---------------------------------
 
  template< int dim, int dimworld, class coord_t >
  class GridFactory< PolyhedralGrid< dim, dimworld, coord_t > >
    : public GridFactoryInterface< PolyhedralGrid< dim, dimworld, coord_t > >
  {
  public:
    typedef PolyhedralGrid< dim, dimworld, coord_t > Grid;
 
    const static int dimension      = Grid::dimension;
    const static int dimensionworld = Grid::dimensionworld;
    typedef typename Grid::ctype ctype;
 
    typedef MPIHelper::MPICommunicator MPICommunicatorType;
    typedef typename Grid::template Codim<0>::Entity Element;
    typedef typename Grid::template Codim<dimension>::Entity Vertex;
 
    typedef Dune::FieldVector<ctype,dimensionworld> CoordinateType;
    typedef CoordinateType  Coordinate;
 
    using UniquePtrType = std::unique_ptr<Grid>;
 
    explicit GridFactory ( const MPICommunicatorType& = MPIHelper::getCommunicator() )
      : nodes_(),
        faces_(),
        cells_()
    {}
 
    virtual void insertVertex(const CoordinateType& pos)
    {
      nodes_.push_back( pos );
    }
 
    virtual void insertElement(const GeometryType& type,
                               const std::vector<unsigned int>& items)
    {
      if( type.isNone() )
      {
        // copy into vector of integers
        std::vector< int > numbers( items.size() );
        std::copy( items.begin(), items.end(), numbers.begin() );
 
        if( type.dim() == dimension-1 )
        {
          faces_.push_back( numbers );
        }
        else if( type.dim() == dimension )
        {
          // note vertices holds the face
          // numbers in this case
          cells_.push_back( numbers );
        }
        else
        {
          DUNE_THROW(Dune::NotImplemented,"insertElement not implemented for type " << type );
        }
      }
      else // use ReferenceElement to insert faces
      {
 
 
      }
    }
 
    void insertBoundarySegment(const std::vector<unsigned int>&)
    {
      DUNE_THROW(NotImplemented,"yet");
    }
 
    UniquePtrType createGrid()
    {
      std::vector< CoordinateType >& nodes = nodes_;
      std::vector< std::vector< int > >& faces = faces_;
      std::vector< std::vector< int > >& cells = cells_;
 
      if( cells.empty() )
      {
        DUNE_THROW( GridError, "No cells found for PolyhedralGrid" );
      }
 
      const auto sumSize = [] ( std::size_t s, const std::vector< int > &v ) { return s + v.size(); };
      const std::size_t numFaceNodes = std::accumulate( faces.begin(), faces.end(), std::size_t( 0 ), sumSize );
      const std::size_t numCellFaces = std::accumulate( cells.begin(), cells.end(), std::size_t( 0 ), sumSize );
 
      typename Grid::UnstructuredGridPtr ug =
        Grid::allocateGrid( cells.size(), faces.size(), numFaceNodes, numCellFaces, nodes.size() );
 
      // copy faces
      {
#ifndef NDEBUG
        std::map< std::vector< int >, std::vector< int > > faceMap;
#endif
 
        const int nFaces = faces.size();
        // set all face_cells values to -2 as default
        std::fill( ug->face_cells, ug->face_cells + 2*nFaces, -1 );
 
        int facepos = 0;
        std::vector< int > faceVertices;
        faceVertices.reserve( 30 );
        for( int face = 0; face < nFaces; ++face )
        {
          //std::cout << "face " << face << ": ";
          faceVertices.clear();
          ug->face_nodepos[ face ] = facepos;
          const int nVertices = faces[ face ].size();
          for( int vx = 0; vx < nVertices; ++vx, ++facepos )
          {
            //std::cout << " " << faces[ face ][ vx ];
            ug->face_nodes[ facepos ] = faces[ face ][ vx ];
            faceVertices.push_back( faces[ face ][ vx ] );
          }
          //std::cout << std::endl;
 
#ifndef NDEBUG
          // sort vertices
          std::sort( faceVertices.begin(), faceVertices.end() );
          // make sure each face only exists once
          faceMap[ faceVertices ].push_back( face );
          assert( faceMap[ faceVertices ].size() == 1 );
#endif
        }
        ug->face_nodepos[ nFaces ] = facepos ;
      }
 
      // copy cells
      {
        const int nCells = cells.size();
        int cellpos = 0;
        for( int cell = 0; cell < nCells; ++cell )
        {
          //std::cout << "Cell " << cell << ": ";
          ug->cell_facepos[ cell ] = cellpos;
          const int nFaces = cells[ cell ].size();
          for( int f = 0; f < nFaces; ++f, ++cellpos )
          {
            const int face = cells[ cell ][ f ];
            // std::cout << " " << face ;
            ug->cell_faces[ cellpos ] = face;
 
            // TODO find cells for each face
            if( ug->face_cells[ 2*face ] == -1 )
            {
              ug->face_cells[ 2*face ] = cell;
            }
            else // if ( ug->face_cells[ 2*face+1 ] == -1 )
            {
              //assert( ug->face_cells[ 2*face+1 ] == -1 );
              ug->face_cells[ 2*face+1 ] = cell;
            }
          }
          //std::cout << std::endl;
        }
        ug->cell_facepos[ nCells ] = cellpos ;
      }
 
      // copy node coordinates
      {
        const int nNodes = nodes.size();
        int nodepos = 0;
        for( int vx = 0 ; vx < nNodes; ++vx )
        {
          for( int d=0; d<dim; ++d, ++nodepos )
            ug->node_coordinates[ nodepos ] = nodes[ vx ][ d ];
        }
      }
 
      /*
      for( int i=0; i<int(faces.size() ); ++i)
      {
        std::cout << "face "<< i<< " connects to " << ug->face_cells[ 2*i ] << " " <<
          ug->face_cells[ 2*i+1] << std::endl;
      }
      */
 
      // free cell face tag since it's not a cartesian grid
      if( ug->cell_facetag )
      {
        std::free( ug->cell_facetag );
        ug->cell_facetag = nullptr ;
        for( int i=0; i<3; ++i ) ug->cartdims[ i ] = 0;
      }
 
      // compute geometric quantities like cell volume and face normals
      Grid::computeGeometry( ug );
 
      // check normal direction
      {
        for( int face = 0 ; face < ug->number_of_faces; ++face )
        {
          const int a = ug->face_cells[ 2*face     ];
          const int b = ug->face_cells[ 2*face + 1 ];
          if( a < 0 || b < 0 )
            continue ;
 
          Coordinate centerDiff( 0 );
          Coordinate normal( 0 );
          //std::cout << "Cell center " << a << " " << b << std::endl;
          for( int d=0; d<dim; ++d )
          {
            //std::cout << ug->cell_centroids[ a*dim + d ] << " " << ug->cell_centroids[ b*dim + d ] << std::endl;
            centerDiff[ d ] = ug->cell_centroids[ b*dim + d ] - ug->cell_centroids[ a*dim + d ];
            normal[ d ] = ug->face_normals[ face*dim + d ];
          }
 
          // if diff and normal point in different direction, flip faces
          if( centerDiff * normal > 0 )
          {
            ug->face_cells[ 2*face     ] = b;
            ug->face_cells[ 2*face + 1 ] = a;
          }
        }
      }
 
      return UniquePtrType( new Grid( std::move( ug ) ));
    }
 
  protected:
    std::vector< CoordinateType > nodes_;
    std::vector< std::vector< int > > faces_;
    std::vector< std::vector< int > > cells_;
  };
 
} // namespace Dune
 
#endif // #ifndef DUNE_POLYHEDRALGRID_DGFPARSER_HH