Adds Verteq (Vertical Equilibrium) Grid

CMakeLists_files.cmake

@@ -58,6 +58,10 @@ list (APPEND MAIN_SOURCE_FILES
   opm/grid/utility/cartesianToCompressed.cpp
   opm/grid/utility/StopWatch.cpp
   opm/grid/utility/WachspressCoord.cpp
+  opm/grid/verteq/topsurf.cpp
+  opm/grid/verteq/nav.cpp
+  opm/grid/verteq/utility/exc.cpp
+  opm/grid/verteq/utility/runlen.cpp
   )
 
 if (opm-common_FOUND)
@@ -173,6 +177,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/grid/common/WellConnections.hpp
   opm/grid/common/ZoltanGraphFunctions.hpp
   opm/grid/common/ZoltanPartition.hpp
+  opm/grid/polyhedralgrid.hh
   opm/grid/polyhedralgrid/capabilities.hh
   opm/grid/polyhedralgrid/cartesianindexmapper.hh
   opm/grid/polyhedralgrid/declaration.hh
@@ -184,7 +189,6 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/grid/polyhedralgrid/gridhelpers.hh
   opm/grid/polyhedralgrid/grid.hh
   opm/grid/polyhedralgrid/gridview.hh
-  opm/grid/polyhedralgrid.hh
   opm/grid/polyhedralgrid/idset.hh
   opm/grid/polyhedralgrid/indexset.hh
   opm/grid/polyhedralgrid/intersection.hh
@@ -221,4 +225,9 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/grid/utility/OpmParserIncludes.hpp
   opm/grid/utility/platform_dependent/disable_warnings.h
   opm/grid/utility/platform_dependent/reenable_warnings.h
+  opm/grid/verteq/topsurf.hpp
+  opm/grid/verteq/nav.hpp
+  opm/grid/verteq/utility/visibility.h
+  opm/grid/verteq/utility/exc.hpp
+  opm/grid/verteq/utility/runlen.hpp
   )

opm/grid/UnstructuredGrid.c

@@ -166,6 +166,67 @@ allocate_grid(size_t ndims     ,
     return G;
 }
 
+#define SWAP(a,b) { int tmp = a; a = b; b =tmp; }
+
+void print_grid(const struct UnstructuredGrid *grd )
+{
+  int dim = grd->dimensions;
+  for(int i=0; i<grd->number_of_nodes; ++i )
+  {
+    printf("Vx (%d) = ( ",i);
+    for( int d=0; d<dim; ++d )
+    {
+      printf(" %2.2f ", grd->node_coordinates[i*dim + d ]);
+    }
+    printf(")\n");
+  }
+
+  for( int f=0; f<grd->number_of_faces; ++f )
+  {
+    //for( int f=grd->cell_facepos[ c ]; f<grd->cell_facepos[ c+1 ]; ++f )
+    printf(" Face (%d) (",f);
+    for( int vx=grd->face_nodepos[ f ]; vx<grd->face_nodepos[ f+1 ]; ++vx )
+    {
+      printf(" %d ",grd->face_nodes[ vx ]);
+    }
+    printf(")");
+    int bnd = grd->face_cells[ 2*f ] < 0 || grd->face_cells[ 2*f+1 ] < 0;
+    if( bnd )
+      printf(" boundary");
+    printf("\n");
+  }
+
+  /*
+  if( dim == 2 )
+  {
+    for( int c=0; c<grd->number_of_cells; ++c )
+    {
+      int f=grd->cell_facepos[ c ];
+      SWAP( grd->cell_faces[ f+1 ], grd->cell_faces[ f+2 ] );
+      SWAP( grd->cell_facetag[ f+1 ], grd->cell_facetag[ f+2 ] );
+    }
+  }
+  */
+
+  for( int c=0; c<grd->number_of_cells; ++c )
+  {
+    printf("Cell %d = (", c);
+    //for( int f=grd->cell_facepos[ c ]; f<grd->cell_facepos[ c+1 ]; ++f )
+    for( int f=grd->cell_facepos[ c ]; f<grd->cell_facepos[ c+1 ]; ++f )
+    {
+      printf(" %d ",grd->cell_faces[ f ]);
+
+      /*
+      for( int vx=grd->face_nodepos[ f ]; vx<grd->face_nodepos[ f+1 ]; ++vx )
+      {
+        printf(" %d ",grd->face_nodes[ vx ]);
+      }
+      printf(")\n");
+      */
+    }
+    printf(")\n");
+  }
+}
 
 #define GRID_NMETA      6
 #define GRID_NDIMS      0

opm/grid/UnstructuredGrid.h

@@ -299,6 +299,8 @@ allocate_grid(size_t ndims     ,
               size_t ncellfaces,
               size_t nnodes    );
 
