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 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 7)
126  using Communication = Dune::Communication<MPICommunicator>;
127  using CollectiveCommunication = Dune::Communication<MPICommunicator>;
128 #else
129  using CollectiveCommunication = Dune::CollectiveCommunication<MPICommunicator>;
130 #endif
131 
132  template< PartitionIteratorType pitype >
133  struct Partition
134  {
135  typedef Dune::GridView< PolyhedralGridViewTraits< dim, dimworld, ctype, pitype > > LeafGridView;
136  typedef Dune::GridView< PolyhedralGridViewTraits< dim, dimworld, ctype, pitype > > LevelGridView;
137  };
138 
139  typedef typename Partition< All_Partition >::LeafGridView LeafGridView;
140  typedef typename Partition< All_Partition >::LevelGridView LevelGridView;
141  };
142  };
143 
144 
145 
146  // PolyhedralGrid
147  // --------------
148 
157  template < int dim, int dimworld, typename coord_t >
160  : public GridDefaultImplementation
161  < dim, dimworld, coord_t, PolyhedralGridFamily< dim, dimworld, coord_t > >
163  {
165 
166  typedef GridDefaultImplementation
167  < dim, dimworld, coord_t, PolyhedralGridFamily< dim, dimworld, coord_t > > Base;
168 
169  public:
171 
172  protected:
174  {
175  inline void operator () ( UnstructuredGridType* grdPtr )
176  {
177  destroy_grid( grdPtr );
178  }
179  };
180 
181  public:
182  typedef std::unique_ptr< UnstructuredGridType, UnstructuredGridDeleter > UnstructuredGridPtr;
183 
184  static UnstructuredGridPtr
185  allocateGrid ( std::size_t nCells, std::size_t nFaces, std::size_t nFaceNodes, std::size_t nCellFaces, std::size_t nNodes )
186  {
187  // Note that we here assign a grid of dimension dimworld in order to obtain global coordinates in the correct dimension
188  UnstructuredGridType *grid = allocate_grid( dimworld, nCells, nFaces, nFaceNodes, nCellFaces, nNodes );
189  if( !grid )
190  DUNE_THROW( GridError, "Unable to allocate grid" );
191  return UnstructuredGridPtr( grid );
192  }
193 
194  static void
195  computeGeometry ( UnstructuredGridPtr& ug )
196  {
197  // get C pointer to UnstructuredGrid
198  UnstructuredGrid* ugPtr = ug.operator ->();
199 
200  // compute geometric quantities like cell volume and face normals
201  compute_geometry( ugPtr );
202  }
203 
205  typedef PolyhedralGridFamily< dim, dimworld, coord_t > GridFamily;
212  typedef typename GridFamily::Traits Traits;
213 
220  template< int codim >
221  struct Codim;
222 
229  typedef typename Traits::HierarchicIterator HierarchicIterator;
231  typedef typename Traits::LeafIntersectionIterator LeafIntersectionIterator;
233  typedef typename Traits::LevelIntersectionIterator LevelIntersectionIterator;
234 
241  template< PartitionIteratorType pitype >
242  struct Partition
243  {
244  typedef typename GridFamily::Traits::template Partition< pitype >::LevelGridView LevelGridView;
245  typedef typename GridFamily::Traits::template Partition< pitype >::LeafGridView LeafGridView;
246  };
247 
249  typedef typename Partition< All_Partition >::LevelGridView LevelGridView;
250  typedef typename Partition< All_Partition >::LeafGridView LeafGridView;
251 
265  typedef typename Traits::LeafIndexSet LeafIndexSet;
266 
275  typedef typename Traits::LevelIndexSet LevelIndexSet;
276 
287  typedef typename Traits::GlobalIdSet GlobalIdSet;
288 
304  typedef typename Traits::LocalIdSet LocalIdSet;
305 
312  typedef typename Traits::ctype ctype;
313 
315 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 7)
316  using Communication = typename Traits::Communication;
317  using CommunicationType = Communication;
318 #else
319  using CollectiveCommunication = typename Traits::CollectiveCommunication;
320  using Communication = CollectiveCommunication;
321  using CommunicationType = CollectiveCommunication;
322 #endif
323 
324  typedef typename Traits :: GlobalCoordinate GlobalCoordinate;
325 
331 #if HAVE_ECL_INPUT
337  explicit PolyhedralGrid ( const Opm::EclipseGrid& inputGrid,
338  const std::vector<double>& poreVolumes = std::vector<double> ())
339  : gridPtr_( createGrid( inputGrid, poreVolumes ) ),
340  grid_( *gridPtr_ ),
341  comm_( MPIHelper::getCommunicator() ),
342  leafIndexSet_( *this ),
343  globalIdSet_( *this ),
344  localIdSet_( *this ),
345  nBndSegments_( 0 )
346  {
347  init();
348  }
349 #endif
350 
356  explicit PolyhedralGrid ( const std::vector< int >& n,
357  const std::vector< double >& dx )
358  : gridPtr_( createGrid( n, dx ) ),
359  grid_( *gridPtr_ ),
360  comm_( MPIHelper::getCommunicator()),
361  leafIndexSet_( *this ),
362  globalIdSet_( *this ),
363  localIdSet_( *this ),
364  nBndSegments_( 0 )
365  {
366  init();
367  }
368 
375  explicit PolyhedralGrid ( UnstructuredGridPtr &&gridPtr )
376  : gridPtr_( std::move( gridPtr ) ),
377  grid_( *gridPtr_ ),
378  comm_( MPIHelper::getCommunicator() ),
379  leafIndexSet_( *this ),
380  globalIdSet_( *this ),
381  localIdSet_( *this ),
382  nBndSegments_( 0 )
383  {
384  init();
385  }
386 
394  explicit PolyhedralGrid ( const UnstructuredGridType& grid )
395  : gridPtr_(),
396  grid_( grid ),
397  comm_( MPIHelper::getCommunicator() ),
398  leafIndexSet_( *this ),
399  globalIdSet_( *this ),
400  localIdSet_( *this ),
401  nBndSegments_( 0 )
402  {
403  init();
404  }
405 
410  operator const UnstructuredGridType& () const { return grid_; }
411 
424  int maxLevel () const
425  {
426  return 1;
427  }
428 
437  int size ( int /* level */, int codim ) const
438  {
439  return size( codim );
440  }
441 
448  int size ( int codim ) const
449  {
450  if( codim == 0 )
451  {
452  return grid_.number_of_cells;
453  }
454  else if ( codim == 1 )
455  {
456  return grid_.number_of_faces;
457  }
458  else if ( codim == dim )
459  {
460  return grid_.number_of_nodes;
461  }
462  else
463  {
464  std::cerr << "Warning: codimension " << codim << " not available in PolyhedralGrid" << std::endl;
465  return 0;
466  }
467  }
468 
477  int size ( int /* level */, GeometryType type ) const
478  {
479  return size( dim - type.dim() );
480  }
481 
486  int size ( GeometryType type ) const
487  {
488  return size( dim - type.dim() );
489  }
490 
497  size_t numBoundarySegments () const
498  {
499  return nBndSegments_;
500  }
503  template< int codim >
504  typename Codim< codim >::LeafIterator leafbegin () const
505  {
506  return leafbegin< codim, All_Partition >();
507  }
508 
509  template< int codim >
510  typename Codim< codim >::LeafIterator leafend () const
511  {
512  return leafend< codim, All_Partition >();
513  }
514 
515  template< int codim, PartitionIteratorType pitype >
516  typename Codim< codim >::template Partition< pitype >::LeafIterator
517  leafbegin () const
518  {
519  typedef typename Traits::template Codim< codim >::template Partition< pitype >::LeafIteratorImpl Impl;
520  return Impl( extraData(), true );
521  }
522 
523  template< int codim, PartitionIteratorType pitype >
524  typename Codim< codim >::template Partition< pitype >::LeafIterator
525  leafend () const
526  {
527  typedef typename Traits::template Codim< codim >::template Partition< pitype >::LeafIteratorImpl Impl;
528  return Impl( extraData(), false );
529  }
530 
531  template< int codim >
532  typename Codim< codim >::LevelIterator lbegin ( const int /* level */ ) const
533  {
534  return leafbegin< codim, All_Partition >();
535  }
536 
537  template< int codim >
538  typename Codim< codim >::LevelIterator lend ( const int /* level */ ) const
539  {
540  return leafend< codim, All_Partition >();
541  }
542 
543  template< int codim, PartitionIteratorType pitype >
544  typename Codim< codim >::template Partition< pitype >::LevelIterator
545  lbegin ( const int /* level */ ) const
546  {
547  return leafbegin< codim, pitype > ();
548  }
549 
550  template< int codim, PartitionIteratorType pitype >
551  typename Codim< codim >::template Partition< pitype >::LevelIterator
552  lend ( const int /* level */ ) const
553  {
554  return leafend< codim, pitype > ();
555  }
556 
557  const GlobalIdSet &globalIdSet () const
558  {
559  return globalIdSet_;
560  }
561 
562  const LocalIdSet &localIdSet () const
563  {
564  return localIdSet_;
565  }
566 
567  const LevelIndexSet &levelIndexSet ( int /* level */ ) const
568  {
569  return leafIndexSet();
570  }
571 
572  const LeafIndexSet &leafIndexSet () const
573  {
574  return leafIndexSet_;
575  }
576 
577  void globalRefine ( int /* refCount */ )
578  {
579  }
580 
581  bool mark ( int /* refCount */, const typename Codim< 0 >::Entity& /* entity */ )
582  {
583  return false;
584  }
585 
586  int getMark ( const typename Codim< 0 >::Entity& /* entity */) const
587  {
588  return false;
589  }
590 
592  bool preAdapt ()
593  {
594  return false;
595  }
596 
598  bool adapt ()
599  {
600  return false ;
601  }
602 
607  template< class DataHandle >
608  bool adapt ( DataHandle & )
609  {
610  return false;
611  }
612 
614  void postAdapt ()
615  {
616  }
617 
625  int overlapSize ( int /* codim */) const
626  {
627  return 0;
628  }
629 
634  int ghostSize( int codim ) const
635  {
636  return (codim == 0 ) ? 1 : 0;
637  }
638 
644  int overlapSize ( int /* level */, int /* codim */ ) const
645  {
646  return 0;
647  }
648 
654  int ghostSize ( int /* level */, int codim ) const
655  {
656  return ghostSize( codim );
657  }
658 
672  template< class DataHandle>
673  void communicate ( DataHandle& /* dataHandle */,
674  InterfaceType /* interface */,
675  CommunicationDirection /* direction */,
676  int /* level */ ) const
677  {
678  OPM_THROW(std::runtime_error, "communicate not implemented for polyhedreal grid!");
679  }
680 
693  template< class DataHandle>
694  void communicate ( DataHandle& /* dataHandle */,
695  InterfaceType /* interface */,
696  CommunicationDirection /* direction */ ) const
697  {
698  OPM_THROW(std::runtime_error, "communicate not implemented for polyhedreal grid!");
699  }
700 
703  {
704  OPM_THROW(std::runtime_error, "switch to global view not implemented for polyhedreal grid!");
705  }
706 
709  {
710  OPM_THROW(std::runtime_error, "switch to distributed view not implemented for polyhedreal grid!");
711  }
712 
721  const CommunicationType &comm () const
722  {
723  return comm_;
724  }
725 
726  // data handle interface different between geo and interface
727 
737  bool loadBalance ()
738  {
739  return false ;
740  }
741 
757  template< class DataHandle, class Data >
758  bool loadBalance ( CommDataHandleIF< DataHandle, Data >& /* datahandle */ )
759  {
760  return false;
761  }
762 
777  template< class DofManager >
778  bool loadBalance ( DofManager& /* dofManager */ )
779  {
780  return false;
781  }
782 
784  template< PartitionIteratorType pitype >
785  typename Partition< pitype >::LevelGridView levelGridView ( int /* level */ ) const
786  {
787  typedef typename Partition< pitype >::LevelGridView View;
788  typedef typename View::GridViewImp ViewImp;
789  return View( ViewImp( *this ) );
790  }
791 
793  template< PartitionIteratorType pitype >
794  typename Partition< pitype >::LeafGridView leafGridView () const
795  {
796  typedef typename Traits::template Partition< pitype >::LeafGridView View;
797  typedef typename View::GridViewImp ViewImp;
798  return View( ViewImp( *this ) );
799  }
800 
802  LevelGridView levelGridView ( int /* level */ ) const
803  {
804  typedef typename LevelGridView::GridViewImp ViewImp;
805  return LevelGridView( ViewImp( *this ) );
806  }
807 
809  LeafGridView leafGridView () const
810  {
811  typedef typename LeafGridView::GridViewImp ViewImp;
812  return LeafGridView( ViewImp( *this ) );
813  }
814 
816  template< class EntitySeed >
817  typename Traits::template Codim< EntitySeed::codimension >::EntityPointer
818  entityPointer ( const EntitySeed &seed ) const
819  {
820  typedef typename Traits::template Codim< EntitySeed::codimension >::EntityPointer EntityPointer;
821  typedef typename Traits::template Codim< EntitySeed::codimension >::EntityPointerImpl EntityPointerImpl;
822  typedef typename Traits::template Codim< EntitySeed::codimension >::EntityImpl EntityImpl;
823  return EntityPointer( EntityPointerImpl( EntityImpl( extraData(), seed ) ) );
824  }
825 
827  template< class EntitySeed >
828  typename Traits::template Codim< EntitySeed::codimension >::Entity
829  entity ( const EntitySeed &seed ) const
830  {
831  typedef typename Traits::template Codim< EntitySeed::codimension >::EntityImpl EntityImpl;
832  return EntityImpl( extraData(), seed );
833  }
834 
848  void update ()
849  {
850  }
851 
854  const std::array<int, 3>& logicalCartesianSize() const
855  {
856  return cartDims_;
857  }
858 
859  const int* globalCell() const
860  {
861  assert( grid_.global_cell != 0 );
862  return grid_.global_cell;
863  }
864 
865  const int* globalCellPtr() const
866  {
867  return grid_.global_cell;
868  }
869 
870  void getIJK(const int c, std::array<int,3>& ijk) const
871  {
872  int gc = globalCell()[c];
873  ijk[0] = gc % logicalCartesianSize()[0]; gc /= logicalCartesianSize()[0];
874  ijk[1] = gc % logicalCartesianSize()[1];
875  ijk[2] = gc / logicalCartesianSize()[1];
876  }
877 
882 
883  template<class DataHandle>
893  void scatterData([[maybe_unused]] DataHandle& handle) const
894  {
895  OPM_THROW(std::runtime_error, "ScatterData not implemented for polyhedral grid!");
896  }
897 
898  protected:
899 #if HAVE_ECL_INPUT
900  UnstructuredGridType* createGrid( const Opm::EclipseGrid& inputGrid, const std::vector< double >& poreVolumes ) const
901  {
902  struct grdecl g;
903 
904  g.dims[0] = inputGrid.getNX();
905  g.dims[1] = inputGrid.getNY();
906  g.dims[2] = inputGrid.getNZ();
907 
908  std::vector<double> coord = inputGrid.getCOORD( );
909  std::vector<double> zcorn = inputGrid.getZCORN( );
910  std::vector<int> actnum = inputGrid.getACTNUM( );
911 
912  g.coord = coord.data();
913  g.zcorn = zcorn.data();
914  g.actnum = actnum.data();
915 
916  if (!poreVolumes.empty() && (inputGrid.getMinpvMode() != Opm::MinpvMode::ModeEnum::Inactive))
917  {
918  Opm::MinpvProcessor mp(g.dims[0], g.dims[1], g.dims[2]);
919  const std::vector<double>& minpvv = inputGrid.getMinpvVector();
920  // Currently the pinchProcessor is not used and only opmfil is supported
921  // The polyhedralgrid only only supports the opmfil option
922  //bool opmfil = inputGrid.getMinpvMode() == Opm::MinpvMode::OpmFIL;
923  bool opmfil = true;
924  const size_t cartGridSize = g.dims[0] * g.dims[1] * g.dims[2];
925  std::vector<double> thickness(cartGridSize);
926  for (size_t i = 0; i < cartGridSize; ++i) {
927  thickness[i] = inputGrid.getCellThickness(i);
928  }
929  const double z_tolerance = inputGrid.isPinchActive() ? inputGrid.getPinchThresholdThickness() : 0.0;
930  mp.process(thickness, z_tolerance, poreVolumes, minpvv, actnum, opmfil, zcorn.data());
931  }
932 
933  /*
934  if (!poreVolumes.empty() && (inputGrid.getMinpvMode() != Opm::MinpvMode::ModeEnum::Inactive)) {
935  Opm::MinpvProcessor mp(g.dims[0], g.dims[1], g.dims[2]);
936  const double minpv_value = inputGrid.getMinpvValue();
937  // Currently the pinchProcessor is not used and only opmfil is supported
938  //bool opmfil = inputGrid.getMinpvMode() == Opm::MinpvMode::OpmFIL;
939  bool opmfil = true;
940  mp.process(poreVolumes, minpv_value, actnum, opmfil, zcorn.data());
941  }
942  */
943 
944  const double z_tolerance = inputGrid.isPinchActive() ?
945  inputGrid.getPinchThresholdThickness() : 0.0;
946  UnstructuredGridType* cgrid = create_grid_cornerpoint(&g, z_tolerance);
947  if (!cgrid) {
948  OPM_THROW(std::runtime_error, "Failed to construct grid.");
949  }
950  return cgrid;
951  }
952 #endif
953 
954  UnstructuredGridType* createGrid( const std::vector< int >& n, const std::vector< double >& dx ) const
955  {
956  UnstructuredGridType* cgrid = nullptr ;
957  assert( int(n.size()) == dim );
958  if( dim == 2 )
959  {
960  cgrid = create_grid_cart2d( n[ 0 ], n[ 1 ], dx[ 0 ], dx[ 1 ] );
961  }
962  else if ( dim == 3 )
963  {
964  cgrid = create_grid_hexa3d( n[ 0 ], n[ 1 ], n[ 2 ], dx[ 0 ], dx[ 1 ], dx[ 2 ] );
965  }
966 
967  //print_grid( cgrid );
968  if (!cgrid) {
969  OPM_THROW(std::runtime_error, "Failed to construct grid.");
970  }
971  return cgrid;
972  }
973 
974  public:
975 #if DUNE_VERSION_LT_REV(DUNE_GRID, 2, 7, 1)
976  using Base::getRealImplementation;
977 #endif
978  typedef typename Traits :: ExtraData ExtraData;
979  ExtraData extraData () const { return this; }
980 
981  template <class EntitySeed>
982  int corners( const EntitySeed& seed ) const
983  {
984  const int codim = EntitySeed :: codimension;
985  const int index = seed.index();
986  if (codim==0)
987  return cellVertices_[ index ].size();
988  if (codim==1)
989  return grid_.face_nodepos[ index+1 ] - grid_.face_nodepos[ index ];
990  if (codim==dim)
991  return 1;
992  return 0;
993  }
994 
995  template <class EntitySeed>
996  GlobalCoordinate
997  corner ( const EntitySeed& seed, const int i ) const
998  {
999  const int codim = EntitySeed :: codimension;
1000  if (codim==0)
1001  {
1002  const int coordIndex = GlobalCoordinate :: dimension * cellVertices_[ seed.index() ][ i ];
1003  return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex );
1004  }
1005  if (codim==1)
1006  {
1007  // for faces we need to swap vertices in 3d since in UnstructuredGrid
1008  // those are ordered counter clockwise, for 2d this does not matter
1009  // TODO: Improve this for performance reasons
1010  const int crners = corners( seed );
1011  const int crner = (crners == 4 && EntitySeed :: dimension == 3 && i > 1 ) ? 5 - i : i;
1012  const int faceVertex = grid_.face_nodes[ grid_.face_nodepos[seed.index() ] + crner ];
1013  return copyToGlobalCoordinate( grid_.node_coordinates + GlobalCoordinate :: dimension * faceVertex );
1014  }
1015  if (codim==dim)
1016  {
1017  const int coordIndex = GlobalCoordinate :: dimension * seed.index();
1018  return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex );
1019  }
1020  return GlobalCoordinate( 0 );
1021  }
1022 
1023  template <class EntitySeed>
1024  int subEntities( const EntitySeed& seed, const int codim ) const
1025  {
1026  const int index = seed.index();
1027  if( seed.codimension == 0 )
1028  {
1029  if (codim==0)
1030  return 1;
1031  if (codim==1)
1032  return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ];
1033  if (codim==dim)
1034  return cellVertices_[ index ].size();
1035  }
1036  else if( seed.codimension == 1 )
1037  {
1038  if (codim==1)
1039  return 1;
1040  if (codim==dim)
1041  return grid_.face_nodepos[ index+1 ] - grid_.face_nodepos[ index ];
1042  }
1043  else if ( seed.codimension == dim )
1044  {
1045  return 1 ;
1046  }
1047 
1048  return 0;
1049  }
1050 
1051  template <int codim, class EntitySeedArg >
1052  typename Codim<codim>::EntitySeed
1053  subEntitySeed( const EntitySeedArg& baseSeed, const int i ) const
1054  {
1055  assert( codim >= EntitySeedArg::codimension );
1056  assert( i>= 0 && i<subEntities( baseSeed, codim ) );
1057  typedef typename Codim<codim>::EntitySeed EntitySeed;
1058 
1059  // if codim equals entity seed codim just return same entity seed.
1060  if( codim == EntitySeedArg::codimension )
1061  {
1062  return EntitySeed( baseSeed.index() );
1063  }
1064 
1065  if( EntitySeedArg::codimension == 0 )
1066  {
1067  if ( codim == 1 )
1068  {
1069  return EntitySeed( grid_.cell_faces[ grid_.cell_facepos[ baseSeed.index() ] + i ] );
1070  }
1071  else if ( codim == dim )
1072  {
1073  return EntitySeed( cellVertices_[ baseSeed.index() ][ i ] );
1074  }
1075  }
1076  else if ( EntitySeedArg::codimension == 1 && codim == dim )
1077  {
1078  return EntitySeed( grid_.face_nodes[ grid_.face_nodepos[ baseSeed.index() + i ] ]);
1079  }
1080 
1081  DUNE_THROW(NotImplemented,"codimension not available");
1082  return EntitySeed();
1083  }
1084 
1085  template <int codim>
1086  typename Codim<codim>::EntitySeed
1087  subEntitySeed( const typename Codim<1>::EntitySeed& faceSeed, const int i ) const
1088  {
1089  assert( i>= 0 && i<subEntities( faceSeed, codim ) );
1090  typedef typename Codim<codim>::EntitySeed EntitySeed;
1091  if ( codim == 1 )
1092  {
1093  return EntitySeed( faceSeed.index() );
1094  }
1095  else if ( codim == dim )
1096  {
1097  return EntitySeed( grid_.face_nodes[ grid_.face_nodepos[ faceSeed.index() ] + i ] );
1098  }
1099  else
1100  {
1101  DUNE_THROW(NotImplemented,"codimension not available");
1102  }
1103  }
1104 
1105  bool hasBoundaryIntersections(const typename Codim<0>::EntitySeed& seed ) const
1106  {
1107  const int faces = subEntities( seed, 1 );
1108  for( int f=0; f<faces; ++f )
1109  {
1110  const auto faceSeed = this->template subEntitySeed<1>( seed, f );
1111  if( isBoundaryFace( faceSeed ) )
1112  return true;
1113  }
1114  return false;
1115  }
1116 
1117  bool isBoundaryFace(const int face ) const
1118  {
1119  assert( face >= 0 && face < grid_.number_of_faces );
1120  const int facePos = 2 * face;
1121  return ((grid_.face_cells[ facePos ] < 0) || (grid_.face_cells[ facePos+1 ] < 0));
1122  }
1123 
1124  bool isBoundaryFace(const typename Codim<1>::EntitySeed& faceSeed ) const
1125  {
1126  assert( faceSeed.isValid() );
1127  return isBoundaryFace( faceSeed.index() );
1128  }
1129 
1130  int boundarySegmentIndex(const typename Codim<0>::EntitySeed& seed, const int face ) const
1131  {
1132  const auto faceSeed = this->template subEntitySeed<1>( seed, face );
1133  assert( faceSeed.isValid() );
1134  const int facePos = 2 * faceSeed.index();
1135  const int idx = std::min( grid_.face_cells[ facePos ], grid_.face_cells[ facePos+1 ]);
1136  // check that this is actually the boundary
1137  assert( idx < 0 );
1138  return -(idx+1); // +1 to include 0 boundary segment index
1139  }
1140 
1141  const std::vector< GeometryType > &geomTypes ( const unsigned int codim ) const
1142  {
1143  static std::vector< GeometryType > emptyDummy;
1144  if (codim < geomTypes_.size())
1145  {
1146  return geomTypes_[codim];
1147  }
1148 
1149  return emptyDummy;
1150  }
1151 
1152  template < class Seed >
1153  GeometryType geometryType( const Seed& seed ) const
1154  {
1155  if( Seed::codimension == 0 )
1156  {
1157  assert(!geomTypes( Seed::codimension ).empty());
1158  return geomTypes( Seed::codimension )[ 0 ];
1159  }
1160  else
1161  {
1162  // 3d faces
1163  if( dim == 3 && Seed::codimension == 1 )
1164  {
1165  GeometryType face;
1166  const int nVx = corners( seed );
1167  if( nVx == 4 ) // quad face
1168  face = Dune::GeometryTypes::cube(2);
1169  else if( nVx == 3 ) // triangle face
1170  face = Dune::GeometryTypes::simplex(2);
1171  else // polygonal face
1172  face = Dune::GeometryTypes::none(2);
1173  return face;
1174  }
1175 
1176  // for codim 2 and codim 3 there is only one option
1177  assert(!geomTypes( Seed::codimension ).empty());
1178  return geomTypes( Seed::codimension )[ 0 ];
1179  }
1180  }
1181 
1182  int indexInInside( const typename Codim<0>::EntitySeed& seed, const int i ) const
1183  {
1184  return ( grid_.cell_facetag ) ? cartesianIndexInInside( seed, i ) : i;
1185  }
1186 
1187  int cartesianIndexInInside( const typename Codim<0>::EntitySeed& seed, const int i ) const
1188  {
1189  assert( i>= 0 && i<subEntities( seed, 1 ) );
1190  return grid_.cell_facetag[ grid_.cell_facepos[ seed.index() ] + i ] ;
1191  }
1192 
1193  typename Codim<0>::EntitySeed
1194  neighbor( const typename Codim<0>::EntitySeed& seed, const int i ) const
1195  {
1196  const int face = this->template subEntitySeed<1>( seed, i ).index();
1197  int nb = grid_.face_cells[ 2 * face ];
1198  if( nb == seed.index() )
1199  {
1200  nb = grid_.face_cells[ 2 * face + 1 ];
1201  }
1202 
1203  typedef typename Codim<0>::EntitySeed EntitySeed;
1204  return EntitySeed( nb );
1205  }
1206 
1207  int
1208  indexInOutside( const typename Codim<0>::EntitySeed& seed, const int i ) const
1209  {
1210  if( grid_.cell_facetag )
1211  {
1212  // if cell_facetag is present we assume pseudo Cartesian corner point case
1213  const int in_inside = cartesianIndexInInside( seed, i );
1214  return in_inside + ((in_inside % 2) ? -1 : 1);
1215  }
1216  else
1217  {
1218  typedef typename Codim<0>::EntitySeed EntitySeed;
1219  EntitySeed nb = neighbor( seed, i );
1220  const int faces = subEntities( seed, 1 );
1221  for( int face = 0; face<faces; ++ face )
1222  {
1223  if( neighbor( nb, face ).equals(seed) )
1224  {
1225  return indexInInside( nb, face );
1226  }
1227  }
1228  DUNE_THROW(InvalidStateException,"inverse intersection not found");
1229  return -1;
1230  }
1231  }
1232 
1233  template <class EntitySeed>
1234  GlobalCoordinate
1235  outerNormal( const EntitySeed& seed, const int i ) const
1236  {
1237  const int face = this->template subEntitySeed<1>( seed, i ).index();
1238  const int normalIdx = face * GlobalCoordinate :: dimension ;
1239  GlobalCoordinate normal = copyToGlobalCoordinate( grid_.face_normals + normalIdx );
1240  const int nb = grid_.face_cells[ 2*face ];
1241  if( nb != seed.index() )
1242  {
1243  normal *= -1.0;
1244  }
1245  return normal;
1246  }
1247 
1248  template <class EntitySeed>
1249  GlobalCoordinate
1250  unitOuterNormal( const EntitySeed& seed, const int i ) const
1251  {
1252  const int face = this->template subEntitySeed<1>( seed, i ).index();
1253  if( seed.index() == grid_.face_cells[ 2*face ] )
1254  {
1255  return unitOuterNormals_[ face ];
1256  }
1257  else
1258  {
1259  GlobalCoordinate normal = unitOuterNormals_[ face ];
1260  normal *= -1.0;
1261  return normal;
1262  }
1263  }
1264 
1265  template <class EntitySeed>
1266  GlobalCoordinate centroids( const EntitySeed& seed ) const
1267  {
1268  if( ! seed.isValid() )
1269  return GlobalCoordinate( 0 );
1270 
1271  const int index = GlobalCoordinate :: dimension * seed.index();
1272  const int codim = EntitySeed::codimension;
1273  assert( index >= 0 && index < size( codim ) * GlobalCoordinate :: dimension );
1274 
1275  if( codim == 0 )
1276  {
1277  return copyToGlobalCoordinate( grid_.cell_centroids + index );
1278  }
1279  else if ( codim == 1 )
1280  {
1281  return copyToGlobalCoordinate( grid_.face_centroids + index );
1282  }
1283  else if( codim == dim )
1284  {
1285  return copyToGlobalCoordinate( grid_.node_coordinates + index );
1286  }
1287  else
1288  {
1289  DUNE_THROW(InvalidStateException,"codimension not implemented");
1290  return GlobalCoordinate( 0 );
1291  }
1292  }
1293 
1294  GlobalCoordinate copyToGlobalCoordinate( const double* coords ) const
1295  {
1296  GlobalCoordinate coordinate;
1297  for( int i=0; i<GlobalCoordinate::dimension; ++i )
1298  {
1299  coordinate[ i ] = coords[ i ];
1300  }
1301  return coordinate;
1302  }
1303 
1304  template <class EntitySeed>
1305  double volumes( const EntitySeed& seed ) const
1306  {
1307  static const int codim = EntitySeed::codimension;
1308  if( codim == dim || ! seed.isValid() )
1309  {
1310  return 1.0;
1311  }
1312  else
1313  {
1314  assert( seed.isValid() );
1315 
1316  if( codim == 0 )
1317  {
1318  return grid_.cell_volumes[ seed.index() ];
1319  }
1320  else if ( codim == 1 )
1321  {
1322  return grid_.face_areas[ seed.index() ];
1323  }
1324  else
1325  {
1326  DUNE_THROW(InvalidStateException,"codimension not implemented");
1327  return 0.0;
1328  }
1329  }
1330  }
1331 
1332  protected:
1333  void init ()
1334  {
1335  // copy Cartesian dimensions
1336  for( int i=0; i<3; ++i )
1337  {
1338  cartDims_[ i ] = grid_.cartdims[ i ];
1339  }
1340 
1341  // setup list of cell vertices
1342  const int numCells = size( 0 );
1343 
1344  cellVertices_.resize( numCells );
1345 
1346  // sort vertices such that they comply with the dune reference cube
1347  if( grid_.cell_facetag )
1348  {
1349  typedef std::array<int, 3> KeyType;
1350  std::map< const KeyType, const int > vertexFaceTags;
1351  const int vertexFacePattern [8][3] = {
1352  { 0, 2, 4 }, // vertex 0
1353  { 1, 2, 4 }, // vertex 1
1354  { 0, 3, 4 }, // vertex 2
1355  { 1, 3, 4 }, // vertex 3
1356  { 0, 2, 5 }, // vertex 4
1357  { 1, 2, 5 }, // vertex 5
1358  { 0, 3, 5 }, // vertex 6
1359  { 1, 3, 5 } // vertex 7
1360  };
1361 
1362  for( int i=0; i<8; ++i )
1363  {
1364  KeyType key; key.fill( 4 ); // default is 4 which is the first z coord (for the 2d case)
1365  for( int j=0; j<dim; ++j )
1366  {
1367  key[ j ] = vertexFacePattern[ i ][ j ];
1368  }
1369 
1370  vertexFaceTags.insert( std::make_pair( key, i ) );
1371  }
1372 
1373  for (int c = 0; c < numCells; ++c)
1374  {
1375  if( dim == 2 )
1376  {
1377  // for 2d Cartesian grids the face ordering is wrong
1378  int f = grid_.cell_facepos[ c ];
1379  std::swap( grid_.cell_faces[ f+1 ], grid_.cell_faces[ f+2 ] );
1380  std::swap( grid_.cell_facetag[ f+1 ], grid_.cell_facetag[ f+2 ] );
1381  }
1382 
1383  typedef std::map<int,int> vertexmap_t;
1384  typedef typename vertexmap_t :: iterator iterator;
1385 
1386  std::vector< vertexmap_t > cell_pts( dim*2 );
1387 
1388  for (int hf=grid_.cell_facepos[ c ]; hf < grid_.cell_facepos[c+1]; ++hf)
1389  {
1390  const int f = grid_.cell_faces[ hf ];
1391  const int faceTag = grid_.cell_facetag[ hf ];
1392 
1393  for( int nodepos=grid_.face_nodepos[f]; nodepos<grid_.face_nodepos[f+1]; ++nodepos )
1394  {
1395  const int node = grid_.face_nodes[ nodepos ];
1396  iterator it = cell_pts[ faceTag ].find( node );
1397  if( it == cell_pts[ faceTag ].end() )
1398  {
1399  cell_pts[ faceTag ].insert( std::make_pair( node, 1 ) );
1400  }
1401  else
1402  {
1403  // increase vertex reference counter
1404  (*it).second++;
1405  }
1406  }
1407  }
1408 
1409  typedef std::map< int, std::set<int> > vertexlist_t;
1410  vertexlist_t vertexList;
1411 
1412  for( int faceTag = 0; faceTag<dim*2; ++faceTag )
1413  {
1414  for( iterator it = cell_pts[ faceTag ].begin(),
1415  end = cell_pts[ faceTag ].end(); it != end; ++it )
1416  {
1417 
1418  // only consider vertices with one appearance
1419  if( (*it).second == 1 )
1420  {
1421  vertexList[ (*it).first ].insert( faceTag );
1422  }
1423  }
1424  }
1425 
1426  assert( int(vertexList.size()) == ( dim == 2 ? 4 : 8) );
1427 
1428  cellVertices_[ c ].resize( vertexList.size() );
1429  for( auto it = vertexList.begin(), end = vertexList.end(); it != end; ++it )
1430  {
1431  assert( (*it).second.size() == dim );
1432  KeyType key; key.fill( 4 ); // fill with 4 which is the first z coord
1433 
1434  std::copy( (*it).second.begin(), (*it).second.end(), key.begin() );
1435  auto vx = vertexFaceTags.find( key );
1436  assert( vx != vertexFaceTags.end() );
1437  if( vx != vertexFaceTags.end() )
1438  {
1439  if( (*vx).second >= int(cellVertices_[ c ].size()) )
1440  cellVertices_[ c ].resize( (*vx).second+1 );
1441  // store node number on correct local position
1442  cellVertices_[ c ][ (*vx).second ] = (*it).first ;
1443  }
1444  }
1445  }
1446 
1447  // if face_tag is available we assume that the elements follow a cube-like structure
1448  geomTypes_.resize(dim + 1);
1449  GeometryType tmp;
1450  for (int codim = 0; codim <= dim; ++codim)
1451  {
1452  tmp = Dune::GeometryTypes::cube(dim - codim);
1453  geomTypes_[codim].push_back(tmp);
1454  }
1455  }
1456  else // if ( grid_.cell_facetag )
1457  {
1458  int maxVx = 0 ;
1459  int minVx = std::numeric_limits<int>::max();
1460 
1461  for (int c = 0; c < numCells; ++c)
1462  {
1463  std::set<int> cell_pts;
1464  for (int hf=grid_.cell_facepos[ c ]; hf < grid_.cell_facepos[c+1]; ++hf)
1465  {
1466  int f = grid_.cell_faces[ hf ];
1467  const int* fnbeg = grid_.face_nodes + grid_.face_nodepos[f];
1468  const int* fnend = grid_.face_nodes + grid_.face_nodepos[f+1];
1469  cell_pts.insert(fnbeg, fnend);
1470  }
1471 
1472  cellVertices_[ c ].resize( cell_pts.size() );
1473  std::copy(cell_pts.begin(), cell_pts.end(), cellVertices_[ c ].begin() );
1474  maxVx = std::max( maxVx, int( cell_pts.size() ) );
1475  minVx = std::min( minVx, int( cell_pts.size() ) );
1476  }
1477 
1478  if( minVx == maxVx && maxVx == 4 )
1479  {
1480  for (int c = 0; c < numCells; ++c)
1481  {
1482  assert( cellVertices_[ c ].size() == 4 );
1483  GlobalCoordinate center( 0 );
1484  GlobalCoordinate p[ dim+1 ];
1485  for( int i=0; i<dim+1; ++i )
1486  {
1487  const int vertex = cellVertices_[ c ][ i ];
1488 
1489  for( int d=0; d<dim; ++d )
1490  {
1491  center[ d ] += grid_.node_coordinates[ vertex*dim + d ];
1492  p[ i ][ d ] = grid_.node_coordinates[ vertex*dim + d ];
1493  }
1494  }
1495  center *= 0.25;
1496  for( int d=0; d<dim; ++d )
1497  {
1498  grid_.cell_centroids[ c*dim + d ] = center[ d ];
1499  }
1500 
1501  Dune::GeometryType simplex;
1502  simplex = Dune::GeometryTypes::simplex(dim);
1503 
1504  typedef Dune::AffineGeometry< ctype, dim, dimworld> AffineGeometryType;
1505  AffineGeometryType geometry( simplex, p );
1506  grid_.cell_volumes[ c ] = geometry.volume();
1507  }
1508  }
1509 
1510  // check face normals
1511  {
1512  const int faces = grid_.number_of_faces;
1513  for( int face = 0 ; face < faces; ++face )
1514  {
1515  const int a = grid_.face_cells[ 2*face ];
1516  const int b = grid_.face_cells[ 2*face + 1 ];
1517 
1518  assert( a >=0 || b >=0 );
1519 
1520  if( grid_.face_areas[ face ] < 0 )
1521  std::abort();
1522 
1523  GlobalCoordinate centerDiff( 0 );
1524  if( b >= 0 )
1525  {
1526  for( int d=0; d<dimworld; ++d )
1527  {
1528  centerDiff[ d ] = grid_.cell_centroids[ b*dimworld + d ];
1529  }
1530  }
1531  else
1532  {
1533  for( int d=0; d<dimworld; ++d )
1534  {
1535  centerDiff[ d ] = grid_.face_centroids[ face*dimworld + d ];
1536  }
1537  }
1538 
1539  if( a >= 0 )
1540  {
1541  for( int d=0; d<dimworld; ++d )
1542  {
1543  centerDiff[ d ] -= grid_.cell_centroids[ a*dimworld + d ];
1544  }
1545  }
1546  else
1547  {
1548  for( int d=0; d<dimworld; ++d )
1549  {
1550  centerDiff[ d ] -= grid_.face_centroids[ face*dimworld + d ];
1551  }
1552  }
1553 
1554  GlobalCoordinate normal( 0 );
1555  for( int d=0; d<dimworld; ++d )
1556  {
1557  normal[ d ] = grid_.face_normals[ face*dimworld + d ];
1558  }
1559 
1560  if( centerDiff.two_norm() < 1e-10 )
1561  std::abort();
1562 
1563  // if diff and normal point in different direction, flip faces
1564  if( centerDiff * normal < 0 )
1565  {
1566  grid_.face_cells[ 2*face ] = b;
1567  grid_.face_cells[ 2*face + 1 ] = a;
1568  }
1569  }
1570  }
1571 
1572  bool allSimplex = true ;
1573  bool allCube = true ;
1574 
1575  for (int c = 0; c < numCells; ++c)
1576  {
1577  const int nVx = cellVertices_[ c ].size();
1578  if( nVx != 4 )
1579  {
1580  allSimplex = false;
1581  }
1582 
1583  if( nVx != 8 )
1584  {
1585  allCube = false;
1586  }
1587  }
1588  // Propogate the cell geometry type to all codimensions
1589  geomTypes_.resize(dim + 1);
1590  GeometryType tmp;
1591  for (int codim = 0; codim <= dim; ++codim)
1592  {
1593  if( allSimplex )
1594  {
1595  tmp = Dune::GeometryTypes::simplex(dim - codim);
1596  geomTypes_[ codim ].push_back( tmp );
1597  }
1598  else if ( allCube )
1599  {
1600  tmp = Dune::GeometryTypes::cube(dim - codim);
1601  geomTypes_[ codim ].push_back( tmp );
1602  }
1603  else
1604  {
1605  tmp = Dune::GeometryTypes::none(dim - codim);
1606  geomTypes_[ codim ].push_back( tmp );
1607  }
1608  }
1609 
1610  } // end else of ( grid_.cell_facetag )
1611 
1612  nBndSegments_ = 0;
1613  unitOuterNormals_.resize( grid_.number_of_faces );
1614  for( int face = 0; face < grid_.number_of_faces; ++face )
1615  {
1616  const int normalIdx = face * GlobalCoordinate :: dimension ;
1617  GlobalCoordinate normal = copyToGlobalCoordinate( grid_.face_normals + normalIdx );
1618  normal /= normal.two_norm();
1619  unitOuterNormals_[ face ] = normal;
1620 
1621  if( isBoundaryFace( face ) )
1622  {
1623  // increase number if boundary segments
1624  ++nBndSegments_;
1625  const int facePos = 2 * face ;
1626  // store negative number to indicate boundary
1627  // the abstract value is the segment index
1628  if( grid_.face_cells[ facePos ] < 0 )
1629  {
1630  grid_.face_cells[ facePos ] = -nBndSegments_;
1631  }
1632  else if ( grid_.face_cells[ facePos+1 ] < 0 )
1633  {
1634  grid_.face_cells[ facePos+1 ] = -nBndSegments_;
1635  }
1636  }
1637  }
1638  }
1639 
1640  void print( std::ostream& out, const UnstructuredGridType& grid ) const
1641  {
1642  const int numCells = grid.number_of_cells;
1643  for( int c=0; c<numCells; ++c )
1644  {
1645  out << "cell " << c << " : faces = " << std::endl;
1646  for (int hf=grid.cell_facepos[ c ]; hf < grid.cell_facepos[c+1]; ++hf)
1647  {
1648  int f = grid_.cell_faces[ hf ];
1649  const int* fnbeg = grid_.face_nodes + grid_.face_nodepos[f];
1650  const int* fnend = grid_.face_nodes + grid_.face_nodepos[f+1];
1651  out << f << " vx = " ;
1652  while( fnbeg != fnend )
1653  {
1654  out << *fnbeg << " ";
1655  ++fnbeg;
1656  }
1657  out << std::endl;
1658  }
1659  out << std::endl;
1660 
1661  const auto& vx = cellVertices_[ c ];
1662  out << "cell " << c << " : vertices = ";
1663  for( size_t i=0; i<vx.size(); ++i )
1664  out << vx[ i ] << " ";
1665  out << std::endl;
1666  }
1667 
1668  }
1669 
1670  protected:
1671  UnstructuredGridPtr gridPtr_;
1672  const UnstructuredGridType& grid_;
1673 
1674  CommunicationType comm_;
1675  std::array< int, 3 > cartDims_;
1676  std::vector< std::vector< GeometryType > > geomTypes_;
1677  std::vector< std::vector< int > > cellVertices_;
1678 
1679  std::vector< GlobalCoordinate > unitOuterNormals_;
1680 
1681  mutable LeafIndexSet leafIndexSet_;
1682  mutable GlobalIdSet globalIdSet_;
1683  mutable LocalIdSet localIdSet_;
1684 
1685  size_t nBndSegments_;
1686 
1687  private:
1688  // no copying
1689  PolyhedralGrid ( const PolyhedralGrid& );
1690  };
1691 
1692 
1693 
1694  // PolyhedralGrid::Codim
1695  // -------------
1696 
1697  template< int dim, int dimworld, typename coord_t >
1698  template< int codim >
1700  : public Base::template Codim< codim >
1701  {
1710  typedef typename Traits::template Codim< codim >::Entity Entity;
1711 
1716  typedef typename Traits::template Codim< codim >::EntityPointer EntityPointer;
1717 
1731  typedef typename Traits::template Codim< codim >::Geometry Geometry;
1732 
1741  typedef typename Traits::template Codim< codim >::LocalGeometry LocalGeometry;
1742 
1748  template< PartitionIteratorType pitype >
1749  struct Partition
1750  {
1751  typedef typename Traits::template Codim< codim >
1752  ::template Partition< pitype >::LeafIterator
1753  LeafIterator;
1754  typedef typename Traits::template Codim< codim >
1755  ::template Partition< pitype >::LevelIterator
1756  LevelIterator;
1757  };
1758 
1766  typedef typename Partition< All_Partition >::LeafIterator LeafIterator;
1767 
1775  typedef typename Partition< All_Partition >::LevelIterator LevelIterator;
1776 
1778  };
1779 
1780 } // namespace Dune
1781 
1782 #include <opm/grid/polyhedralgrid/persistentcontainer.hh>
1783 #include <opm/grid/polyhedralgrid/cartesianindexmapper.hh>
1784 #include <opm/grid/polyhedralgrid/gridhelpers.hh>
1785 
1786 #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:250
Definition: idset.hh:17
Definition: indexset.hh:25
Definition: intersectioniterator.hh:16
Definition: intersection.hh:20
Definition: iterator.hh:21
Definition: geometry.hh:269
identical grid wrapper
Definition: grid.hh:163
void postAdapt()
Definition: grid.hh:614
Traits::ctype ctype
type of vector coordinates (e.g., double)
Definition: grid.hh:312
int ghostSize(int, int codim) const
obtain size of ghost region for a grid level
Definition: grid.hh:654
bool loadBalance(CommDataHandleIF< DataHandle, Data > &)
rebalance the load each process has to handle
Definition: grid.hh:758
Traits::template Codim< EntitySeed::codimension >::EntityPointer entityPointer(const EntitySeed &seed) const
obtain EntityPointer from EntitySeed.
Definition: grid.hh:818
const CommunicationType & comm() const
obtain CollectiveCommunication object
Definition: grid.hh:721
Traits::template Codim< EntitySeed::codimension >::Entity entity(const EntitySeed &seed) const
obtain EntityPointer from EntitySeed.
Definition: grid.hh:829
bool preAdapt()
Definition: grid.hh:592
typename Traits::CollectiveCommunication CollectiveCommunication
communicator with all other processes having some part of the grid
Definition: grid.hh:319
int ghostSize(int codim) const
obtain size of ghost region for the leaf grid
Definition: grid.hh:634
bool loadBalance(DofManager &)
rebalance the load each process has to handle
Definition: grid.hh:778
Traits::GlobalIdSet GlobalIdSet
type of global id set
Definition: grid.hh:287
int maxLevel() const
obtain maximal grid level
Definition: grid.hh:424
Traits::LevelIndexSet LevelIndexSet
type of level index set
Definition: grid.hh:275
LevelGridView levelGridView(int) const
View for a grid level for All_Partition.
Definition: grid.hh:802
PolyhedralGrid(const std::vector< int > &n, const std::vector< double > &dx)
constructor
Definition: grid.hh:356
Partition< pitype >::LeafGridView leafGridView() const
View for the leaf grid.
Definition: grid.hh:794
Traits::HierarchicIterator HierarchicIterator
iterator over the grid hierarchy
Definition: grid.hh:229
void scatterData([[maybe_unused]] DataHandle &handle) const
Moves data from the global (all data on process) view to the distributed view.
Definition: grid.hh:893
int size(int, int codim) const
obtain number of entites on a level
Definition: grid.hh:437
Partition< pitype >::LevelGridView levelGridView(int) const
View for a grid level.
Definition: grid.hh:785
Partition< All_Partition >::LevelGridView LevelGridView
View types for All_Partition.
Definition: grid.hh:249
void update()
update grid caches
Definition: grid.hh:848
void switchToDistributedView()
Switch to the distributed view.
Definition: grid.hh:708
int overlapSize(int, int) const
obtain size of overlap region for a grid level
Definition: grid.hh:644
int size(int codim) const
obtain number of leaf entities
Definition: grid.hh:448
Traits::LocalIdSet LocalIdSet
type of local id set
Definition: grid.hh:304
int overlapSize(int) const
obtain size of overlap region for the leaf grid
Definition: grid.hh:625
Traits::LevelIntersectionIterator LevelIntersectionIterator
iterator over intersections with other entities on the same level
Definition: grid.hh:233
GridFamily::Traits Traits
type of the grid traits
Definition: grid.hh:212
PolyhedralGrid(UnstructuredGridPtr &&gridPtr)
constructor
Definition: grid.hh:375
Traits::LeafIntersectionIterator LeafIntersectionIterator
iterator over intersections with other entities on the leaf level
Definition: grid.hh:231
bool adapt(DataHandle &)
Definition: grid.hh:608
int size(int, GeometryType type) const
obtain number of entites on a level
Definition: grid.hh:477
void communicate(DataHandle &, InterfaceType, CommunicationDirection) const
communicate information on leaf entities
Definition: grid.hh:694
size_t numBoundarySegments() const
obtain number of leaf entities
Definition: grid.hh:497
void switchToGlobalView()
Switch to the global view.
Definition: grid.hh:702
bool loadBalance()
rebalance the load each process has to handle
Definition: grid.hh:737
PolyhedralGrid(const UnstructuredGridType &grid)
constructor
Definition: grid.hh:394
LeafGridView leafGridView() const
View for the leaf grid for All_Partition.
Definition: grid.hh:809
int size(GeometryType type) const
returns the number of boundary segments within the macro grid
Definition: grid.hh:486
void communicate(DataHandle &, InterfaceType, CommunicationDirection, int) const
communicate information on a grid level
Definition: grid.hh:673
bool adapt()
Definition: grid.hh:598
Traits::LeafIndexSet LeafIndexSet
type of leaf index set
Definition: grid.hh:265
Transform a corner-point grid ZCORN field to account for MINPV processing.
Definition: MinpvProcessor.hpp:37
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:1701
Partition< All_Partition >::LevelIterator LevelIterator
type of leaf iterator
Definition: grid.hh:1775
Traits::template Codim< codim >::Entity Entity
type of entity
Definition: grid.hh:1710
Traits::template Codim< codim >::LocalGeometry LocalGeometry
type of local geometry
Definition: grid.hh:1741
Traits::template Codim< codim >::EntityPointer EntityPointer
type of entity pointer
Definition: grid.hh:1716
Partition< All_Partition >::LeafIterator LeafIterator
type of level iterator
Definition: grid.hh:1766
Traits::template Codim< codim >::Geometry Geometry
type of world geometry
Definition: grid.hh:1731
Types for GridView.
Definition: grid.hh:243
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