My Project
grid.hh
1 // -*- mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=2 sw=2 sts=2:
3 #ifndef DUNE_POLYHEDRALGRID_GRID_HH
4 #define DUNE_POLYHEDRALGRID_GRID_HH
5 
6 #include <set>
7 #include <vector>
8 
9 // Warning suppression for Dune includes.
10 #include <opm/grid/utility/platform_dependent/disable_warnings.h>
11 
12 //- dune-common includes
13 #include <dune/common/version.hh>
14 #include <dune/common/parallel/mpihelper.hh>
15 
16 //- dune-grid includes
17 #include <dune/grid/common/grid.hh>
18 
19 #if DUNE_VERSION_GTE(DUNE_COMMON, 2, 7)
20 #include <dune/common/parallel/communication.hh>
21 #else
22 #include <dune/common/parallel/collectivecommunication.hh>
23 #endif
24 
25 //- polyhedralgrid includes
26 #include <opm/grid/polyhedralgrid/capabilities.hh>
27 #include <opm/grid/polyhedralgrid/declaration.hh>
28 #include <opm/grid/polyhedralgrid/entity.hh>
29 #include <opm/grid/polyhedralgrid/entityseed.hh>
30 #include <opm/grid/polyhedralgrid/geometry.hh>
31 #include <opm/grid/polyhedralgrid/gridview.hh>
32 #include <opm/grid/polyhedralgrid/idset.hh>
33 
34 // Re-enable warnings.
35 #include <opm/grid/utility/platform_dependent/reenable_warnings.h>
36 
37 #include <opm/grid/utility/ErrorMacros.hpp>
38 
40 #include <opm/grid/cart_grid.h>
42 #include <opm/grid/GridManager.hpp>
44 #include <opm/grid/MinpvProcessor.hpp>
45 
46 namespace Dune
47 {
48 
49 
50  // PolyhedralGridFamily
51  // ------------
52 
53  template< int dim, int dimworld, typename coord_t >
55  {
56  struct Traits
57  {
59 
60  typedef coord_t ctype;
61 
62  // type of data passed to entities, intersections, and iterators
63  // for PolyhedralGrid this is just an empty place holder
64  typedef const Grid* ExtraData;
65 
66  typedef int Index ;
67 
68  static const int dimension = dim;
69  static const int dimensionworld = dimworld;
70 
71  typedef Dune::FieldVector< ctype, dimensionworld > GlobalCoordinate ;
72 
77 
78  typedef Dune::Intersection< const Grid, LeafIntersectionImpl > LeafIntersection;
79  typedef Dune::Intersection< const Grid, LevelIntersectionImpl > LevelIntersection;
80 
81  typedef Dune::IntersectionIterator< const Grid, LeafIntersectionIteratorImpl, LeafIntersectionImpl > LeafIntersectionIterator;
82  typedef Dune::IntersectionIterator< const Grid, LevelIntersectionIteratorImpl, LevelIntersectionImpl > LevelIntersectionIterator;
83 
85  typedef Dune::EntityIterator< 0, const Grid, HierarchicIteratorImpl > HierarchicIterator;
86 
87  template< int codim >
88  struct Codim
89  {
90  typedef PolyhedralGridGeometry<dimension-codim, dimensionworld, const Grid> GeometryImpl;
91  typedef Dune::Geometry< dimension-codim, dimensionworld, const Grid, PolyhedralGridGeometry > Geometry;
92 
93  typedef PolyhedralGridLocalGeometry< dimension-codim, dimensionworld, const Grid> LocalGeometryImpl;
94  typedef Dune::Geometry< dimension-codim, dimensionworld, const Grid, PolyhedralGridLocalGeometry > LocalGeometry;
95 
97  typedef Dune::Entity< codim, dimension, const Grid, PolyhedralGridEntity > Entity;
98 
100  typedef Entity EntityPointer;
101 
102  //typedef Dune::EntitySeed< const Grid, PolyhedralGridEntitySeed< codim, const Grid > > EntitySeed;
104 
105  template< PartitionIteratorType pitype >
106  struct Partition
107  {
109  typedef Dune::EntityIterator< codim, const Grid, LeafIteratorImpl > LeafIterator;
110 
111  typedef LeafIterator LevelIterator;
112  };
113 
114  typedef typename Partition< All_Partition >::LeafIterator LeafIterator;
115  typedef typename Partition< All_Partition >::LevelIterator LevelIterator;
116  };
117 
120 
122  typedef GlobalIdSet LocalIdSet;
123 
124  typedef Dune::MPIHelper::MPICommunicator MPICommunicator;
125  typedef Dune::CollectiveCommunication<MPICommunicator> CollectiveCommunication;
126 
127  template< PartitionIteratorType pitype >
128  struct Partition
129  {
130  typedef Dune::GridView< PolyhedralGridViewTraits< dim, dimworld, ctype, pitype > > LeafGridView;
131  typedef Dune::GridView< PolyhedralGridViewTraits< dim, dimworld, ctype, pitype > > LevelGridView;
132  };
133 
134  typedef typename Partition< All_Partition >::LeafGridView LeafGridView;
135  typedef typename Partition< All_Partition >::LevelGridView LevelGridView;
136  };
137  };
138 
139 
140 
141  // PolyhedralGrid
142  // --------------
143 
152  template < int dim, int dimworld, typename coord_t >
155  : public GridDefaultImplementation
156  < dim, dimworld, coord_t, PolyhedralGridFamily< dim, dimworld, coord_t > >
158  {
160 
161  typedef GridDefaultImplementation
162  < dim, dimworld, coord_t, PolyhedralGridFamily< dim, dimworld, coord_t > > Base;
163 
164  public:
166 
167  protected:
169  {
170  inline void operator () ( UnstructuredGridType* grdPtr )
171  {
172  destroy_grid( grdPtr );
173  }
174  };
175 
176  public:
177  typedef std::unique_ptr< UnstructuredGridType, UnstructuredGridDeleter > UnstructuredGridPtr;
178 
179  static UnstructuredGridPtr
180  allocateGrid ( std::size_t nCells, std::size_t nFaces, std::size_t nFaceNodes, std::size_t nCellFaces, std::size_t nNodes )
181  {
182  // Note that we here assign a grid of dimension dimworld in order to obtain global coordinates in the correct dimension
183  UnstructuredGridType *grid = allocate_grid( dimworld, nCells, nFaces, nFaceNodes, nCellFaces, nNodes );
184  if( !grid )
185  DUNE_THROW( GridError, "Unable to allocate grid" );
186  return UnstructuredGridPtr( grid );
187  }
188 
189  static void
190  computeGeometry ( UnstructuredGridPtr& ug )
191  {
192  // get C pointer to UnstructuredGrid
193  UnstructuredGrid* ugPtr = ug.operator ->();
194 
195  // compute geometric quantities like cell volume and face normals
196  compute_geometry( ugPtr );
197  }
198 
200  typedef PolyhedralGridFamily< dim, dimworld, coord_t > GridFamily;
207  typedef typename GridFamily::Traits Traits;
208 
215  template< int codim >
216  struct Codim;
217 
224  typedef typename Traits::HierarchicIterator HierarchicIterator;
226  typedef typename Traits::LeafIntersectionIterator LeafIntersectionIterator;
228  typedef typename Traits::LevelIntersectionIterator LevelIntersectionIterator;
229 
236  template< PartitionIteratorType pitype >
237  struct Partition
238  {
239  typedef typename GridFamily::Traits::template Partition< pitype >::LevelGridView LevelGridView;
240  typedef typename GridFamily::Traits::template Partition< pitype >::LeafGridView LeafGridView;
241  };
242 
244  typedef typename Partition< All_Partition >::LevelGridView LevelGridView;
245  typedef typename Partition< All_Partition >::LeafGridView LeafGridView;
246 
260  typedef typename Traits::LeafIndexSet LeafIndexSet;
261 
270  typedef typename Traits::LevelIndexSet LevelIndexSet;
271 
282  typedef typename Traits::GlobalIdSet GlobalIdSet;
283 
299  typedef typename Traits::LocalIdSet LocalIdSet;
300 
307  typedef typename Traits::ctype ctype;
308 
310  typedef typename Traits::CollectiveCommunication CollectiveCommunication;
311 
312  typedef typename Traits :: GlobalCoordinate GlobalCoordinate;
313 
319 #if HAVE_ECL_INPUT
325  explicit PolyhedralGrid ( const Opm::EclipseGrid& inputGrid,
326  const std::vector<double>& poreVolumes = std::vector<double> ())
327  : gridPtr_( createGrid( inputGrid, poreVolumes ) ),
328  grid_( *gridPtr_ ),
329  comm_( MPIHelper::getCommunicator() ),
330  leafIndexSet_( *this ),
331  globalIdSet_( *this ),
332  localIdSet_( *this ),
333  nBndSegments_( 0 )
334  {
335  init();
336  }
337 #endif
338 
344  explicit PolyhedralGrid ( const std::vector< int >& n,
345  const std::vector< double >& dx )
346  : gridPtr_( createGrid( n, dx ) ),
347  grid_( *gridPtr_ ),
348  comm_( MPIHelper::getCommunicator()),
349  leafIndexSet_( *this ),
350  globalIdSet_( *this ),
351  localIdSet_( *this ),
352  nBndSegments_( 0 )
353  {
354  init();
355  }
356 
363  explicit PolyhedralGrid ( UnstructuredGridPtr &&gridPtr )
364  : gridPtr_( std::move( gridPtr ) ),
365  grid_( *gridPtr_ ),
366  comm_( MPIHelper::getCommunicator() ),
367  leafIndexSet_( *this ),
368  globalIdSet_( *this ),
369  localIdSet_( *this ),
370  nBndSegments_( 0 )
371  {
372  init();
373  }
374 
382  explicit PolyhedralGrid ( const UnstructuredGridType& grid )
383  : gridPtr_(),
384  grid_( grid ),
385  comm_( MPIHelper::getCommunicator() ),
386  leafIndexSet_( *this ),
387  globalIdSet_( *this ),
388  localIdSet_( *this ),
389  nBndSegments_( 0 )
390  {
391  init();
392  }
393 
398  operator const UnstructuredGridType& () const { return grid_; }
399 
412  int maxLevel () const
413  {
414  return 1;
415  }
416 
425  int size ( int /* level */, int codim ) const
426  {
427  return size( codim );
428  }
429 
436  int size ( int codim ) const
437  {
438  if( codim == 0 )
439  {
440  return grid_.number_of_cells;
441  }
442  else if ( codim == 1 )
443  {
444  return grid_.number_of_faces;
445  }
446  else if ( codim == dim )
447  {
448  return grid_.number_of_nodes;
449  }
450  else
451  {
452  std::cerr << "Warning: codimension " << codim << " not available in PolyhedralGrid" << std::endl;
453  return 0;
454  }
455  }
456 
465  int size ( int /* level */, GeometryType type ) const
466  {
467  return size( dim - type.dim() );
468  }
469 
474  int size ( GeometryType type ) const
475  {
476  return size( dim - type.dim() );
477  }
478 
485  size_t numBoundarySegments () const
486  {
487  return nBndSegments_;
488  }
491  template< int codim >
492  typename Codim< codim >::LeafIterator leafbegin () const
493  {
494  return leafbegin< codim, All_Partition >();
495  }
496 
497  template< int codim >
498  typename Codim< codim >::LeafIterator leafend () const
499  {
500  return leafend< codim, All_Partition >();
501  }
502 
503  template< int codim, PartitionIteratorType pitype >
504  typename Codim< codim >::template Partition< pitype >::LeafIterator
505  leafbegin () const
506  {
507  typedef typename Traits::template Codim< codim >::template Partition< pitype >::LeafIteratorImpl Impl;
508  return Impl( extraData(), true );
509  }
510 
511  template< int codim, PartitionIteratorType pitype >
512  typename Codim< codim >::template Partition< pitype >::LeafIterator
513  leafend () const
514  {
515  typedef typename Traits::template Codim< codim >::template Partition< pitype >::LeafIteratorImpl Impl;
516  return Impl( extraData(), false );
517  }
518 
519  template< int codim >
520  typename Codim< codim >::LevelIterator lbegin ( const int /* level */ ) const
521  {
522  return leafbegin< codim, All_Partition >();
523  }
524 
525  template< int codim >
526  typename Codim< codim >::LevelIterator lend ( const int /* level */ ) const
527  {
528  return leafend< codim, All_Partition >();
529  }
530 
531  template< int codim, PartitionIteratorType pitype >
532  typename Codim< codim >::template Partition< pitype >::LevelIterator
533  lbegin ( const int /* level */ ) const
534  {
535  return leafbegin< codim, pitype > ();
536  }
537 
538  template< int codim, PartitionIteratorType pitype >
539  typename Codim< codim >::template Partition< pitype >::LevelIterator
540  lend ( const int /* level */ ) const
541  {
542  return leafend< codim, pitype > ();
543  }
544 
545  const GlobalIdSet &globalIdSet () const
546  {
547  return globalIdSet_;
548  }
549 
550  const LocalIdSet &localIdSet () const
551  {
552  return localIdSet_;
553  }
554 
555  const LevelIndexSet &levelIndexSet ( int /* level */ ) const
556  {
557  return leafIndexSet();
558  }
559 
560  const LeafIndexSet &leafIndexSet () const
561  {
562  return leafIndexSet_;
563  }
564 
565  void globalRefine ( int /* refCount */ )
566  {
567  }
568 
569  bool mark ( int /* refCount */, const typename Codim< 0 >::Entity& /* entity */ )
570  {
571  return false;
572  }
573 
574  int getMark ( const typename Codim< 0 >::Entity& /* entity */) const
575  {
576  return false;
577  }
578 
580  bool preAdapt ()
581  {
582  return false;
583  }
584 
586  bool adapt ()
587  {
588  return false ;
589  }
590 
595  template< class DataHandle >
596  bool adapt ( DataHandle & )
597  {
598  return false;
599  }
600 
602  void postAdapt ()
603  {
604  }
605 
613  int overlapSize ( int /* codim */) const
614  {
615  return 0;
616  }
617 
622  int ghostSize( int codim ) const
623  {
624  return (codim == 0 ) ? 1 : 0;
625  }
626 
632  int overlapSize ( int /* level */, int /* codim */ ) const
633  {
634  return 0;
635  }
636 
642  int ghostSize ( int /* level */, int codim ) const
643  {
644  return ghostSize( codim );
645  }
646 
660  template< class DataHandle>
661  void communicate ( DataHandle& /* dataHandle */,
662  InterfaceType /* interface */,
663  CommunicationDirection /* direction */,
664  int /* level */ ) const
665  {
666  OPM_THROW(std::runtime_error, "communicate not implemented for polyhedreal grid!");
667  }
668 
681  template< class DataHandle>
682  void communicate ( DataHandle& /* dataHandle */,
683  InterfaceType /* interface */,
684  CommunicationDirection /* direction */ ) const
685  {
686  OPM_THROW(std::runtime_error, "communicate not implemented for polyhedreal grid!");
687  }
688 
691  {
692  OPM_THROW(std::runtime_error, "switch to global view not implemented for polyhedreal grid!");
693  }
694 
697  {
698  OPM_THROW(std::runtime_error, "switch to distributed view not implemented for polyhedreal grid!");
699  }
700 
710  {
711  return comm_;
712  }
713 
714  // data handle interface different between geo and interface
715 
725  bool loadBalance ()
726  {
727  return false ;
728  }
729 
745  template< class DataHandle, class Data >
746  bool loadBalance ( CommDataHandleIF< DataHandle, Data >& /* datahandle */ )
747  {
748  return false;
749  }
750 
765  template< class DofManager >
766  bool loadBalance ( DofManager& /* dofManager */ )
767  {
768  return false;
769  }
770 
772  template< PartitionIteratorType pitype >
773  typename Partition< pitype >::LevelGridView levelGridView ( int /* level */ ) const
774  {
775  typedef typename Partition< pitype >::LevelGridView View;
776  typedef typename View::GridViewImp ViewImp;
777  return View( ViewImp( *this ) );
778  }
779 
781  template< PartitionIteratorType pitype >
782  typename Partition< pitype >::LeafGridView leafGridView () const
783  {
784  typedef typename Traits::template Partition< pitype >::LeafGridView View;
785  typedef typename View::GridViewImp ViewImp;
786  return View( ViewImp( *this ) );
787  }
788 
790  LevelGridView levelGridView ( int /* level */ ) const
791  {
792  typedef typename LevelGridView::GridViewImp ViewImp;
793  return LevelGridView( ViewImp( *this ) );
794  }
795 
797  LeafGridView leafGridView () const
798  {
799  typedef typename LeafGridView::GridViewImp ViewImp;
800  return LeafGridView( ViewImp( *this ) );
801  }
802 
804  template< class EntitySeed >
805  typename Traits::template Codim< EntitySeed::codimension >::EntityPointer
806  entityPointer ( const EntitySeed &seed ) const
807  {
808  typedef typename Traits::template Codim< EntitySeed::codimension >::EntityPointer EntityPointer;
809  typedef typename Traits::template Codim< EntitySeed::codimension >::EntityPointerImpl EntityPointerImpl;
810  typedef typename Traits::template Codim< EntitySeed::codimension >::EntityImpl EntityImpl;
811  return EntityPointer( EntityPointerImpl( EntityImpl( extraData(), seed ) ) );
812  }
813 
815  template< class EntitySeed >
816  typename Traits::template Codim< EntitySeed::codimension >::Entity
817  entity ( const EntitySeed &seed ) const
818  {
819  typedef typename Traits::template Codim< EntitySeed::codimension >::EntityImpl EntityImpl;
820  return EntityImpl( extraData(), seed );
821  }
822 
836  void update ()
837  {
838  }
839 
842  const std::array<int, 3>& logicalCartesianSize() const
843  {
844  return cartDims_;
845  }
846 
847  const int* globalCell() const
848  {
849  assert( grid_.global_cell != 0 );
850  return grid_.global_cell;
851  }
852 
853  const int* globalCellPtr() const
854  {
855  return grid_.global_cell;
856  }
857 
858  void getIJK(const int c, std::array<int,3>& ijk) const
859  {
860  int gc = globalCell()[c];
861  ijk[0] = gc % logicalCartesianSize()[0]; gc /= logicalCartesianSize()[0];
862  ijk[1] = gc % logicalCartesianSize()[1];
863  ijk[2] = gc / logicalCartesianSize()[1];
864  }
865 
870 
871  template<class DataHandle>
881  void scatterData([[maybe_unused]] DataHandle& handle) const
882  {
883  OPM_THROW(std::runtime_error, "ScatterData not implemented for polyhedral grid!");
884  }
885 
886  protected:
887 #if HAVE_ECL_INPUT
888  UnstructuredGridType* createGrid( const Opm::EclipseGrid& inputGrid, const std::vector< double >& poreVolumes ) const
889  {
890  struct grdecl g;
891 
892  g.dims[0] = inputGrid.getNX();
893  g.dims[1] = inputGrid.getNY();
894  g.dims[2] = inputGrid.getNZ();
895 
896  std::vector<double> coord = inputGrid.getCOORD( );
897  std::vector<double> zcorn = inputGrid.getZCORN( );
898  std::vector<int> actnum = inputGrid.getACTNUM( );
899 
900  g.coord = coord.data();
901  g.zcorn = zcorn.data();
902  g.actnum = actnum.data();
903 
904  if (!poreVolumes.empty() && (inputGrid.getMinpvMode() != Opm::MinpvMode::ModeEnum::Inactive))
905  {
906  Opm::MinpvProcessor mp(g.dims[0], g.dims[1], g.dims[2]);
907  const std::vector<double>& minpvv = inputGrid.getMinpvVector();
908  // Currently the pinchProcessor is not used and only opmfil is supported
909  // The polyhedralgrid only only supports the opmfil option
910  //bool opmfil = inputGrid.getMinpvMode() == Opm::MinpvMode::OpmFIL;
911  bool opmfil = true;
912  const size_t cartGridSize = g.dims[0] * g.dims[1] * g.dims[2];
913  std::vector<double> thickness(cartGridSize);
914  for (size_t i = 0; i < cartGridSize; ++i) {
915  thickness[i] = inputGrid.getCellThickness(i);
916  }
917  const double z_tolerance = inputGrid.isPinchActive() ? inputGrid.getPinchThresholdThickness() : 0.0;
918  mp.process(thickness, z_tolerance, poreVolumes, minpvv, actnum, opmfil, zcorn.data());
919  }
920 
921  /*
922  if (!poreVolumes.empty() && (inputGrid.getMinpvMode() != Opm::MinpvMode::ModeEnum::Inactive)) {
923  Opm::MinpvProcessor mp(g.dims[0], g.dims[1], g.dims[2]);
924  const double minpv_value = inputGrid.getMinpvValue();
925  // Currently the pinchProcessor is not used and only opmfil is supported
926  //bool opmfil = inputGrid.getMinpvMode() == Opm::MinpvMode::OpmFIL;
927  bool opmfil = true;
928  mp.process(poreVolumes, minpv_value, actnum, opmfil, zcorn.data());
929  }
930  */
931 
932  const double z_tolerance = inputGrid.isPinchActive() ?
933  inputGrid.getPinchThresholdThickness() : 0.0;
934  UnstructuredGridType* cgrid = create_grid_cornerpoint(&g, z_tolerance);
935  if (!cgrid) {
936  OPM_THROW(std::runtime_error, "Failed to construct grid.");
937  }
938  return cgrid;
939  }
940 #endif
941 
942  UnstructuredGridType* createGrid( const std::vector< int >& n, const std::vector< double >& dx ) const
943  {
944  UnstructuredGridType* cgrid = nullptr ;
945  assert( int(n.size()) == dim );
946  if( dim == 2 )
947  {
948  cgrid = create_grid_cart2d( n[ 0 ], n[ 1 ], dx[ 0 ], dx[ 1 ] );
949  }
950  else if ( dim == 3 )
951  {
952  cgrid = create_grid_hexa3d( n[ 0 ], n[ 1 ], n[ 2 ], dx[ 0 ], dx[ 1 ], dx[ 2 ] );
953  }
954 
955  //print_grid( cgrid );
956  if (!cgrid) {
957  OPM_THROW(std::runtime_error, "Failed to construct grid.");
958  }
959  return cgrid;
960  }
961 
962  public:
963 #if DUNE_VERSION_LT_REV(DUNE_GRID, 2, 7, 1)
964  using Base::getRealImplementation;
965 #endif
966  typedef typename Traits :: ExtraData ExtraData;
967  ExtraData extraData () const { return this; }
968 
969  template <class EntitySeed>
970  int corners( const EntitySeed& seed ) const
971  {
972  const int codim = EntitySeed :: codimension;
973  const int index = seed.index();
974  if (codim==0)
975  return cellVertices_[ index ].size();
976  if (codim==1)
977  return grid_.face_nodepos[ index+1 ] - grid_.face_nodepos[ index ];
978  if (codim==dim)
979  return 1;
980  return 0;
981  }
982 
983  template <class EntitySeed>
984  GlobalCoordinate
985  corner ( const EntitySeed& seed, const int i ) const
986  {
987  const int codim = EntitySeed :: codimension;
988  if (codim==0)
989  {
990  const int coordIndex = GlobalCoordinate :: dimension * cellVertices_[ seed.index() ][ i ];
991  return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex );
992  }
993  if (codim==1)
994  {
995  // for faces we need to swap vertices in 3d since in UnstructuredGrid
996  // those are ordered counter clockwise, for 2d this does not matter
997  // TODO: Improve this for performance reasons
998  const int crners = corners( seed );
999  const int crner = (crners == 4 && EntitySeed :: dimension == 3 && i > 1 ) ? 5 - i : i;
1000  const int faceVertex = grid_.face_nodes[ grid_.face_nodepos[seed.index() ] + crner ];
1001  return copyToGlobalCoordinate( grid_.node_coordinates + GlobalCoordinate :: dimension * faceVertex );
1002  }
1003  if (codim==dim)
1004  {
1005  const int coordIndex = GlobalCoordinate :: dimension * seed.index();
1006  return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex );
1007  }
1008  return GlobalCoordinate( 0 );
1009  }
1010 
1011  template <class EntitySeed>
1012  int subEntities( const EntitySeed& seed, const int codim ) const
1013  {
1014  const int index = seed.index();
1015  if( seed.codimension == 0 )
1016  {
1017  if (codim==0)
1018  return 1;
1019  if (codim==1)
1020  return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ];
1021  if (codim==dim)
1022  return cellVertices_[ index ].size();
1023  }
1024  else if( seed.codimension == 1 )
1025  {
1026  if (codim==1)
1027  return 1;
1028  if (codim==dim)
1029  return grid_.face_nodepos[ index+1 ] - grid_.face_nodepos[ index ];
1030  }
1031  else if ( seed.codimension == dim )
1032  {
1033  return 1 ;
1034  }
1035 
1036  return 0;
1037  }
1038 
1039  template <int codim, class EntitySeedArg >
1040  typename Codim<codim>::EntitySeed
1041  subEntitySeed( const EntitySeedArg& baseSeed, const int i ) const
1042  {
1043  assert( codim >= EntitySeedArg::codimension );
1044  assert( i>= 0 && i<subEntities( baseSeed, codim ) );
1045  typedef typename Codim<codim>::EntitySeed EntitySeed;
1046 
1047  // if codim equals entity seed codim just return same entity seed.
1048  if( codim == EntitySeedArg::codimension )
1049  {
1050  return EntitySeed( baseSeed.index() );
1051  }
1052 
1053  if( EntitySeedArg::codimension == 0 )
1054  {
1055  if ( codim == 1 )
1056  {
1057  return EntitySeed( grid_.cell_faces[ grid_.cell_facepos[ baseSeed.index() ] + i ] );
1058  }
1059  else if ( codim == dim )
1060  {
1061  return EntitySeed( cellVertices_[ baseSeed.index() ][ i ] );
1062  }
1063  }
1064  else if ( EntitySeedArg::codimension == 1 && codim == dim )
1065  {
1066  return EntitySeed( grid_.face_nodes[ grid_.face_nodepos[ baseSeed.index() + i ] ]);
1067  }
1068 
1069  DUNE_THROW(NotImplemented,"codimension not available");
1070  return EntitySeed();
1071  }
1072 
1073  template <int codim>
1074  typename Codim<codim>::EntitySeed
1075  subEntitySeed( const typename Codim<1>::EntitySeed& faceSeed, const int i ) const
1076  {
1077  assert( i>= 0 && i<subEntities( faceSeed, codim ) );
1078  typedef typename Codim<codim>::EntitySeed EntitySeed;
1079  if ( codim == 1 )
1080  {
1081  return EntitySeed( faceSeed.index() );
1082  }
1083  else if ( codim == dim )
1084  {
1085  return EntitySeed( grid_.face_nodes[ grid_.face_nodepos[ faceSeed.index() ] + i ] );
1086  }
1087  else
1088  {
1089  DUNE_THROW(NotImplemented,"codimension not available");
1090  }
1091  }
1092 
1093  bool hasBoundaryIntersections(const typename Codim<0>::EntitySeed& seed ) const
1094  {
1095  const int faces = subEntities( seed, 1 );
1096  for( int f=0; f<faces; ++f )
1097  {
1098  const auto faceSeed = this->template subEntitySeed<1>( seed, f );
1099  if( isBoundaryFace( faceSeed ) )
1100  return true;
1101  }
1102  return false;
1103  }
1104 
1105  bool isBoundaryFace(const int face ) const
1106  {
1107  assert( face >= 0 && face < grid_.number_of_faces );
1108  const int facePos = 2 * face;
1109  return ((grid_.face_cells[ facePos ] < 0) || (grid_.face_cells[ facePos+1 ] < 0));
1110  }
1111 
1112  bool isBoundaryFace(const typename Codim<1>::EntitySeed& faceSeed ) const
1113  {
1114  assert( faceSeed.isValid() );
1115  return isBoundaryFace( faceSeed.index() );
1116  }
1117 
1118  int boundarySegmentIndex(const typename Codim<0>::EntitySeed& seed, const int face ) const
1119  {
1120  const auto faceSeed = this->template subEntitySeed<1>( seed, face );
1121  assert( faceSeed.isValid() );
1122  const int facePos = 2 * faceSeed.index();
1123  const int idx = std::min( grid_.face_cells[ facePos ], grid_.face_cells[ facePos+1 ]);
1124  // check that this is actually the boundary
1125  assert( idx < 0 );
1126  return -(idx+1); // +1 to include 0 boundary segment index
1127  }
1128 
1129  const std::vector< GeometryType > &geomTypes ( const unsigned int codim ) const
1130  {
1131  static std::vector< GeometryType > emptyDummy;
1132  if (codim < geomTypes_.size())
1133  {
1134  return geomTypes_[codim];
1135  }
1136 
1137  return emptyDummy;
1138  }
1139 
1140  template < class Seed >
1141  GeometryType geometryType( const Seed& seed ) const
1142  {
1143  if( Seed::codimension == 0 )
1144  {
1145  if( cellGeomTypes_.empty() )
1146  {
1147  assert( geomTypes( Seed::codimension ).size() == 1 );
1148  return geomTypes( Seed::codimension )[ 0 ];
1149  }
1150  else
1151  {
1152  assert( static_cast<size_t>(seed.index()) < cellGeomTypes_.size() );
1153  return cellGeomTypes_[ seed.index() ];
1154  }
1155  }
1156  else
1157  {
1158  // 3d faces
1159  if( dim == 3 && Seed::codimension == 1 )
1160  {
1161  GeometryType face;
1162  const int nVx = corners( seed );
1163  if( nVx == 4 ) // quad face
1164  face = Dune::GeometryTypes::cube(2);
1165  else if( nVx == 3 ) // triangle face
1166  face = Dune::GeometryTypes::simplex(2);
1167  else // polygonal face
1168  face = Dune::GeometryTypes::none(2);
1169  return face;
1170  }
1171 
1172  // for codim 2 and codim 3 there is only one option
1173  return geomTypes( Seed::codimension )[ 0 ];
1174  }
1175  }
1176 
1177  int indexInInside( const typename Codim<0>::EntitySeed& seed, const int i ) const
1178  {
1179  return ( grid_.cell_facetag ) ? cartesianIndexInInside( seed, i ) : i;
1180  }
1181 
1182  int cartesianIndexInInside( const typename Codim<0>::EntitySeed& seed, const int i ) const
1183  {
1184  assert( i>= 0 && i<subEntities( seed, 1 ) );
1185  return grid_.cell_facetag[ grid_.cell_facepos[ seed.index() ] + i ] ;
1186  }
1187 
1188  typename Codim<0>::EntitySeed
1189  neighbor( const typename Codim<0>::EntitySeed& seed, const int i ) const
1190  {
1191  const int face = this->template subEntitySeed<1>( seed, i ).index();
1192  int nb = grid_.face_cells[ 2 * face ];
1193  if( nb == seed.index() )
1194  {
1195  nb = grid_.face_cells[ 2 * face + 1 ];
1196  }
1197 
1198  typedef typename Codim<0>::EntitySeed EntitySeed;
1199  return EntitySeed( nb );
1200  }
1201 
1202  int
1203  indexInOutside( const typename Codim<0>::EntitySeed& seed, const int i ) const
1204  {
1205  if( grid_.cell_facetag )
1206  {
1207  // if cell_facetag is present we assume pseudo Cartesian corner point case
1208  const int in_inside = cartesianIndexInInside( seed, i );
1209  return in_inside + ((in_inside % 2) ? -1 : 1);
1210  }
1211  else
1212  {
1213  typedef typename Codim<0>::EntitySeed EntitySeed;
1214  EntitySeed nb = neighbor( seed, i );
1215  const int faces = subEntities( seed, 1 );
1216  for( int face = 0; face<faces; ++ face )
1217  {
1218  if( neighbor( nb, face ).equals(seed) )
1219  {
1220  return indexInInside( nb, face );
1221  }
1222  }
1223  DUNE_THROW(InvalidStateException,"inverse intersection not found");
1224  return -1;
1225  }
1226  }
1227 
1228  template <class EntitySeed>
1229  GlobalCoordinate
1230  outerNormal( const EntitySeed& seed, const int i ) const
1231  {
1232  const int face = this->template subEntitySeed<1>( seed, i ).index();
1233  const int normalIdx = face * GlobalCoordinate :: dimension ;
1234  GlobalCoordinate normal = copyToGlobalCoordinate( grid_.face_normals + normalIdx );
1235  const int nb = grid_.face_cells[ 2*face ];
1236  if( nb != seed.index() )
1237  {
1238  normal *= -1.0;
1239  }
1240  return normal;
1241  }
1242 
1243  template <class EntitySeed>
1244  GlobalCoordinate
1245  unitOuterNormal( const EntitySeed& seed, const int i ) const
1246  {
1247  const int face = this->template subEntitySeed<1>( seed, i ).index();
1248  if( seed.index() == grid_.face_cells[ 2*face ] )
1249  {
1250  return unitOuterNormals_[ face ];
1251  }
1252  else
1253  {
1254  GlobalCoordinate normal = unitOuterNormals_[ face ];
1255  normal *= -1.0;
1256  return normal;
1257  }
1258  }
1259 
1260  template <class EntitySeed>
1261  GlobalCoordinate centroids( const EntitySeed& seed ) const
1262  {
1263  if( ! seed.isValid() )
1264  return GlobalCoordinate( 0 );
1265 
1266  const int index = GlobalCoordinate :: dimension * seed.index();
1267  const int codim = EntitySeed::codimension;
1268  assert( index >= 0 && index < size( codim ) * GlobalCoordinate :: dimension );
1269 
1270  if( codim == 0 )
1271  {
1272  return copyToGlobalCoordinate( grid_.cell_centroids + index );
1273  }
1274  else if ( codim == 1 )
1275  {
1276  return copyToGlobalCoordinate( grid_.face_centroids + index );
1277  }
1278  else if( codim == dim )
1279  {
1280  return copyToGlobalCoordinate( grid_.node_coordinates + index );
1281  }
1282  else
1283  {
1284  DUNE_THROW(InvalidStateException,"codimension not implemented");
1285  return GlobalCoordinate( 0 );
1286  }
1287  }
1288 
1289  GlobalCoordinate copyToGlobalCoordinate( const double* coords ) const
1290  {
1291  GlobalCoordinate coordinate;
1292  for( int i=0; i<GlobalCoordinate::dimension; ++i )
1293  {
1294  coordinate[ i ] = coords[ i ];
1295  }
1296  return coordinate;
1297  }
1298 
1299  template <class EntitySeed>
1300  double volumes( const EntitySeed& seed ) const
1301  {
1302  static const int codim = EntitySeed::codimension;
1303  if( codim == dim || ! seed.isValid() )
1304  {
1305  return 1.0;
1306  }
1307  else
1308  {
1309  assert( seed.isValid() );
1310 
1311  if( codim == 0 )
1312  {
1313  return grid_.cell_volumes[ seed.index() ];
1314  }
1315  else if ( codim == 1 )
1316  {
1317  return grid_.face_areas[ seed.index() ];
1318  }
1319  else
1320  {
1321  DUNE_THROW(InvalidStateException,"codimension not implemented");
1322  return 0.0;
1323  }
1324  }
1325  }
1326 
1327  protected:
1328  void init ()
1329  {
1330  // copy Cartesian dimensions
1331  for( int i=0; i<3; ++i )
1332  {
1333  cartDims_[ i ] = grid_.cartdims[ i ];
1334  }
1335 
1336  cellGeomTypes_.clear();
1337 
1338  // setup list of cell vertices
1339  const int numCells = size( 0 );
1340 
1341  cellVertices_.resize( numCells );
1342 
1343  // sort vertices such that they comply with the dune reference cube
1344  if( grid_.cell_facetag )
1345  {
1346  typedef std::array<int, 3> KeyType;
1347  std::map< const KeyType, const int > vertexFaceTags;
1348  const int vertexFacePattern [8][3] = {
1349  { 0, 2, 4 }, // vertex 0
1350  { 1, 2, 4 }, // vertex 1
1351  { 0, 3, 4 }, // vertex 2
1352  { 1, 3, 4 }, // vertex 3
1353  { 0, 2, 5 }, // vertex 4
1354  { 1, 2, 5 }, // vertex 5
1355  { 0, 3, 5 }, // vertex 6
1356  { 1, 3, 5 } // vertex 7
1357  };
1358 
1359  for( int i=0; i<8; ++i )
1360  {
1361  KeyType key; key.fill( 4 ); // default is 4 which is the first z coord (for the 2d case)
1362  for( int j=0; j<dim; ++j )
1363  {
1364  key[ j ] = vertexFacePattern[ i ][ j ];
1365  }
1366 
1367  vertexFaceTags.insert( std::make_pair( key, i ) );
1368  }
1369 
1370  for (int c = 0; c < numCells; ++c)
1371  {
1372  if( dim == 2 )
1373  {
1374  // for 2d Cartesian grids the face ordering is wrong
1375  int f = grid_.cell_facepos[ c ];
1376  std::swap( grid_.cell_faces[ f+1 ], grid_.cell_faces[ f+2 ] );
1377  std::swap( grid_.cell_facetag[ f+1 ], grid_.cell_facetag[ f+2 ] );
1378  }
1379 
1380  typedef std::map<int,int> vertexmap_t;
1381  typedef typename vertexmap_t :: iterator iterator;
1382 
1383  std::vector< vertexmap_t > cell_pts( dim*2 );
1384 
1385  for (int hf=grid_.cell_facepos[ c ]; hf < grid_.cell_facepos[c+1]; ++hf)
1386  {
1387  const int f = grid_.cell_faces[ hf ];
1388  const int faceTag = grid_.cell_facetag[ hf ];
1389 
1390  for( int nodepos=grid_.face_nodepos[f]; nodepos<grid_.face_nodepos[f+1]; ++nodepos )
1391  {
1392  const int node = grid_.face_nodes[ nodepos ];
1393  iterator it = cell_pts[ faceTag ].find( node );
1394  if( it == cell_pts[ faceTag ].end() )
1395  {
1396  cell_pts[ faceTag ].insert( std::make_pair( node, 1 ) );
1397  }
1398  else
1399  {
1400  // increase vertex reference counter
1401  (*it).second++;
1402  }
1403  }
1404  }
1405 
1406  typedef std::map< int, std::set<int> > vertexlist_t;
1407  vertexlist_t vertexList;
1408 
1409  for( int faceTag = 0; faceTag<dim*2; ++faceTag )
1410  {
1411  for( iterator it = cell_pts[ faceTag ].begin(),
1412  end = cell_pts[ faceTag ].end(); it != end; ++it )
1413  {
1414 
1415  // only consider vertices with one appearance
1416  if( (*it).second == 1 )
1417  {
1418  vertexList[ (*it).first ].insert( faceTag );
1419  }
1420  }
1421  }
1422 
1423  assert( int(vertexList.size()) == ( dim == 2 ? 4 : 8) );
1424 
1425  cellVertices_[ c ].resize( vertexList.size() );
1426  for( auto it = vertexList.begin(), end = vertexList.end(); it != end; ++it )
1427  {
1428  assert( (*it).second.size() == dim );
1429  KeyType key; key.fill( 4 ); // fill with 4 which is the first z coord
1430 
1431  std::copy( (*it).second.begin(), (*it).second.end(), key.begin() );
1432  auto vx = vertexFaceTags.find( key );
1433  assert( vx != vertexFaceTags.end() );
1434  if( vx != vertexFaceTags.end() )
1435  {
1436  if( (*vx).second >= int(cellVertices_[ c ].size()) )
1437  cellVertices_[ c ].resize( (*vx).second+1 );
1438  // store node number on correct local position
1439  cellVertices_[ c ][ (*vx).second ] = (*it).first ;
1440  }
1441  }
1442  }
1443 
1444  // if face_tag is available we assume that the elements follow a cube-like structure
1445  geomTypes_.resize(dim + 1);
1446  GeometryType tmp;
1447  for (int codim = 0; codim <= dim; ++codim)
1448  {
1449  tmp = Dune::GeometryTypes::cube(dim - codim);
1450  geomTypes_[codim].push_back(tmp);
1451  }
1452  }
1453  else // if ( grid_.cell_facetag )
1454  {
1455  int maxVx = 0 ;
1456  int minVx = std::numeric_limits<int>::max();
1457 
1458  for (int c = 0; c < numCells; ++c)
1459  {
1460  std::set<int> cell_pts;
1461  for (int hf=grid_.cell_facepos[ c ]; hf < grid_.cell_facepos[c+1]; ++hf)
1462  {
1463  int f = grid_.cell_faces[ hf ];
1464  const int* fnbeg = grid_.face_nodes + grid_.face_nodepos[f];
1465  const int* fnend = grid_.face_nodes + grid_.face_nodepos[f+1];
1466  cell_pts.insert(fnbeg, fnend);
1467  }
1468 
1469  cellVertices_[ c ].resize( cell_pts.size() );
1470  std::copy(cell_pts.begin(), cell_pts.end(), cellVertices_[ c ].begin() );
1471  maxVx = std::max( maxVx, int( cell_pts.size() ) );
1472  minVx = std::min( minVx, int( cell_pts.size() ) );
1473  }
1474 
1475  if( minVx == maxVx && maxVx == 4 )
1476  {
1477  for (int c = 0; c < numCells; ++c)
1478  {
1479  assert( cellVertices_[ c ].size() == 4 );
1480  GlobalCoordinate center( 0 );
1481  GlobalCoordinate p[ dim+1 ];
1482  for( int i=0; i<dim+1; ++i )
1483  {
1484  const int vertex = cellVertices_[ c ][ i ];
1485 
1486  for( int d=0; d<dim; ++d )
1487  {
1488  center[ d ] += grid_.node_coordinates[ vertex*dim + d ];
1489  p[ i ][ d ] = grid_.node_coordinates[ vertex*dim + d ];
1490  }
1491  }
1492  center *= 0.25;
1493  for( int d=0; d<dim; ++d )
1494  {
1495  grid_.cell_centroids[ c*dim + d ] = center[ d ];
1496  }
1497 
1498  Dune::GeometryType simplex;
1499  simplex = Dune::GeometryTypes::simplex(dim);
1500 
1501  typedef Dune::AffineGeometry< ctype, dim, dimworld> AffineGeometryType;
1502  AffineGeometryType geometry( simplex, p );
1503  grid_.cell_volumes[ c ] = geometry.volume();
1504  }
1505  }
1506 
1507  // check face normals
1508  {
1509  const int faces = grid_.number_of_faces;
1510  for( int face = 0 ; face < faces; ++face )
1511  {
1512  const int a = grid_.face_cells[ 2*face ];
1513  const int b = grid_.face_cells[ 2*face + 1 ];
1514 
1515  assert( a >=0 || b >=0 );
1516 
1517  if( grid_.face_areas[ face ] < 0 )
1518  std::abort();
1519 
1520  GlobalCoordinate centerDiff( 0 );
1521  if( b >= 0 )
1522  {
1523  for( int d=0; d<dimworld; ++d )
1524  {
1525  centerDiff[ d ] = grid_.cell_centroids[ b*dimworld + d ];
1526  }
1527  }
1528  else
1529  {
1530  for( int d=0; d<dimworld; ++d )
1531  {
1532  centerDiff[ d ] = grid_.face_centroids[ face*dimworld + d ];
1533  }
1534  }
1535 
1536  if( a >= 0 )
1537  {
1538  for( int d=0; d<dimworld; ++d )
1539  {
1540  centerDiff[ d ] -= grid_.cell_centroids[ a*dimworld + d ];
1541  }
1542  }
1543  else
1544  {
1545  for( int d=0; d<dimworld; ++d )
1546  {
1547  centerDiff[ d ] -= grid_.face_centroids[ face*dimworld + d ];
1548  }
1549  }
1550 
1551  GlobalCoordinate normal( 0 );
1552  for( int d=0; d<dimworld; ++d )
1553  {
1554  normal[ d ] = grid_.face_normals[ face*dimworld + d ];
1555  }
1556 
1557  if( centerDiff.two_norm() < 1e-10 )
1558  std::abort();
1559 
1560  // if diff and normal point in different direction, flip faces
1561  if( centerDiff * normal < 0 )
1562  {
1563  grid_.face_cells[ 2*face ] = b;
1564  grid_.face_cells[ 2*face + 1 ] = a;
1565  }
1566  }
1567  }
1568 
1569  // if no face_tag is available we set the reference element based
1570  // on the number of nodes in a cell.
1571  // By default set all types to None. This corresponds to hasPolygon
1572  GeometryType tmp;
1573  tmp = Dune::GeometryTypes::none(dim);
1574  cellGeomTypes_.resize( numCells );
1575  std::fill( cellGeomTypes_.begin(), cellGeomTypes_.end(), tmp );
1576 
1577  bool hasSimplex = false ;
1578  bool hasCube = false ;
1579  bool hasPolyhedron = false;
1580 
1581  for (int c = 0; c < numCells; ++c)
1582  {
1583  const int nVx = cellVertices_[ c ].size();
1584  if( nVx == 4 )
1585  {
1586  cellGeomTypes_[ c ] = Dune::GeometryTypes::simplex(dim);
1587  hasSimplex = true;
1588  }
1589  else if( nVx == 8 )
1590  {
1591  cellGeomTypes_[ c ] = Dune::GeometryTypes::cube(dim);
1592 
1593  hasCube = true;
1594  }
1595  else
1596  {
1597  hasPolyhedron = true;
1598  }
1599  }
1600  // Propogate the cell geometry type to all codimensions
1601  geomTypes_.resize(dim + 1);
1602  for (int codim = 0; codim <= dim; ++codim)
1603  {
1604  if( hasSimplex )
1605  {
1606  tmp = Dune::GeometryTypes::simplex(dim - codim);
1607  geomTypes_[ codim ].push_back( tmp );
1608  }
1609  else if ( hasCube )
1610  {
1611  tmp = Dune::GeometryTypes::cube(dim - codim);
1612  geomTypes_[ codim ].push_back( tmp );
1613  }
1614  else if (hasPolyhedron)
1615  {
1616  tmp = Dune::GeometryTypes::none(dim - codim);
1617  geomTypes_[ codim ].push_back( tmp );
1618  }
1619  else
1620  {
1621  OPM_THROW(std::runtime_error, "Grid error, unkown geometry type.");
1622  }
1623  }
1624 
1625  } // end else of ( grid_.cell_facetag )
1626 
1627  nBndSegments_ = 0;
1628  unitOuterNormals_.resize( grid_.number_of_faces );
1629  for( int face = 0; face < grid_.number_of_faces; ++face )
1630  {
1631  const int normalIdx = face * GlobalCoordinate :: dimension ;
1632  GlobalCoordinate normal = copyToGlobalCoordinate( grid_.face_normals + normalIdx );
1633  normal /= normal.two_norm();
1634  unitOuterNormals_[ face ] = normal;
1635 
1636  if( isBoundaryFace( face ) )
1637  {
1638  // increase number if boundary segments
1639  ++nBndSegments_;
1640  const int facePos = 2 * face ;
1641  // store negative number to indicate boundary
1642  // the abstract value is the segment index
1643  if( grid_.face_cells[ facePos ] < 0 )
1644  {
1645  grid_.face_cells[ facePos ] = -nBndSegments_;
1646  }
1647  else if ( grid_.face_cells[ facePos+1 ] < 0 )
1648  {
1649  grid_.face_cells[ facePos+1 ] = -nBndSegments_;
1650  }
1651  }
1652  }
1653  }
1654 
1655  void print( std::ostream& out, const UnstructuredGridType& grid ) const
1656  {
1657  const int numCells = grid.number_of_cells;
1658  for( int c=0; c<numCells; ++c )
1659  {
1660  out << "cell " << c << " : faces = " << std::endl;
1661  for (int hf=grid.cell_facepos[ c ]; hf < grid.cell_facepos[c+1]; ++hf)
1662  {
1663  int f = grid_.cell_faces[ hf ];
1664  const int* fnbeg = grid_.face_nodes + grid_.face_nodepos[f];
1665  const int* fnend = grid_.face_nodes + grid_.face_nodepos[f+1];
1666  out << f << " vx = " ;
1667  while( fnbeg != fnend )
1668  {
1669  out << *fnbeg << " ";
1670  ++fnbeg;
1671  }
1672  out << std::endl;
1673  }
1674  out << std::endl;
1675 
1676  const auto& vx = cellVertices_[ c ];
1677  out << "cell " << c << " : vertices = ";
1678  for( size_t i=0; i<vx.size(); ++i )
1679  out << vx[ i ] << " ";
1680  out << std::endl;
1681  }
1682 
1683  }
1684 
1685  protected:
1686  UnstructuredGridPtr gridPtr_;
1687  const UnstructuredGridType& grid_;
1688 
1690  std::array< int, 3 > cartDims_;
1691  std::vector< std::vector< GeometryType > > geomTypes_;
1692  std::vector< std::vector< int > > cellVertices_;
1693 
1694  std::vector< GlobalCoordinate > unitOuterNormals_;
1695 
1696  // geometry type of each cell if existing (if not then all are polyhedral)
1697  std::vector< GeometryType > cellGeomTypes_;
1698 
1699  mutable LeafIndexSet leafIndexSet_;
1700  mutable GlobalIdSet globalIdSet_;
1701  mutable LocalIdSet localIdSet_;
1702 
1703  size_t nBndSegments_;
1704 
1705  private:
1706  // no copying
1707  PolyhedralGrid ( const PolyhedralGrid& );
1708  };
1709 
1710 
1711 
1712  // PolyhedralGrid::Codim
1713  // -------------
1714 
1715  template< int dim, int dimworld, typename coord_t >
1716  template< int codim >
1718  : public Base::template Codim< codim >
1719  {
1728  typedef typename Traits::template Codim< codim >::Entity Entity;
1729 
1734  typedef typename Traits::template Codim< codim >::EntityPointer EntityPointer;
1735 
1749  typedef typename Traits::template Codim< codim >::Geometry Geometry;
1750 
1759  typedef typename Traits::template Codim< codim >::LocalGeometry LocalGeometry;
1760 
1766  template< PartitionIteratorType pitype >
1767  struct Partition
1768  {
1769  typedef typename Traits::template Codim< codim >
1770  ::template Partition< pitype >::LeafIterator
1771  LeafIterator;
1772  typedef typename Traits::template Codim< codim >
1773  ::template Partition< pitype >::LevelIterator
1774  LevelIterator;
1775  };
1776 
1784  typedef typename Partition< All_Partition >::LeafIterator LeafIterator;
1785 
1793  typedef typename Partition< All_Partition >::LevelIterator LevelIterator;
1794 
1796  };
1797 
1798 } // namespace Dune
1799 
1800 #include <opm/grid/polyhedralgrid/persistentcontainer.hh>
1801 #include <opm/grid/polyhedralgrid/cartesianindexmapper.hh>
1802 #include <opm/grid/polyhedralgrid/gridhelpers.hh>
1803 
1804 #endif // #ifndef DUNE_POLYHEDRALGRID_GRID_HH
Main OPM-Core grid data structure along with helper functions for construction, destruction and readi...
void destroy_grid(struct UnstructuredGrid *g)
Destroy and deallocate an UnstructuredGrid and all its data.
Definition: UnstructuredGrid.c:32
struct UnstructuredGrid * allocate_grid(size_t ndims, size_t ncells, size_t nfaces, size_t nfacenodes, size_t ncellfaces, size_t nnodes)
Allocate and initialise an UnstructuredGrid where pointers are set to location with correct size.
Definition: UnstructuredGrid.c:87
Routines to construct fully formed grid structures from a simple Cartesian (i.e., tensor product) des...
struct UnstructuredGrid * create_grid_cart2d(int nx, int ny, double dx, double dy)
Form geometrically Cartesian grid in two space dimensions with equally sized cells.
Definition: cart_grid.c:95
struct UnstructuredGrid * create_grid_hexa3d(int nx, int ny, int nz, double dx, double dy, double dz)
Form geometrically Cartesian grid in three space dimensions with equally sized cells.
Definition: cart_grid.c:59
Definition: entityseed.hh:16
Definition: entity.hh:152
Definition: geometry.hh:230
Definition: idset.hh:17
Definition: indexset.hh:24
Definition: intersectioniterator.hh:16
Definition: intersection.hh:20
Definition: iterator.hh:21
Definition: geometry.hh:249
identical grid wrapper
Definition: grid.hh:158
void postAdapt()
Definition: grid.hh:602
Traits::ctype ctype
type of vector coordinates (e.g., double)
Definition: grid.hh:307
int ghostSize(int, int codim) const
obtain size of ghost region for a grid level
Definition: grid.hh:642
bool loadBalance(CommDataHandleIF< DataHandle, Data > &)
rebalance the load each process has to handle
Definition: grid.hh:746
const CollectiveCommunication & comm() const
obtain CollectiveCommunication object
Definition: grid.hh:709
Traits::template Codim< EntitySeed::codimension >::EntityPointer entityPointer(const EntitySeed &seed) const
obtain EntityPointer from EntitySeed.
Definition: grid.hh:806
Traits::template Codim< EntitySeed::codimension >::Entity entity(const EntitySeed &seed) const
obtain EntityPointer from EntitySeed.
Definition: grid.hh:817
bool preAdapt()
Definition: grid.hh:580
int ghostSize(int codim) const
obtain size of ghost region for the leaf grid
Definition: grid.hh:622
bool loadBalance(DofManager &)
rebalance the load each process has to handle
Definition: grid.hh:766
Traits::GlobalIdSet GlobalIdSet
type of global id set
Definition: grid.hh:282
int maxLevel() const
obtain maximal grid level
Definition: grid.hh:412
Traits::LevelIndexSet LevelIndexSet
type of level index set
Definition: grid.hh:270
LevelGridView levelGridView(int) const
View for a grid level for All_Partition.
Definition: grid.hh:790
PolyhedralGrid(const std::vector< int > &n, const std::vector< double > &dx)
constructor
Definition: grid.hh:344
Partition< pitype >::LeafGridView leafGridView() const
View for the leaf grid.
Definition: grid.hh:782
Traits::HierarchicIterator HierarchicIterator
iterator over the grid hierarchy
Definition: grid.hh:224
void scatterData([[maybe_unused]] DataHandle &handle) const
Moves data from the global (all data on process) view to the distributed view.
Definition: grid.hh:881
int size(int, int codim) const
obtain number of entites on a level
Definition: grid.hh:425
Partition< pitype >::LevelGridView levelGridView(int) const
View for a grid level.
Definition: grid.hh:773
Partition< All_Partition >::LevelGridView LevelGridView
View types for All_Partition.
Definition: grid.hh:244
void update()
update grid caches
Definition: grid.hh:836
void switchToDistributedView()
Switch to the distributed view.
Definition: grid.hh:696
int overlapSize(int, int) const
obtain size of overlap region for a grid level
Definition: grid.hh:632
int size(int codim) const
obtain number of leaf entities
Definition: grid.hh:436
Traits::CollectiveCommunication CollectiveCommunication
communicator with all other processes having some part of the grid
Definition: grid.hh:310
Traits::LocalIdSet LocalIdSet
type of local id set
Definition: grid.hh:299
int overlapSize(int) const
obtain size of overlap region for the leaf grid
Definition: grid.hh:613
Traits::LevelIntersectionIterator LevelIntersectionIterator
iterator over intersections with other entities on the same level
Definition: grid.hh:228
GridFamily::Traits Traits
type of the grid traits
Definition: grid.hh:207
PolyhedralGrid(UnstructuredGridPtr &&gridPtr)
constructor
Definition: grid.hh:363
Traits::LeafIntersectionIterator LeafIntersectionIterator
iterator over intersections with other entities on the leaf level
Definition: grid.hh:226
bool adapt(DataHandle &)
Definition: grid.hh:596
int size(int, GeometryType type) const
obtain number of entites on a level
Definition: grid.hh:465
void communicate(DataHandle &, InterfaceType, CommunicationDirection) const
communicate information on leaf entities
Definition: grid.hh:682
size_t numBoundarySegments() const
obtain number of leaf entities
Definition: grid.hh:485
void switchToGlobalView()
Switch to the global view.
Definition: grid.hh:690
bool loadBalance()
rebalance the load each process has to handle
Definition: grid.hh:725
PolyhedralGrid(const UnstructuredGridType &grid)
constructor
Definition: grid.hh:382
LeafGridView leafGridView() const
View for the leaf grid for All_Partition.
Definition: grid.hh:797
int size(GeometryType type) const
returns the number of boundary segments within the macro grid
Definition: grid.hh:474
void communicate(DataHandle &, InterfaceType, CommunicationDirection, int) const
communicate information on a grid level
Definition: grid.hh:661
bool adapt()
Definition: grid.hh:586
Traits::LeafIndexSet LeafIndexSet
type of leaf index set
Definition: grid.hh:260
Transform a corner-point grid ZCORN field to account for MINPV processing.
Definition: MinpvProcessor.hpp:34
Routines to form a complete UnstructuredGrid from a corner-point specification.
struct UnstructuredGrid * create_grid_cornerpoint(const struct grdecl *in, double tol)
Construct grid representation from corner-point specification of a particular geological model.
Definition: cornerpoint_grid.c:164
void compute_geometry(struct UnstructuredGrid *g)
Compute derived geometric primitives in a grid.
Definition: cornerpoint_grid.c:137
Copyright 2019 Equinor AS.
Definition: CartesianIndexMapper.hpp:10
Low-level corner-point processing routines and supporting data structures.
Definition: grid.hh:55
traits structure containing types for a codimension
Definition: grid.hh:1719
Partition< All_Partition >::LevelIterator LevelIterator
type of leaf iterator
Definition: grid.hh:1793
Traits::template Codim< codim >::Entity Entity
type of entity
Definition: grid.hh:1728
Traits::template Codim< codim >::LocalGeometry LocalGeometry
type of local geometry
Definition: grid.hh:1759
Traits::template Codim< codim >::EntityPointer EntityPointer
type of entity pointer
Definition: grid.hh:1734
Partition< All_Partition >::LeafIterator LeafIterator
type of level iterator
Definition: grid.hh:1784
Traits::template Codim< codim >::Geometry Geometry
type of world geometry
Definition: grid.hh:1749
Types for GridView.
Definition: grid.hh:238
Data structure for an unstructured grid, unstructured meaning that any cell may have an arbitrary num...
Definition: UnstructuredGrid.h:99
int * face_nodes
Contains for each face, the indices of its adjacent nodes.
Definition: UnstructuredGrid.h:121
int * face_nodepos
For a face f, face_nodepos[f] contains the starting index for f's nodes in the face_nodes array.
Definition: UnstructuredGrid.h:127
int number_of_cells
The number of cells in the grid.
Definition: UnstructuredGrid.h:109
double * face_areas
Exact or approximate face areas.
Definition: UnstructuredGrid.h:173
int * cell_faces
Contains for each cell, the indices of its adjacent faces.
Definition: UnstructuredGrid.h:146
int number_of_faces
The number of faces in the grid.
Definition: UnstructuredGrid.h:111
double * cell_centroids
Exact or approximate cell centroids, stored consecutively for each cell.
Definition: UnstructuredGrid.h:192
int * cell_facepos
For a cell c, cell_facepos[c] contains the starting index for c's faces in the cell_faces array.
Definition: UnstructuredGrid.h:152
int * cell_facetag
If non-null, this array contains a number for cell-face adjacency indicating the face's position with...
Definition: UnstructuredGrid.h:244
double * cell_volumes
Exact or approximate cell volumes.
Definition: UnstructuredGrid.h:197
int * face_cells
For a face f, face_cells[2*f] and face_cells[2*f + 1] contain the cell indices of the cells adjacent ...
Definition: UnstructuredGrid.h:138
double * face_centroids
Exact or approximate face centroids, stored consecutively for each face.
Definition: UnstructuredGrid.h:168
int number_of_nodes
The number of nodes in the grid.
Definition: UnstructuredGrid.h:113
int cartdims[3]
Contains the size of the logical cartesian structure (if any) of the grid.
Definition: UnstructuredGrid.h:227
int * global_cell
If non-null, this array contains the logical cartesian indices (in a lexicographic ordering) of each ...
Definition: UnstructuredGrid.h:214
double * face_normals
Exact or approximate face normals, stored consecutively for each face.
Definition: UnstructuredGrid.h:184
double * node_coordinates
Node coordinates, stored consecutively for each node.
Definition: UnstructuredGrid.h:160
Raw corner-point specification of a particular geological model.
Definition: preprocess.h:56
const double * coord
Pillar end-points.
Definition: preprocess.h:58
int dims[3]
Cartesian box dimensions.
Definition: preprocess.h:57
const double * zcorn
Corner-point depths.
Definition: preprocess.h:59
const int * actnum
Explicit "active" map.
Definition: preprocess.h:60