+void print_grid(const struct UnstructuredGrid *);
+
 
 /**
    Will allocate storage internally in the grid object to hold a copy

opm/grid/common/VerteqColumnUtility.hpp

@@ -0,0 +1,56 @@
+#ifndef OPM_VERTEQCOLUMNUTILITY_HEADER
+#define OPM_VERTEQCOLUMNUTILITY_HEADER
+
+#include <array>
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+#include <opm/grid/verteq/topsurf.hpp>
+
+
+namespace Dune
+{
+
+    class ColumnCell
+    {
+        typedef Opm::TopSurf  TopSurfaceGridType;
+        const TopSurfaceGridType& topSurf_;
+        const int index_;
+
+    public:
+        ColumnCell( const TopSurfaceGridType& topSurf, const int index )
+          : topSurf_( topSurf ),
+            index_( index )
+        {}
+
+        int index () const { return index_; }
+
+        double dz () const { return topSurf_.dz[ index_ ]; }
+        double h () const  { return topSurf_.h [ index_ ]; }
+
+        int maxVertRes() const { return topSurf_.max_vert_res; }
+
+        // this can be used to access an entity in the original 3D grid
+        int fineCellIndex() const { return topSurf_.fine_col[ index_ ]; }
+    };
+
+
+    /** \brief Interface class to access the logical Cartesian grid as used in industry
+               standard simulator decks.
+               */
+    template< class Grid >
+    class VerteqColumnUtility
+    {
+    public:
+        /** \brief dimension of the grid */
+        static const int dimension = Grid :: dimension ;
+
+        /** \brief constructor taking grid */
+        explicit VerteqColumnUtility( const Grid& )
+        {
+            DUNE_THROW(InvalidStateException,"CartesianIndexMapper not specialized for given grid");
+        }
+    };
+
+} // end namespace Opm
+#endif

opm/grid/polyhedralgrid/geometry.hh

@@ -191,6 +191,7 @@ namespace Dune
       return storage_.volume();
     }
 
+
     JacobianTransposed jacobianTransposed ( const LocalCoordinate & local ) const
     {
       if( geometryImpl_ )

opm/grid/polyhedralgrid/grid.hh

@@ -30,6 +30,7 @@
 #include <opm/grid/polyhedralgrid/geometry.hh>
 #include <opm/grid/polyhedralgrid/gridview.hh>
 #include <opm/grid/polyhedralgrid/idset.hh>
+#include <opm/grid/polyhedralgrid/polyhedralmesh.hh>
 
 // Re-enable warnings.
 #include <opm/grid/utility/platform_dependent/reenable_warnings.h>
@@ -43,6 +44,10 @@
 #include <opm/grid/cornerpoint_grid.h>
 #include <opm/grid/MinpvProcessor.hpp>
 
+#include <opm/grid/verteq/topsurf.hpp>
+
+#include <opm/grid/verteq/topsurf.hpp>
+
 namespace Dune
 {
 
@@ -163,6 +168,7 @@ namespace Dune
 
   public:
     typedef UnstructuredGrid  UnstructuredGridType;
+    typedef Opm::TopSurf      TopSurfaceGridType;
 
   protected:
     struct UnstructuredGridDeleter
@@ -174,6 +180,7 @@ namespace Dune
     };
 
   public:
+    typedef PolyhedralMesh< dim, dimworld, coord_t, int > PolyhedralMeshType;
     typedef std::unique_ptr< UnstructuredGridType, UnstructuredGridDeleter > UnstructuredGridPtr;
 
     static UnstructuredGridPtr
@@ -344,6 +351,7 @@ namespace Dune
                               const std::vector< double >& dx )
     : gridPtr_( createGrid( n, dx ) ),
       grid_( *gridPtr_ ),
+      topSurfaceGrid_( nullptr ),
       comm_( MPIHelper::getCommunicator()),
       leafIndexSet_( *this ),
       globalIdSet_( *this ),
@@ -362,6 +370,8 @@ namespace Dune
     explicit PolyhedralGrid ( UnstructuredGridPtr &&gridPtr )
     : gridPtr_( std::move( gridPtr ) ),
       grid_( *gridPtr_ ),
+      topSurfaceGrid_( nullptr ),
+      //polyhedralMesh_( grid_ ),
       comm_( MPIHelper::getCommunicator() ),
       leafIndexSet_( *this ),
       globalIdSet_( *this ),
@@ -381,6 +391,7 @@ namespace Dune
     explicit PolyhedralGrid ( const UnstructuredGridType& grid )
     : gridPtr_(),
       grid_( grid ),
+      topSurfaceGrid_( nullptr ),
       comm_( MPIHelper::getCommunicator() ),
       leafIndexSet_( *this ),
       globalIdSet_( *this ),
