All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Groups Pages
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/array.hh>
15 
16 //- dune-grid includes
17 #include <dune/grid/common/grid.hh>
18 #if DUNE_VERSION_NEWER(DUNE_COMMON,2,3)
19 #include <dune/common/parallel/collectivecommunication.hh>
20 #else
21 #include <dune/common/collectivecommunication.hh>
22 #endif
23 
24 //- polyhedralgrid includes
25 #include <dune/grid/polyhedralgrid/capabilities.hh>
26 #include <dune/grid/polyhedralgrid/declaration.hh>
27 #include <dune/grid/polyhedralgrid/entity.hh>
28 #include <dune/grid/polyhedralgrid/entityseed.hh>
29 #include <dune/grid/polyhedralgrid/geometry.hh>
30 #include <dune/grid/polyhedralgrid/gridview.hh>
31 #include <dune/grid/polyhedralgrid/idset.hh>
32 
33 // Re-enable warnings.
34 #include <opm/grid/utility/platform_dependent/reenable_warnings.h>
35 
36 #include <opm/grid/utility/ErrorMacros.hpp>
37 #include <opm/core/grid.h>
39 #include <opm/core/grid/GridManager.hpp>
41 #include <opm/core/grid/MinpvProcessor.hpp>
42 
43 namespace Dune
44 {
45 
46 
47  // PolyhedralGridFamily
48  // ------------
49 
50  template< int dim, int dimworld >
52  {
53  struct Traits
54  {
56 
57  typedef double ctype;
58 
59  // type of data passed to entities, intersections, and iterators
60  // for PolyhedralGrid this is just an empty place holder
61  typedef const Grid* ExtraData;
62 
63  typedef int Index ;
64 
65  static const int dimension = dim;
66  static const int dimensionworld = dimworld;
67 
68  typedef Dune::FieldVector< ctype, dimensionworld > GlobalCoordinate ;
69 
74 
75 #if DUNE_VERSION_NEWER(DUNE_GRID,2,3)
76  typedef Dune::Intersection< const Grid, LeafIntersectionImpl > LeafIntersection;
77  typedef Dune::Intersection< const Grid, LevelIntersectionImpl > LevelIntersection;
78 
79  typedef Dune::IntersectionIterator< const Grid, LeafIntersectionIteratorImpl, LeafIntersectionImpl > LeafIntersectionIterator;
80  typedef Dune::IntersectionIterator< const Grid, LevelIntersectionIteratorImpl, LevelIntersectionImpl > LevelIntersectionIterator;
81 #else
82  typedef Dune::Intersection< const Grid, PolyhedralGridIntersection > LeafIntersection;
83  typedef Dune::Intersection< const Grid, PolyhedralGridIntersection > LevelIntersection;
84 
85  typedef Dune::IntersectionIterator< const Grid, PolyhedralGridIntersectionIterator, PolyhedralGridIntersection > LeafIntersectionIterator;
86  typedef Dune::IntersectionIterator< const Grid, PolyhedralGridIntersectionIterator, PolyhedralGridIntersection > LevelIntersectionIterator;
87 #endif
88 
90  typedef Dune::EntityIterator< 0, const Grid, HierarchicIteratorImpl > HierarchicIterator;
91 
92  template< int codim >
93  struct Codim
94  {
95  typedef PolyhedralGridGeometry<dimension-codim, dimensionworld, const Grid> GeometryImpl;
96  typedef Dune::Geometry< dimension-codim, dimensionworld, const Grid, PolyhedralGridGeometry > Geometry;
97 
98  typedef PolyhedralGridLocalGeometry< dimension-codim, dimensionworld, const Grid> LocalGeometryImpl;
99  typedef Dune::Geometry< dimension-codim, dimension, const Grid, PolyhedralGridLocalGeometry > LocalGeometry;
100 
102  typedef Dune::Entity< codim, dimension, const Grid, PolyhedralGridEntity > Entity;
103 
104 #if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
106  typedef Entity EntityPointer;
107 #else
109  typedef Dune::EntityPointer< const Grid, EntityPointerImpl > EntityPointer;
110 #endif
111 
112  //typedef Dune::EntitySeed< const Grid, PolyhedralGridEntitySeed< codim, const Grid > > EntitySeed;
114 
115  template< PartitionIteratorType pitype >
116  struct Partition
117  {
119  typedef Dune::EntityIterator< codim, const Grid, LeafIteratorImpl > LeafIterator;
120 
121  typedef LeafIterator LevelIterator;
122  };
123 
124  typedef typename Partition< All_Partition >::LeafIterator LeafIterator;
125  typedef typename Partition< All_Partition >::LevelIterator LevelIterator;
126  };
127 
130 
132  typedef GlobalIdSet LocalIdSet;
133 
134  typedef Dune::CollectiveCommunication< Grid > CollectiveCommunication;
135 
136  template< PartitionIteratorType pitype >
137  struct Partition
138  {
139  typedef Dune::GridView< PolyhedralGridViewTraits< dim, dimworld, pitype > > LeafGridView;
140  typedef Dune::GridView< PolyhedralGridViewTraits< dim, dimworld, pitype > > LevelGridView;
141  };
142 
143  typedef typename Partition<All_Partition>::LevelGridView LevelGridView;
144  typedef typename Partition<All_Partition>::LeafGridView LeafGridView;
145 
146  };
147  };
148 
149 
150 
151  // PolyhedralGrid
152  // ------
153 
162  template < int dim, int dimworld >
163  class PolyhedralGrid
165  : public GridDefaultImplementation
166  < dim, dimworld, double, PolyhedralGridFamily< dim, dimworld > >
168  {
169  typedef PolyhedralGrid< dim, dimworld > Grid;
170 
171  typedef GridDefaultImplementation
172  < dim, dimworld, double, PolyhedralGridFamily< dim, dimworld > > Base;
173 
174  typedef UnstructuredGrid UnstructuredGridType;
175 
176  struct UnstructuredGridDeleter
177  {
178  inline void operator () ( UnstructuredGridType* grdPtr )
179  {
180  destroy_grid( grdPtr );
181  }
182  };
183  public:
185  typedef PolyhedralGridFamily< dim, dimworld > GridFamily;
191  typedef typename GridFamily::Traits Traits;
193 
200  template< int codim >
201  struct Codim;
202 
208  typedef typename Traits::HierarchicIterator HierarchicIterator;
211  typedef typename Traits::LeafIntersectionIterator LeafIntersectionIterator;
213  typedef typename Traits::LevelIntersectionIterator LevelIntersectionIterator;
214 
221  template< PartitionIteratorType pitype >
222  struct Partition
223  {
224  typedef typename GridFamily::Traits::template Partition< pitype >::LevelGridView LevelGridView;
225  typedef typename GridFamily::Traits::template Partition< pitype >::LeafGridView LeafGridView;
226  };
227 
229  typedef typename Partition< All_Partition >::LevelGridView LevelGridView;
230  typedef typename Partition< All_Partition >::LeafGridView LeafGridView;
231 
245  typedef typename Traits::LeafIndexSet LeafIndexSet;
246 
255  typedef typename Traits::LevelIndexSet LevelIndexSet;
256 
267  typedef typename Traits::GlobalIdSet GlobalIdSet;
268 
284  typedef typename Traits::LocalIdSet LocalIdSet;
285 
291  typedef typename Traits::ctype ctype;
293 
295  typedef typename Traits::CollectiveCommunication CollectiveCommunication;
296 
297  typedef typename Traits :: GlobalCoordinate GlobalCoordinate;
298 
304 #if HAVE_OPM_PARSER
305 
310  explicit PolyhedralGrid ( const Opm::Deck& deck,
311  const std::vector<double>& poreVolumes = std::vector<double> ())
312  : gridPtr_( createGrid( deck, poreVolumes ) ),
313  grid_( *gridPtr_ ),
314  comm_( *this ),
315  leafIndexSet_( *this ),
316  globalIdSet_( *this ),
317  localIdSet_( *this )
318  {
319  init();
320  }
321 #endif
322 
330  explicit PolyhedralGrid ( const UnstructuredGridType& grid )
331  : gridPtr_(),
332  grid_( grid ),
333  comm_( *this ),
334  leafIndexSet_( *this ),
335  globalIdSet_( *this ),
336  localIdSet_( *this )
337  {
338  init();
339  }
340 
345  operator const UnstructuredGridType& () const { return grid_; }
346 
359  int maxLevel () const
360  {
361  return 1;
362  }
363 
372  int size ( int /* level */, int codim ) const
373  {
374  return size( codim );
375  }
376 
383  int size ( int codim ) const
384  {
385  if( codim == 0 )
386  {
387  return grid_.number_of_cells;
388  }
389  else if ( codim == 1 )
390  {
391  return grid_.number_of_faces;
392  }
393  else if ( codim == dim )
394  {
395  return grid_.number_of_nodes;
396  }
397  else
398  {
399  std::cerr << "Warning: codimension " << codim << " not available in PolyhedralGrid" << std::endl;
400  return 0;
401  }
402  }
403 
412  int size ( int /* level */, GeometryType type ) const
413  {
414  return size( dim - type.dim() );
415  }
416 
421  int size ( GeometryType type ) const
422  {
423  return size( dim - type.dim() );
424  }
425 
432  size_t numBoundarySegments () const
433  {
434  return 0;
435  }
438  template< int codim >
439  typename Codim< codim >::LeafIterator leafbegin () const
440  {
441  return leafbegin< codim, All_Partition >();
442  }
443 
444  template< int codim >
445  typename Codim< codim >::LeafIterator leafend () const
446  {
447  return leafend< codim, All_Partition >();
448  }
449 
450  template< int codim, PartitionIteratorType pitype >
451  typename Codim< codim >::template Partition< pitype >::LeafIterator
452  leafbegin () const
453  {
454  typedef typename Traits::template Codim< codim >::template Partition< pitype >::LeafIteratorImpl Impl;
455  return Impl( extraData(), true );
456  }
457 
458  template< int codim, PartitionIteratorType pitype >
459  typename Codim< codim >::template Partition< pitype >::LeafIterator
460  leafend () const
461  {
462  typedef typename Traits::template Codim< codim >::template Partition< pitype >::LeafIteratorImpl Impl;
463  return Impl( extraData(), false );
464  }
465 
466  template< int codim >
467  typename Codim< codim >::LevelIterator lbegin ( const int /* level */ ) const
468  {
469  return leafbegin< codim, All_Partition >();
470  }
471 
472  template< int codim >
473  typename Codim< codim >::LevelIterator lend ( const int /* level */ ) const
474  {
475  return leafend< codim, All_Partition >();
476  }
477 
478  template< int codim, PartitionIteratorType pitype >
479  typename Codim< codim >::template Partition< pitype >::LevelIterator
480  lbegin ( const int /* level */ ) const
481  {
482  return leafbegin< codim, pitype > ();
483  }
484 
485  template< int codim, PartitionIteratorType pitype >
486  typename Codim< codim >::template Partition< pitype >::LevelIterator
487  lend ( const int /* level */ ) const
488  {
489  return leafend< codim, pitype > ();
490  }
491 
492  const GlobalIdSet &globalIdSet () const
493  {
494  return globalIdSet_;
495  }
496 
497  const LocalIdSet &localIdSet () const
498  {
499  return localIdSet_;
500  }
501 
502  const LevelIndexSet &levelIndexSet ( int /* level */ ) const
503  {
504  return leafIndexSet();
505  }
506 
507  const LeafIndexSet &leafIndexSet () const
508  {
509  return leafIndexSet_;
510  }
511 
512  void globalRefine ( int /* refCount */ )
513  {
514  }
515 
516  bool mark ( int /* refCount */, const typename Codim< 0 >::Entity& /* entity */ )
517  {
518  return false;
519  }
520 
521  int getMark ( const typename Codim< 0 >::Entity& /* entity */) const
522  {
523  return false;
524  }
525 
527  bool preAdapt ()
528  {
529  return false;
530  }
531 
533  bool adapt ()
534  {
535  return false ;
536  }
537 
542  template< class DataHandle >
543  bool adapt ( DataHandle & )
544  {
545  return false;
546  }
547 
549  void postAdapt ()
550  {
551  }
552 
560  int overlapSize ( int /* codim */) const
561  {
562  return 0;
563  }
564 
569  int ghostSize( int codim ) const
570  {
571  return (codim == 0 ) ? 1 : 0;
572  }
573 
579  int overlapSize ( int /* level */, int /* codim */ ) const
580  {
581  return 0;
582  }
583 
589  int ghostSize ( int /* level */, int codim ) const
590  {
591  return ghostSize( codim );
592  }
593 
607  template< class DataHandle, class Data >
608  void communicate ( CommDataHandleIF< DataHandle, Data >& /* dataHandle */,
609  InterfaceType /* interface */,
610  CommunicationDirection /* direction */,
611  int /* level */ ) const
612  {
613  //levelGridView( level ).communicate( dataHandle, interface, direction );
614  }
615 
628  template< class DataHandle, class Data >
629  void communicate ( CommDataHandleIF< DataHandle, Data >& /* dataHandle */,
630  InterfaceType /* interface */,
631  CommunicationDirection /* direction */ ) const
632  {
633  //leafGridView().communicate( dataHandle, interface, direction );
634  }
635 
645  {
646  return comm_;
647  }
648 
649  // data handle interface different between geo and interface
650 
660  bool loadBalance ()
661  {
662  return false ;
663  }
664 
680  template< class DataHandle, class Data >
681  bool loadBalance ( CommDataHandleIF< DataHandle, Data >& /* datahandle */ )
682  {
683  return false;
684  }
685 
700  template< class DofManager >
701  bool loadBalance ( DofManager& /* dofManager */ )
702  {
703  return false;
704  }
705 
707  template< PartitionIteratorType pitype >
708  typename Partition< pitype >::LevelGridView levelGridView ( int /* level */ ) const
709  {
710  typedef typename Partition< pitype >::LevelGridView View;
711  typedef typename View::GridViewImp ViewImp;
712  return View( ViewImp( *this ) );
713  }
714 
716  template< PartitionIteratorType pitype >
717  typename Partition< pitype >::LeafGridView leafGridView () const
718  {
719  typedef typename Traits::template Partition< pitype >::LeafGridView View;
720  typedef typename View::GridViewImp ViewImp;
721  return View( ViewImp( *this ) );
722  }
723 
725  LevelGridView levelGridView ( int /* level */ ) const
726  {
727  typedef typename LevelGridView::GridViewImp ViewImp;
728  return LevelGridView( ViewImp( *this ) );
729  }
730 
732  LeafGridView leafGridView () const
733  {
734  typedef typename LeafGridView::GridViewImp ViewImp;
735  return LeafGridView( ViewImp( *this ) );
736  }
737 
739  template< class EntitySeed >
740  typename Traits::template Codim< EntitySeed::codimension >::EntityPointer
741  entityPointer ( const EntitySeed &seed ) const
742  {
743  typedef typename Traits::template Codim< EntitySeed::codimension >::EntityPointer EntityPointer;
744  typedef typename Traits::template Codim< EntitySeed::codimension >::EntityPointerImpl EntityPointerImpl;
745  typedef typename Traits::template Codim< EntitySeed::codimension >::EntityImpl EntityImpl;
746  return EntityPointer( EntityPointerImpl( EntityImpl( extraData(), seed ) ) );
747  }
748 
750  template< class EntitySeed >
751  typename Traits::template Codim< EntitySeed::codimension >::Entity
752  entity ( const EntitySeed &seed ) const
753  {
754  typedef typename Traits::template Codim< EntitySeed::codimension >::EntityImpl EntityImpl;
755  return EntityImpl( extraData(), seed );
756  }
757 
771  void update ()
772  {
773  }
774 
777  const std::array<int, 3>& logicalCartesianSize() const
778  {
779  return cartDims_;
780  }
781 
782  const int* globalCell() const
783  {
784  assert( grid_.global_cell != 0 );
785  return grid_.global_cell;
786  }
787 
788  void getIJK(const int c, std::array<int,3>& ijk) const
789  {
790  int gc = globalCell()[c];
791  ijk[0] = gc % logicalCartesianSize()[0]; gc /= logicalCartesianSize()[0];
792  ijk[1] = gc % logicalCartesianSize()[1];
793  ijk[2] = gc / logicalCartesianSize()[1];
794  }
795 
796  protected:
797 
798 #if HAVE_OPM_PARSER
799  UnstructuredGridType* createGrid( const Opm::Deck& deck, const std::vector< double >& poreVolumes ) const
800  {
801  const int* rawactnum = deck.hasKeyword("ACTNUM")
802  ? deck.getKeyword("ACTNUM").getIntData().data()
803  : nullptr;
804  const auto eclipseGrid = std::make_shared<Opm::EclipseGrid>(deck, rawactnum);
805 
806  struct grdecl g;
807  std::vector<int> actnum;
808  std::vector<double> coord;
809  std::vector<double> zcorn;
810  std::vector<double> mapaxes;
811 
812  g.dims[0] = eclipseGrid->getNX();
813  g.dims[1] = eclipseGrid->getNY();
814  g.dims[2] = eclipseGrid->getNZ();
815 
816  eclipseGrid->exportMAPAXES( mapaxes );
817  eclipseGrid->exportCOORD( coord );
818  eclipseGrid->exportZCORN( zcorn );
819  eclipseGrid->exportACTNUM( actnum );
820 
821  g.coord = coord.data();
822  g.zcorn = zcorn.data();
823  g.actnum = actnum.data();
824  g.mapaxes = mapaxes.data();
825 
826  if (!poreVolumes.empty() && (eclipseGrid->getMinpvMode() != Opm::MinpvMode::ModeEnum::Inactive)) {
827  Opm::MinpvProcessor mp(g.dims[0], g.dims[1], g.dims[2]);
828  const double minpv_value = eclipseGrid->getMinpvValue();
829  // Currently the pinchProcessor is not used and only opmfil is supported
830  //bool opmfil = eclipseGrid->getMinpvMode() == Opm::MinpvMode::OpmFIL;
831  bool opmfil = true;
832  mp.process(poreVolumes, minpv_value, actnum, opmfil, zcorn.data());
833  }
834 
835  const double z_tolerance = eclipseGrid->isPinchActive() ?
836  eclipseGrid->getPinchThresholdThickness() : 0.0;
837  UnstructuredGridType* cgrid = create_grid_cornerpoint(&g, z_tolerance);
838  if (!cgrid) {
839  OPM_THROW(std::runtime_error, "Failed to construct grid.");
840  }
841  return cgrid;
842  }
843 #endif
844 
845  public:
846  using Base::getRealImplementation;
847 
848  typedef typename Traits :: ExtraData ExtraData;
849  ExtraData extraData () const { return this; }
850 
851  template <class EntitySeed>
852  int corners( const EntitySeed& seed ) const
853  {
854  const int codim = EntitySeed :: codimension;
855  const int index = seed.index();
856  switch (codim)
857  {
858  case 0:
859  {
860  return cellVertices_[ index ].size();
861  }
862  case 1:
863  {
864  return 0;//grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ];
865  }
866  case dim:
867  {
868  return 1;
869  }
870  }
871  return 0;
872  }
873 
874  template <class EntitySeed>
875  GlobalCoordinate
876  corner( const EntitySeed& seed, const int i ) const
877  {
878  const int codim = EntitySeed :: codimension;
879  switch (codim)
880  {
881  case 0:
882  {
883  const int coordIndex = GlobalCoordinate :: dimension * cellVertices_[ seed.index() ][ i ];
884  return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex );
885  }
886  case 1:
887  {
888  const int faceVertex = grid_.face_nodes[grid_.face_nodepos[seed.index()] + i];
889  return copyToGlobalCoordinate( grid_.node_coordinates + GlobalCoordinate :: dimension * faceVertex );
890  }
891  case dim:
892  {
893  const int coordIndex = GlobalCoordinate :: dimension * seed.index();
894  return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex );
895  }
896  }
897  return GlobalCoordinate( 0 );
898  }
899 
900  template <class EntitySeed>
901  int subEntities( const EntitySeed& seed, const int codim ) const
902  {
903  const int index = seed.index();
904  switch (codim)
905  {
906  case 0:
907  return 1;
908  case 1:
909  return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ];
910  case dim:
911  return cellVertices_[ index ].size();
912  }
913  return 0;
914  }
915 
916  template <int codim, class EntitySeedArg >
917  typename Codim<codim>::EntitySeed
918  subEntitySeed( const EntitySeedArg& baseSeed, const int i ) const
919  {
920  assert( codim >= EntitySeedArg::codimension );
921  assert( i>= 0 && i<subEntities( baseSeed, codim ) );
922  typedef typename Codim<codim>::EntitySeed EntitySeed;
923 
924  // if codim equals entity seed codim just return same entity seed.
925  if( codim == EntitySeedArg::codimension )
926  {
927  return EntitySeed( baseSeed.index() );
928  }
929 
930  if( EntitySeedArg::codimension == 0 )
931  {
932  if ( codim == 1 )
933  {
934  return EntitySeed( grid_.cell_faces[ grid_.cell_facepos[ baseSeed.index() ] + i ] );
935  }
936  else if ( codim == dim )
937  {
938  return EntitySeed( cellVertices_[ baseSeed.index() ][ i ] );
939  }
940  }
941  else if ( EntitySeedArg::codimension == 1 && codim == dim )
942  {
943  return EntitySeed( grid_.face_nodes[ grid_.face_nodepos[ baseSeed.index() + i ] ]);
944  }
945 
946  DUNE_THROW(NotImplemented,"codimension not available");
947  return EntitySeed();
948  }
949 
950  const std::vector< GeometryType > &geomTypes ( const unsigned int codim ) const
951  {
952  static std::vector< GeometryType > emptyDummy;
953  if (0 <= codim && codim < geomTypes_.size())
954  {
955  return geomTypes_[codim];
956  }
957 
958  return emptyDummy;
959  }
960 
961  int indexInInside( const typename Codim<0>::EntitySeed& seed, const int i ) const
962  {
963  return ( grid_.cell_facetag ) ? cartesianIndexInInside( seed, i ) : i;
964  }
965 
966  int cartesianIndexInInside( const typename Codim<0>::EntitySeed& seed, const int i ) const
967  {
968  assert( i>= 0 && i<subEntities( seed, 1 ) );
969  return grid_.cell_facetag[ grid_.cell_facepos[ seed.index() ] + i ] ;
970  }
971 
972  typename Codim<0>::EntitySeed
973  neighbor( const typename Codim<0>::EntitySeed& seed, const int i ) const
974  {
975  const int face = this->template subEntitySeed<1>( seed, i ).index();
976  int nb = grid_.face_cells[ 2 * face ];
977  if( nb == seed.index() )
978  {
979  nb = grid_.face_cells[ 2 * face + 1 ];
980  }
981 
982  typedef typename Codim<0>::EntitySeed EntitySeed;
983  return EntitySeed( nb );
984  }
985 
986  int
987  indexInOutside( const typename Codim<0>::EntitySeed& seed, const int i ) const
988  {
989  if( grid_.cell_facetag )
990  {
991  // if cell_facetag is present we assume pseudo Cartesian corner point case
992  const int in_inside = cartesianIndexInInside( seed, i );
993  return in_inside + ((in_inside % 2) ? -1 : 1);
994  }
995  else
996  {
997  typedef typename Codim<0>::EntitySeed EntitySeed;
998  EntitySeed nb = neighbor( seed, i );
999  const int faces = subEntities( seed, 1 );
1000  for( int face = 0; face<faces; ++ face )
1001  {
1002  if( neighbor( nb, face ).equals(seed) )
1003  {
1004  return indexInInside( nb, face );
1005  }
1006  }
1007  DUNE_THROW(InvalidStateException,"inverse intersection not found");
1008  return -1;
1009  }
1010  }
1011 
1012  template <class EntitySeed>
1013  GlobalCoordinate
1014  outerNormal( const EntitySeed& seed, const int i ) const
1015  {
1016  const int face = this->template subEntitySeed<1>( seed, i ).index();
1017  const int normalIdx = face * GlobalCoordinate :: dimension ;
1018  GlobalCoordinate normal = copyToGlobalCoordinate( grid_.face_normals + normalIdx );
1019  const int nb = grid_.face_cells[ 2*face ];
1020  if( nb != seed.index() )
1021  {
1022  normal *= -1.0;
1023  }
1024  return normal;
1025  }
1026 
1027  template <class EntitySeed>
1028  GlobalCoordinate
1029  unitOuterNormal( const EntitySeed& seed, const int i ) const
1030  {
1031  const int face = this->template subEntitySeed<1>( seed, i ).index();
1032  if( seed.index() == grid_.face_cells[ 2*face ] )
1033  {
1034  return unitOuterNormals_[ face ];
1035  }
1036  else
1037  {
1038  GlobalCoordinate normal = unitOuterNormals_[ face ];
1039  normal *= -1.0;
1040  return normal;
1041  }
1042  }
1043 
1044  template <class EntitySeed>
1045  GlobalCoordinate centroids( const EntitySeed& seed ) const
1046  {
1047  const int index = GlobalCoordinate :: dimension * seed.index();
1048  const int codim = EntitySeed::codimension;
1049  assert( index >= 0 && index < size( codim ) * GlobalCoordinate :: dimension );
1050 
1051  if( codim == 0 )
1052  {
1053  return copyToGlobalCoordinate( grid_.cell_centroids + index );
1054  }
1055  else if ( codim == 1 )
1056  {
1057  return copyToGlobalCoordinate( grid_.face_centroids + index );
1058  }
1059  else if( codim == dim )
1060  {
1061  return copyToGlobalCoordinate( grid_.node_coordinates + index );
1062  }
1063  else
1064  {
1065  DUNE_THROW(InvalidStateException,"codimension not implemented");
1066  return GlobalCoordinate( 0 );
1067  }
1068  }
1069 
1070  GlobalCoordinate copyToGlobalCoordinate( const double* coords ) const
1071  {
1072  GlobalCoordinate coordinate;
1073  for( int i=0; i<GlobalCoordinate::dimension; ++i )
1074  {
1075  coordinate[ i ] = coords[ i ];
1076  }
1077  return coordinate;
1078  }
1079 
1080  template <class EntitySeed>
1081  double volumes( const EntitySeed& seed ) const
1082  {
1083  const int index = seed.index();
1084  const int codim = EntitySeed::codimension;
1085  if( codim == 0 )
1086  {
1087  return grid_.cell_volumes[ index ];
1088  }
1089  else if ( codim == 1 )
1090  {
1091  return grid_.face_areas[ index ];
1092  }
1093  else if ( codim == dim )
1094  {
1095  return 1.0;
1096  }
1097  else
1098  {
1099  DUNE_THROW(InvalidStateException,"codimension not implemented");
1100  return 0.0;
1101  }
1102  }
1103 
1104  void init()
1105  {
1106  // copy Cartesian dimensions
1107  for( int i=0; i<3; ++i )
1108  {
1109  cartDims_[ i ] = grid_.cartdims[ i ];
1110  }
1111 
1112  // setup list of cell vertices
1113  const int numCells = size( 0 );
1114  cellVertices_.resize( numCells );
1115 
1116  // sort vertices such that they comply with the dune cube reference element
1117  if( grid_.cell_facetag )
1118  {
1119  typedef Dune::array<int, 3> KeyType;
1120  std::map< const KeyType, const int > vertexFaceTags;
1121  const int vertexFacePattern [8][3] = {
1122  { 0, 2, 4 }, // vertex 0
1123  { 1, 2, 4 }, // vertex 1
1124  { 0, 3, 4 }, // vertex 2
1125  { 1, 3, 4 }, // vertex 3
1126  { 0, 2, 5 }, // vertex 4
1127  { 1, 2, 5 }, // vertex 5
1128  { 0, 3, 5 }, // vertex 6
1129  { 1, 3, 5 } // vertex 7
1130  };
1131 
1132  for( int i=0; i<8; ++i )
1133  {
1134  KeyType key; key.fill( 4 ); // default is 4 which is the first z coord (for the 2d case)
1135  for( int j=0; j<dim; ++j )
1136  {
1137  key[ j ] = vertexFacePattern[ i ][ j ];
1138  }
1139 
1140  vertexFaceTags.insert( std::make_pair( key, i ) );
1141  }
1142 
1143  for (int c = 0; c < numCells; ++c)
1144  {
1145  typedef std::map<int,int> vertexmap_t;
1146  typedef typename vertexmap_t :: iterator iterator;
1147 
1148  std::vector< vertexmap_t > cell_pts( dim*2 );
1149 
1150  for (int hf=grid_.cell_facepos[ c ]; hf < grid_.cell_facepos[c+1]; ++hf)
1151  {
1152  const int f = grid_.cell_faces[ hf ];
1153  const int faceTag = grid_.cell_facetag[ hf ];
1154 
1155  for( int nodepos=grid_.face_nodepos[f]; nodepos<grid_.face_nodepos[f+1]; ++nodepos )
1156  {
1157  const int node = grid_.face_nodes[ nodepos ];
1158  iterator it = cell_pts[ faceTag ].find( node );
1159  if( it == cell_pts[ faceTag ].end() )
1160  {
1161  cell_pts[ faceTag ].insert( std::make_pair( node, 1 ) );
1162  }
1163  else
1164  {
1165  // increase vertex reference counter
1166  (*it).second++;
1167  }
1168  }
1169  }
1170 
1171  typedef std::map< int, std::set<int> > vertexlist_t;
1172  vertexlist_t vertexList;
1173 
1174  for( int faceTag = 0; faceTag<6; ++faceTag )
1175  {
1176  for( iterator it = cell_pts[ faceTag ].begin(),
1177  end = cell_pts[ faceTag ].end(); it != end; ++it )
1178  {
1179  // only consider vertices with one appearance
1180  if( (*it).second == 1 )
1181  {
1182  vertexList[ (*it).first ].insert( faceTag );
1183  }
1184  }
1185  }
1186 
1187  assert( int(vertexList.size()) == ( dim == 2 ) ? 4 : 8 );
1188 
1189  cellVertices_[ c ].resize( vertexList.size() );
1190  for( auto it = vertexList.begin(), end = vertexList.end(); it != end; ++it )
1191  {
1192  assert( (*it).second.size() == dim );
1193  KeyType key; key.fill( 4 ); // fill with 4 which is the first z coord
1194 
1195  std::copy( (*it).second.begin(), (*it).second.end(), key.begin() );
1196  auto vx = vertexFaceTags.find( key );
1197  assert( vx != vertexFaceTags.end() );
1198  if( vx != vertexFaceTags.end() )
1199  {
1200  if( (*vx).second >= int(cellVertices_[ c ].size()) )
1201  cellVertices_[ c ].resize( (*vx).second+1 );
1202  // store node number on correct local position
1203  cellVertices_[ c ][ (*vx).second ] = (*it).first ;
1204  }
1205  }
1206  }
1207  // if face_tag is available we assume that the elements follow a cube-like structure
1208  geomTypes_.resize(dim + 1);
1209  GeometryType tmp;
1210  for (int codim = 0; codim <= dim; ++codim)
1211  {
1212  tmp.makeCube(dim - codim);
1213  geomTypes_[codim].push_back(tmp);
1214  }
1215  }
1216  else // if ( grid_.cell_facetag )
1217  {
1218  for (int c = 0; c < numCells; ++c)
1219  {
1220  std::set<int> cell_pts;
1221  for (int hf=grid_.cell_facepos[ c ]; hf < grid_.cell_facepos[c+1]; ++hf)
1222  {
1223  int f = grid_.cell_faces[ hf ];
1224  const int* fnbeg = grid_.face_nodes + grid_.face_nodepos[f];
1225  const int* fnend = grid_.face_nodes + grid_.face_nodepos[f+1];
1226  cell_pts.insert(fnbeg, fnend);
1227  }
1228 
1229  cellVertices_[ c ].resize( cell_pts.size() );
1230  std::copy(cell_pts.begin(), cell_pts.end(), cellVertices_[ c ].begin() );
1231  }
1232  // if no face_tag is available we assume that no reference element can be
1233  // assigned to the elements
1234  geomTypes_.resize(dim + 1);
1235  GeometryType tmp;
1236  for (int codim = 0; codim <= dim; ++codim)
1237  {
1238  if( codim == dim )
1239  {
1240  tmp.makeCube(dim - codim);
1241  }
1242  else
1243  {
1244  tmp.makeNone(dim - codim);
1245  }
1246  geomTypes_[codim].push_back(tmp);
1247  }
1248  } // end else of ( grid_.cell_facetag )
1249 
1250  unitOuterNormals_.resize( grid_.number_of_faces );
1251  for( int face = 0; face < grid_.number_of_faces; ++face )
1252  {
1253  const int normalIdx = face * GlobalCoordinate :: dimension ;
1254  GlobalCoordinate normal = copyToGlobalCoordinate( grid_.face_normals + normalIdx );
1255  normal /= normal.two_norm();
1256 
1257  unitOuterNormals_[ face ] = normal;
1258  }
1259  }
1260 
1261  protected:
1262  std::unique_ptr< UnstructuredGridType, UnstructuredGridDeleter > gridPtr_;
1263  const UnstructuredGridType& grid_;
1264 
1266  std::array< int, 3 > cartDims_;
1267  std::vector< std::vector< GeometryType > > geomTypes_;
1268  std::vector< std::vector< int > > cellVertices_;
1269 
1270  std::vector< GlobalCoordinate > unitOuterNormals_;
1271 
1272  mutable LeafIndexSet leafIndexSet_;
1273  mutable GlobalIdSet globalIdSet_;
1274  mutable LocalIdSet localIdSet_;
1275 
1276  private:
1277  // no copying
1278  PolyhedralGrid ( const PolyhedralGrid& );
1279  };
1280 
1281 
1282 
1283  // PolyhedralGrid::Codim
1284  // -------------
1285 
1286  template< int dim, int dimworld >
1287  template< int codim >
1288  struct PolyhedralGrid< dim, dimworld >::Codim
1289  : public Base::template Codim< codim >
1290  {
1299  typedef typename Traits::template Codim< codim >::Entity Entity;
1300 
1305  typedef typename Traits::template Codim< codim >::EntityPointer EntityPointer;
1306 
1320  typedef typename Traits::template Codim< codim >::Geometry Geometry;
1321 
1330  typedef typename Traits::template Codim< codim >::LocalGeometry LocalGeometry;
1331 
1337  template< PartitionIteratorType pitype >
1338  struct Partition
1339  {
1340  typedef typename Traits::template Codim< codim >
1341  ::template Partition< pitype >::LeafIterator
1342  LeafIterator;
1343  typedef typename Traits::template Codim< codim >
1344  ::template Partition< pitype >::LevelIterator
1345  LevelIterator;
1346  };
1347 
1355  typedef typename Partition< All_Partition >::LeafIterator LeafIterator;
1356 
1364  typedef typename Partition< All_Partition >::LevelIterator LevelIterator;
1365 
1367  };
1368 
1369 } // namespace Dune
1370 
1371 #include <dune/grid/polyhedralgrid/persistentcontainer.hh>
1372 #include <dune/grid/polyhedralgrid/cartesianindexmapper.hh>
1373 #include <dune/grid/polyhedralgrid/gridhelpers.hh>
1374 
1375 #endif // #ifndef DUNE_POLYHEDRALGRID_GRID_HH
int * cell_facetag
If non-null, this array contains a number for cell-face adjacency indicating the face&#39;s position with...
Definition: grid.h:244
void communicate(CommDataHandleIF< DataHandle, Data > &, InterfaceType, CommunicationDirection, int) const
communicate information on a grid level
Definition: grid.hh:608
Definition: idset.hh:15
int * global_cell
If non-null, this array contains the logical cartesian indices (in a lexicographic ordering) of each ...
Definition: grid.h:214
int size(int codim) const
obtain number of leaf entities
Definition: grid.hh:383
Definition: intersection.hh:19
Definition: geometry.hh:24
Raw corner-point specification of a particular geological model.
Definition: preprocess.h:56
int size(int, GeometryType type) const
obtain number of entites on a level
Definition: grid.hh:412
size_t numBoundarySegments() const
obtain number of leaf entities
Definition: grid.hh:432
Traits::template Codim< codim >::EntityPointer EntityPointer
type of entity pointer
Definition: grid.hh:1305
const int * actnum
Explicit &quot;active&quot; map.
Definition: preprocess.h:60
double * cell_volumes
Exact or approximate cell volumes.
Definition: grid.h:197
Definition: entitypointer.hh:18
double * node_coordinates
Node coordinates, stored consecutively for each node.
Definition: grid.h:160
Traits::LevelIntersectionIterator LevelIntersectionIterator
iterator over intersections with other entities on the same level
Definition: grid.hh:213
double * face_centroids
Exact or approximate face centroids, stored consecutively for each face.
Definition: grid.h:168
Partition< pitype >::LeafGridView leafGridView() const
View for the leaf grid.
Definition: grid.hh:717
int number_of_cells
The number of cells in the grid.
Definition: grid.h:109
int * cell_facepos
For a cell c, cell_facepos[c] contains the starting index for c&#39;s faces in the cell_faces array...
Definition: grid.h:152
Partition< pitype >::LevelGridView levelGridView(int) const
View for a grid level.
Definition: grid.hh:708
int overlapSize(int) const
obtain size of overlap region for the leaf grid
Definition: grid.hh:560
PolyhedralGrid(const UnstructuredGridType &grid)
constructor
Definition: grid.hh:330
Traits::HierarchicIterator HierarchicIterator
iterator over the grid hierarchy
Definition: grid.hh:209
bool preAdapt()
Definition: grid.hh:527
Partition< All_Partition >::LevelIterator LevelIterator
type of leaf iterator
Definition: grid.hh:1364
Definition: iterator.hh:19
double * cell_centroids
Exact or approximate cell centroids, stored consecutively for each cell.
Definition: grid.h:192
const CollectiveCommunication & comm() const
obtain CollectiveCommunication object
Definition: grid.hh:644
Traits::LocalIdSet LocalIdSet
type of local id set
Definition: grid.hh:284
double * face_normals
Exact or approximate face normals, stored consecutively for each face.
Definition: grid.h:184
LeafGridView leafGridView() const
View for the leaf grid for All_Partition.
Definition: grid.hh:732
int size(int, int codim) const
obtain number of entites on a level
Definition: grid.hh:372
const double * coord
Pillar end-points.
Definition: preprocess.h:58
const double * mapaxes
6 Element rotation vector - can be NULL.
Definition: preprocess.h:61
Traits::LeafIndexSet LeafIndexSet
type of leaf index set
Definition: grid.hh:245
Traits::LeafIntersectionIterator LeafIntersectionIterator
iterator over intersections with other entities on the leaf level
Definition: grid.hh:211
int number_of_faces
The number of faces in the grid.
Definition: grid.h:111
bool loadBalance(DofManager &)
rebalance the load each process has to handle
Definition: grid.hh:701
Types for GridView.
Definition: grid.hh:222
traits structure containing types for a codimension
Definition: grid.hh:201
bool adapt()
Definition: grid.hh:533
int overlapSize(int, int) const
obtain size of overlap region for a grid level
Definition: grid.hh:579
Traits::template Codim< EntitySeed::codimension >::EntityPointer entityPointer(const EntitySeed &seed) const
obtain EntityPointer from EntitySeed.
Definition: grid.hh:741
bool adapt(DataHandle &)
Definition: grid.hh:543
Traits::template Codim< codim >::LocalGeometry LocalGeometry
type of local geometry
Definition: grid.hh:1330
bool loadBalance(CommDataHandleIF< DataHandle, Data > &)
rebalance the load each process has to handle
Definition: grid.hh:681
int maxLevel() const
obtain maximal grid level
Definition: grid.hh:359
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
Definition: entityseed.hh:17
int size(GeometryType type) const
returns the number of boundary segments within the macro grid
Definition: grid.hh:421
Definition: intersectioniterator.hh:15
identical grid wrapper
Definition: declaration.hh:10
int * face_nodepos
For a face f, face_nodepos[f] contains the starting index for f&#39;s nodes in the face_nodes array...
Definition: grid.h:127
Traits::ctype ctype
type of vector coordinates (e.g., double)
Definition: grid.hh:292
Traits::CollectiveCommunication CollectiveCommunication
communicator with all other processes having some part of the grid
Definition: grid.hh:295
const double * zcorn
Corner-point depths.
Definition: preprocess.h:59
int ghostSize(int, int codim) const
obtain size of ghost region for a grid level
Definition: grid.hh:589
Main OPM-Core grid data structure along with helper functions for construction, destruction and readi...
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: grid.h:138
void update()
update grid caches
Definition: grid.hh:771
int * cell_faces
Contains for each cell, the indices of its adjacent faces.
Definition: grid.h:146
Definition: indexset.hh:22
int ghostSize(int codim) const
obtain size of ghost region for the leaf grid
Definition: grid.hh:569
LevelGridView levelGridView(int) const
View for a grid level for All_Partition.
Definition: grid.hh:725
Traits::GlobalIdSet GlobalIdSet
type of global id set
Definition: grid.hh:267
double * face_areas
Exact or approximate face areas.
Definition: grid.h:173
Traits::template Codim< codim >::Entity Entity
type of entity
Definition: grid.hh:1299
void postAdapt()
Definition: grid.hh:549
GridFamily::Traits Traits
type of the grid traits
Definition: grid.hh:192
Partition< All_Partition >::LevelGridView LevelGridView
View types for All_Partition.
Definition: grid.hh:229
Definition: entity.hh:151
void communicate(CommDataHandleIF< DataHandle, Data > &, InterfaceType, CommunicationDirection) const
communicate information on leaf entities
Definition: grid.hh:629
Low-level corner-point processing routines and supporting data structures.
Definition: geometry.hh:25
bool loadBalance()
rebalance the load each process has to handle
Definition: grid.hh:660
Partition< All_Partition >::LeafIterator LeafIterator
type of level iterator
Definition: grid.hh:1355
int * face_nodes
Contains for each face, the indices of its adjacent nodes.
Definition: grid.h:121
Routines to form a complete UnstructuredGrid from a corner-point specification.
Data structure for an unstructured grid, unstructured meaning that any cell may have an arbitrary num...
Definition: grid.h:98
Traits::template Codim< EntitySeed::codimension >::Entity entity(const EntitySeed &seed) const
obtain EntityPointer from EntitySeed.
Definition: grid.hh:752
Traits::LevelIndexSet LevelIndexSet
type of level index set
Definition: grid.hh:255
int number_of_nodes
The number of nodes in the grid.
Definition: grid.h:113
Transform a corner-point grid ZCORN field to account for MINPV processing.
Definition: MinpvProcessor.hpp:32
Definition: grid.hh:51
void destroy_grid(struct UnstructuredGrid *g)
Destroy and deallocate an UnstructuredGrid and all its data.
Definition: grid.c:32
Traits::template Codim< codim >::Geometry Geometry
type of world geometry
Definition: grid.hh:1320
int cartdims[3]
Contains the size of the logical cartesian structure (if any) of the grid.
Definition: grid.h:227