@@ -390,6 +401,27 @@ namespace Dune
       init();
     }
 
+    /** \brief constructor
+     *
+     *  The references to ug are stored in the grid.
+     *  Therefore, they must remain valid until the grid is destroyed.
+     *
+     *  \param[in]  ug    UnstructuredGrid reference
+     */
+    explicit PolyhedralGrid ( const TopSurfaceGridType& topSurf )
+    : gridPtr_(),
+      grid_( static_cast< UnstructuredGridType > ( topSurf ) ),
+      topSurfaceGrid_( &topSurf ),
+      comm_( *this ),
+      leafIndexSet_( *this ),
+      globalIdSet_( *this ),
+      localIdSet_( *this ),
+      nBndSegments_( 0 )
+    {
+      // std::cout << "Creating TopSurfaceGrid" <<  topSurfaceGrid_ << std::endl;
+      init();
+    }
+
     /** \} */
 
     /** \name Casting operators
@@ -398,6 +430,8 @@ namespace Dune
 
     /** \} */
 
+    const TopSurfaceGridType* topSurfaceGrid() const { return topSurfaceGrid_; }
+
     /** \name Size Methods
      *  \{ */
 
@@ -434,6 +468,7 @@ namespace Dune
      */
     int size ( int codim ) const
     {
+      //return polyhedralMesh_.size( codim );
       if( codim == 0 )
       {
         return grid_.number_of_cells;
@@ -1353,6 +1388,8 @@ namespace Dune
   protected:
     void init ()
     {
+      //std::cout << "PolyhedralGrid init" << std::endl;
+
       // copy Cartesian dimensions
       for( int i=0; i<3; ++i )
       {
@@ -1465,6 +1502,11 @@ namespace Dune
               cellVertices_[ c ][ (*vx).second ] = (*it).first ;
             }
           }
+
+          for( int c=0; c<numCells; ++c )
+          {
+            // sort face_nodes according to reference element
+          }
         }
 
         // if face_tag is available we assume that the elements follow a cube-like structure
@@ -1617,6 +1659,15 @@ namespace Dune
 
             hasCube = true;
           }
+          else if ( codim == 0 )
+          {
+            //if( minVx == maxVx && maxVx == 8 )
+            // tmp.makeCube(dim);
+            //if( minVx == maxVx && maxVx == 4 )
+            //  tmp.makeSimplex(dim);
+            //else
+            tmp.makeNone( dim );
+          }
           else
           {
             hasPolyhedron = true;
@@ -1691,17 +1742,6 @@ namespace Dune
           }
         }
       }
-
-      /*
-      for( int i=0; i<= dim ; ++ i )
-      {
-        for( const auto& geomType : geomTypes_[ i ] )
-        {
-          std::cout << "Codim " << i << "  type = " << geomType << std::endl;
-        }
-      }
-      */
-      //print( std::cout, grid_ );
     }
 
     void print( std::ostream& out, const UnstructuredGridType& grid ) const
@@ -1731,13 +1771,16 @@ namespace Dune
           out << vx[ i ] << " ";
         out << std::endl;
       }
-
     }
 
   protected:
     UnstructuredGridPtr gridPtr_;
     const UnstructuredGridType& grid_;
 
+    const TopSurfaceGridType* topSurfaceGrid_;
+
+    PolyhedralMeshType polyhedralMesh_;
+
     CollectiveCommunication comm_;
     std::array< int, 3 > cartDims_;
     std::vector< std::vector< GeometryType > > geomTypes_;
@@ -1852,5 +1895,6 @@ namespace Dune
 #include <opm/grid/polyhedralgrid/persistentcontainer.hh>
 #include <opm/grid/polyhedralgrid/cartesianindexmapper.hh>
 #include <opm/grid/polyhedralgrid/gridhelpers.hh>
+#include <opm/grid/polyhedralgrid/verteqcolumnutility.hh>
 
 #endif // #ifndef DUNE_POLYHEDRALGRID_GRID_HH

opm/grid/polyhedralgrid/gridfactory.hh

@@ -18,7 +18,6 @@ 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 > >

opm/grid/polyhedralgrid/idset.hh

@@ -48,6 +48,15 @@ namespace Dune
       }
     }
 
+#if ! DUNE_VERSION_NEWER(DUNE_GRID,2,4)
+    //! id meethod for entity and specific codim
+    template< int codim >
+    IdType id ( const typename Traits::template Codim< codim >::EntityPointer &entityPointer ) const
+    {
+      return id( *entityPointer );
+    }
+#endif
+
     //! id method of all entities
     template< class Entity >
     IdType id ( const Entity &entity ) const