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  assert(!geomTypes( Seed::codimension ).empty());
1146  return geomTypes( Seed::codimension )[ 0 ];
1147  }
1148  else
1149  {
1150  // 3d faces
1151  if( dim == 3 && Seed::codimension == 1 )
1152  {
1153  GeometryType face;
1154  const int nVx = corners( seed );
1155  if( nVx == 4 ) // quad face
1156  face = Dune::GeometryTypes::cube(2);
1157  else if( nVx == 3 ) // triangle face
1158  face = Dune::GeometryTypes::simplex(2);
1159  else // polygonal face
1160  face = Dune::GeometryTypes::none(2);
1161  return face;
1162  }
1163 
1164  // for codim 2 and codim 3 there is only one option
1165  assert(!geomTypes( Seed::codimension ).empty());
1166  return geomTypes( Seed::codimension )[ 0 ];
1167  }
1168  }
1169 
1170  int indexInInside( const typename Codim<0>::EntitySeed& seed, const int i ) const
1171  {
1172  return ( grid_.cell_facetag ) ? cartesianIndexInInside( seed, i ) : i;
1173  }
1174 
1175  int cartesianIndexInInside( const typename Codim<0>::EntitySeed& seed, const int i ) const
1176  {
1177  assert( i>= 0 && i<subEntities( seed, 1 ) );
1178  return grid_.cell_facetag[ grid_.cell_facepos[ seed.index() ] + i ] ;
1179  }
1180 
1181  typename Codim<0>::EntitySeed
1182  neighbor( const typename Codim<0>::EntitySeed& seed, const int i ) const
1183  {
1184  const int face = this->template subEntitySeed<1>( seed, i ).index();
1185  int nb = grid_.face_cells[ 2 * face ];
1186  if( nb == seed.index() )
1187  {
1188  nb = grid_.face_cells[ 2 * face + 1 ];
1189  }
1190 
1191  typedef typename Codim<0>::EntitySeed EntitySeed;
1192  return EntitySeed( nb );
1193  }
1194 
1195  int
1196  indexInOutside( const typename Codim<0>::EntitySeed& seed, const int i ) const
1197  {
1198  if( grid_.cell_facetag )
1199  {
1200  // if cell_facetag is present we assume pseudo Cartesian corner point case
1201  const int in_inside = cartesianIndexInInside( seed, i );
1202  return in_inside + ((in_inside % 2) ? -1 : 1);
1203  }
1204  else
1205  {
1206  typedef typename Codim<0>::EntitySeed EntitySeed;
1207  EntitySeed nb = neighbor( seed, i );
1208  const int faces = subEntities( seed, 1 );
1209  for( int face = 0; face<faces; ++ face )
1210  {
1211  if( neighbor( nb, face ).equals(seed) )
1212  {
1213  return indexInInside( nb, face );
1214  }
1215  }
1216  DUNE_THROW(InvalidStateException,"inverse intersection not found");
1217  return -1;
1218  }
1219  }
1220 
1221  template <class EntitySeed>
1222  GlobalCoordinate
1223  outerNormal( const EntitySeed& seed, const int i ) const
1224  {
1225  const int face = this->template subEntitySeed<1>( seed, i ).index();
1226  const int normalIdx = face * GlobalCoordinate :: dimension ;
1227  GlobalCoordinate normal = copyToGlobalCoordinate( grid_.face_normals + normalIdx );
1228  const int nb = grid_.face_cells[ 2*face ];
1229  if( nb != seed.index() )
1230  {
1231  normal *= -1.0;
1232  }
1233  return normal;
1234  }
1235 
1236  template <class EntitySeed>
1237  GlobalCoordinate
1238  unitOuterNormal( const EntitySeed& seed, const int i ) const
1239  {
1240  const int face = this->template subEntitySeed<1>( seed, i ).index();
1241  if( seed.index() == grid_.face_cells[ 2*face ] )
1242  {
1243  return unitOuterNormals_[ face ];
1244  }
1245  else
1246  {
1247  GlobalCoordinate normal = unitOuterNormals_[ face ];
1248  normal *= -1.0;
1249  return normal;
1250  }
1251  }
1252 
1253  template <class EntitySeed>
1254  GlobalCoordinate centroids( const EntitySeed& seed ) const
1255  {
1256  if( ! seed.isValid() )
1257  return GlobalCoordinate( 0 );
1258 
1259  const int index = GlobalCoordinate :: dimension * seed.index();
1260  const int codim = EntitySeed::codimension;
1261  assert( index >= 0 && index < size( codim ) * GlobalCoordinate :: dimension );
1262 
1263  if( codim == 0 )
1264  {
1265  return copyToGlobalCoordinate( grid_.cell_centroids + index );
1266  }
1267  else if ( codim == 1 )
1268  {
1269  return copyToGlobalCoordinate( grid_.face_centroids + index );
1270  }
1271  else if( codim == dim )
1272  {
1273  return copyToGlobalCoordinate( grid_.node_coordinates + index );
1274  }
1275  else
1276  {
1277  DUNE_THROW(InvalidStateException,"codimension not implemented");
1278  return GlobalCoordinate( 0 );
1279  }
1280  }
1281 
1282  GlobalCoordinate copyToGlobalCoordinate( const double* coords ) const
1283  {
1284  GlobalCoordinate coordinate;
1285  for( int i=0; i<GlobalCoordinate::dimension; ++i )
1286  {
1287  coordinate[ i ] = coords[ i ];
1288  }
1289  return coordinate;
1290  }
1291 
1292  template <class EntitySeed>
1293  double volumes( const EntitySeed& seed ) const
1294  {
1295  static const int codim = EntitySeed::codimension;
1296  if( codim == dim || ! seed.isValid() )
1297  {
1298  return 1.0;
1299  }
1300  else
1301  {
1302  assert( seed.isValid() );
1303 
1304  if( codim == 0 )
1305  {
1306  return grid_.cell_volumes[ seed.index() ];
1307  }
1308  else if ( codim == 1 )
1309  {
1310  return grid_.face_areas[ seed.index() ];
1311  }
1312  else
1313  {
1314  DUNE_THROW(InvalidStateException,"codimension not implemented");
1315  return 0.0;
1316  }
1317  }
1318  }
1319 
1320  protected:
1321  void init ()
1322  {
1323  // copy Cartesian dimensions
1324  for( int i=0; i<3; ++i )
1325  {
1326  cartDims_[ i ] = grid_.cartdims[ i ];
1327  }
1328 
1329  // setup list of cell vertices
1330  const int numCells = size( 0 );
1331 
1332  cellVertices_.resize( numCells );
1333 
1334  // sort vertices such that they comply with the dune reference cube
1335  if( grid_.cell_facetag )
1336  {
1337  typedef std::array<int, 3> KeyType;
1338  std::map< const KeyType, const int > vertexFaceTags;
1339  const int vertexFacePattern [8][3] = {
1340  { 0, 2, 4 }, // vertex 0
1341  { 1, 2, 4 }, // vertex 1
1342  { 0, 3, 4 }, // vertex 2
1343  { 1, 3, 4 }, // vertex 3
1344  { 0, 2, 5 }, // vertex 4
1345  { 1, 2, 5 }, // vertex 5
1346  { 0, 3, 5 }, // vertex 6
1347  { 1, 3, 5 } // vertex 7
1348  };
1349 
1350  for( int i=0; i<8; ++i )
1351  {
1352  KeyType key; key.fill( 4 ); // default is 4 which is the first z coord (for the 2d case)
1353  for( int j=0; j<dim; ++j )
1354  {
1355  key[ j ] = vertexFacePattern[ i ][ j ];
1356  }
1357 
1358  vertexFaceTags.insert( std::make_pair( key, i ) );
1359  }
1360 
1361  for (int c = 0; c < numCells; ++c)
1362  {
1363  if( dim == 2 )
1364  {
1365  // for 2d Cartesian grids the face ordering is wrong
1366  int f = grid_.cell_facepos[ c ];
1367  std::swap( grid_.cell_faces[ f+1 ], grid_.cell_faces[ f+2 ] );
1368  std::swap( grid_.cell_facetag[ f+1 ], grid_.cell_facetag[ f+2 ] );
1369  }
1370 
1371  typedef std::map<int,int> vertexmap_t;
1372  typedef typename vertexmap_t :: iterator iterator;
1373 
1374  std::vector< vertexmap_t > cell_pts( dim*2 );
1375 
1376  for (int hf=grid_.cell_facepos[ c ]; hf < grid_.cell_facepos[c+1]; ++hf)
1377  {
1378  const int f = grid_.cell_faces[ hf ];
1379  const int faceTag = grid_.cell_facetag[ hf ];
1380 
1381  for( int nodepos=grid_.face_nodepos[f]; nodepos<grid_.face_nodepos[f+1]; ++nodepos )
1382  {
1383  const int node = grid_.face_nodes[ nodepos ];
1384  iterator it = cell_pts[ faceTag ].find( node );
1385  if( it == cell_pts[ faceTag ].end() )
1386  {
1387  cell_pts[ faceTag ].insert( std::make_pair( node, 1 ) );
1388  }
1389  else
1390  {
1391  // increase vertex reference counter
1392  (*it).second++;
1393  }
1394  }
1395  }
1396 
1397  typedef std::map< int, std::set<int> > vertexlist_t;
1398  vertexlist_t vertexList;
1399 
1400  for( int faceTag = 0; faceTag<dim*2; ++faceTag )
1401  {
1402  for( iterator it = cell_pts[ faceTag ].begin(),
1403  end = cell_pts[ faceTag ].end(); it != end; ++it )
1404  {
1405 
1406  // only consider vertices with one appearance
1407  if( (*it).second == 1 )
1408  {
1409  vertexList[ (*it).first ].insert( faceTag );
1410  }
1411  }
1412  }
1413 
1414  assert( int(vertexList.size()) == ( dim == 2 ? 4 : 8) );
1415 
1416  cellVertices_[ c ].resize( vertexList.size() );
1417  for( auto it = vertexList.begin(), end = vertexList.end(); it != end; ++it )
1418  {
1419  assert( (*it).second.size() == dim );
1420  KeyType key; key.fill( 4 ); // fill with 4 which is the first z coord
1421 
1422  std::copy( (*it).second.begin(), (*it).second.end(), key.begin() );
1423  auto vx = vertexFaceTags.find( key );
1424  assert( vx != vertexFaceTags.end() );
1425  if( vx != vertexFaceTags.end() )
1426  {
1427  if( (*vx).second >= int(cellVertices_[ c ].size()) )
1428  cellVertices_[ c ].resize( (*vx).second+1 );
1429  // store node number on correct local position
1430  cellVertices_[ c ][ (*vx).second ] = (*it).first ;
1431  }
1432  }
1433  }
1434 
1435  // if face_tag is available we assume that the elements follow a cube-like structure
1436  geomTypes_.resize(dim + 1);
1437  GeometryType tmp;
1438  for (int codim = 0; codim <= dim; ++codim)
1439  {
1440  tmp = Dune::GeometryTypes::cube(dim - codim);
1441  geomTypes_[codim].push_back(tmp);
1442  }
1443  }
1444  else // if ( grid_.cell_facetag )
1445  {
1446  int maxVx = 0 ;
1447  int minVx = std::numeric_limits<int>::max();
1448 
1449  for (int c = 0; c < numCells; ++c)
1450  {
1451  std::set<int> cell_pts;
1452  for (int hf=grid_.cell_facepos[ c ]; hf < grid_.cell_facepos[c+1]; ++hf)
1453  {
1454  int f = grid_.cell_faces[ hf ];
1455  const int* fnbeg = grid_.face_nodes + grid_.face_nodepos[f];
1456  const int* fnend = grid_.face_nodes + grid_.face_nodepos[f+1];
1457  cell_pts.insert(fnbeg, fnend);
1458  }
1459 
1460  cellVertices_[ c ].resize( cell_pts.size() );
1461  std::copy(cell_pts.begin(), cell_pts.end(), cellVertices_[ c ].begin() );
1462  maxVx = std::max( maxVx, int( cell_pts.size() ) );
1463  minVx = std::min( minVx, int( cell_pts.size() ) );
1464  }
1465 
1466  if( minVx == maxVx && maxVx == 4 )
1467  {
1468  for (int c = 0; c < numCells; ++c)
1469  {
1470  assert( cellVertices_[ c ].size() == 4 );
1471  GlobalCoordinate center( 0 );
1472  GlobalCoordinate p[ dim+1 ];
1473  for( int i=0; i<dim+1; ++i )
1474  {
1475  const int vertex = cellVertices_[ c ][ i ];
1476 
1477  for( int d=0; d<dim; ++d )
1478  {
1479  center[ d ] += grid_.node_coordinates[ vertex*dim + d ];
1480  p[ i ][ d ] = grid_.node_coordinates[ vertex*dim + d ];
1481  }
1482  }
1483  center *= 0.25;
1484  for( int d=0; d<dim; ++d )
1485  {
1486  grid_.cell_centroids[ c*dim + d ] = center[ d ];
1487  }
1488 
1489  Dune::GeometryType simplex;
1490  simplex = Dune::GeometryTypes::simplex(dim);
1491 
1492  typedef Dune::AffineGeometry< ctype, dim, dimworld> AffineGeometryType;
1493  AffineGeometryType geometry( simplex, p );
1494  grid_.cell_volumes[ c ] = geometry.volume();
1495  }
1496  }
1497 
1498  // check face normals
1499  {
1500  const int faces = grid_.number_of_faces;
1501  for( int face = 0 ; face < faces; ++face )
1502  {
1503  const int a = grid_.face_cells[ 2*face ];
1504  const int b = grid_.face_cells[ 2*face + 1 ];
1505 
1506  assert( a >=0 || b >=0 );
1507 
1508  if( grid_.face_areas[ face ] < 0 )
1509  std::abort();
1510 
1511  GlobalCoordinate centerDiff( 0 );
1512  if( b >= 0 )
1513  {
1514  for( int d=0; d<dimworld; ++d )
1515  {
1516  centerDiff[ d ] = grid_.cell_centroids[ b*dimworld + d ];
1517  }
1518  }
1519  else
1520  {
1521  for( int d=0; d<dimworld; ++d )
1522  {
1523  centerDiff[ d ] = grid_.face_centroids[ face*dimworld + d ];
1524  }
1525  }
1526 
1527  if( a >= 0 )
1528  {
1529  for( int d=0; d<dimworld; ++d )
1530  {
1531  centerDiff[ d ] -= grid_.cell_centroids[ a*dimworld + d ];
1532  }
1533  }
1534  else
1535  {
1536  for( int d=0; d<dimworld; ++d )
1537  {
1538  centerDiff[ d ] -= grid_.face_centroids[ face*dimworld + d ];
1539  }
1540  }
1541 
1542  GlobalCoordinate normal( 0 );
1543  for( int d=0; d<dimworld; ++d )
1544  {
1545  normal[ d ] = grid_.face_normals[ face*dimworld + d ];
1546  }
1547 
1548  if( centerDiff.two_norm() < 1e-10 )
1549  std::abort();
1550 
1551  // if diff and normal point in different direction, flip faces
1552  if( centerDiff * normal < 0 )
1553  {
1554  grid_.face_cells[ 2*face ] = b;
1555  grid_.face_cells[ 2*face + 1 ] = a;
1556  }
1557  }
1558  }
1559 
1560  bool allSimplex = true ;
1561  bool allCube = true ;
1562 
1563  for (int c = 0; c < numCells; ++c)
1564  {
1565  const int nVx = cellVertices_[ c ].size();
1566  if( nVx != 4 )
1567  {
1568  allSimplex = false;
1569  }
1570 
1571  if( nVx != 8 )
1572  {
1573  allCube = false;
1574  }
1575  }
1576  // Propogate the cell geometry type to all codimensions
1577  geomTypes_.resize(dim + 1);
1578  GeometryType tmp;
1579  for (int codim = 0; codim <= dim; ++codim)
1580  {
1581  if( allSimplex )
1582  {
1583  tmp = Dune::GeometryTypes::simplex(dim - codim);
1584  geomTypes_[ codim ].push_back( tmp );
1585  }
1586  else if ( allCube )
1587  {
1588  tmp = Dune::GeometryTypes::cube(dim - codim);
1589  geomTypes_[ codim ].push_back( tmp );
1590  }
1591  else
1592  {
1593  tmp = Dune::GeometryTypes::none(dim - codim);
1594  geomTypes_[ codim ].push_back( tmp );
1595  }
1596  }
1597 
1598  } // end else of ( grid_.cell_facetag )
1599 
1600  nBndSegments_ = 0;
1601  unitOuterNormals_.resize( grid_.number_of_faces );
1602  for( int face = 0; face < grid_.number_of_faces; ++face )
1603  {
1604  const int normalIdx = face * GlobalCoordinate :: dimension ;
1605  GlobalCoordinate normal = copyToGlobalCoordinate( grid_.face_normals + normalIdx );
1606  normal /= normal.two_norm();
1607  unitOuterNormals_[ face ] = normal;
1608 
1609  if( isBoundaryFace( face ) )
1610  {
1611  // increase number if boundary segments
1612  ++nBndSegments_;
1613  const int facePos = 2 * face ;
1614  // store negative number to indicate boundary
1615  // the abstract value is the segment index
1616  if( grid_.face_cells[ facePos ] < 0 )
1617  {
1618  grid_.face_cells[ facePos ] = -nBndSegments_;
1619  }
1620  else if ( grid_.face_cells[ facePos+1 ] < 0 )
1621  {
1622  grid_.face_cells[ facePos+1 ] = -nBndSegments_;
1623  }
1624  }
1625  }
1626  }
1627 
1628  void print( std::ostream& out, const UnstructuredGridType& grid ) const
1629  {
1630  const int numCells = grid.number_of_cells;
1631  for( int c=0; c<numCells; ++c )
1632  {
1633  out << "cell " << c << " : faces = " << std::endl;
1634  for (int hf=grid.cell_facepos[ c ]; hf < grid.cell_facepos[c+1]; ++hf)
1635  {
1636  int f = grid_.cell_faces[ hf ];
1637  const int* fnbeg = grid_.face_nodes + grid_.face_nodepos[f];
1638  const int* fnend = grid_.face_nodes + grid_.face_nodepos[f+1];
1639  out << f << " vx = " ;
1640  while( fnbeg != fnend )
1641  {
1642  out << *fnbeg << " ";
1643  ++fnbeg;
1644  }
1645  out << std::endl;
1646  }
1647  out << std::endl;
1648 
1649  const auto& vx = cellVertices_[ c ];
1650  out << "cell " << c << " : vertices = ";
1651  for( size_t i=0; i<vx.size(); ++i )
1652  out << vx[ i ] << " ";
1653  out << std::endl;
1654  }
1655 
1656  }
1657 
1658  protected:
1659  UnstructuredGridPtr gridPtr_;
1660  const UnstructuredGridType& grid_;
1661 
1663  std::array< int, 3 > cartDims_;
1664  std::vector< std::vector< GeometryType > > geomTypes_;
1665  std::vector< std::vector< int > > cellVertices_;
1666 
1667  std::vector< GlobalCoordinate > unitOuterNormals_;
1668 
1669  mutable LeafIndexSet leafIndexSet_;
1670  mutable GlobalIdSet globalIdSet_;
1671  mutable LocalIdSet localIdSet_;
1672 
1673  size_t nBndSegments_;
1674 
1675  private:
1676  // no copying
1677  PolyhedralGrid ( const PolyhedralGrid& );
1678  };
1679 
1680 
1681 
1682  // PolyhedralGrid::Codim
1683  // -------------
1684 
1685  template< int dim, int dimworld, typename coord_t >
1686  template< int codim >
1688  : public Base::template Codim< codim >
1689  {
1698  typedef typename Traits::template Codim< codim >::Entity Entity;
1699 
1704  typedef typename Traits::template Codim< codim >::EntityPointer EntityPointer;
1705 
1719  typedef typename Traits::template Codim< codim >::Geometry Geometry;
1720 
1729  typedef typename Traits::template Codim< codim >::LocalGeometry LocalGeometry;
1730 
1736  template< PartitionIteratorType pitype >
1737  struct Partition
1738  {
1739  typedef typename Traits::template Codim< codim >
1740  ::template Partition< pitype >::LeafIterator
1741  LeafIterator;
1742  typedef typename Traits::template Codim< codim >
1743  ::template Partition< pitype >::LevelIterator
1744  LevelIterator;
1745  };
1746 
1754  typedef typename Partition< All_Partition >::LeafIterator LeafIterator;
1755 
1763  typedef typename Partition< All_Partition >::LevelIterator LevelIterator;
1764 
1766  };
1767 
1768 } // namespace Dune
1769 
1770 #include <opm/grid/polyhedralgrid/persistentcontainer.hh>
1771 #include <opm/grid/polyhedralgrid/cartesianindexmapper.hh>
1772 #include <opm/grid/polyhedralgrid/gridhelpers.hh>
1773 
1774 #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:1689
Partition< All_Partition >::LevelIterator LevelIterator
type of leaf iterator
Definition: grid.hh:1763
Traits::template Codim< codim >::Entity Entity
type of entity
Definition: grid.hh:1698
Traits::template Codim< codim >::LocalGeometry LocalGeometry
type of local geometry
Definition: grid.hh:1729
Traits::template Codim< codim >::EntityPointer EntityPointer
type of entity pointer
Definition: grid.hh:1704
Partition< All_Partition >::LeafIterator LeafIterator
type of level iterator
Definition: grid.hh:1754
Traits::template Codim< codim >::Geometry Geometry
type of world geometry
Definition: grid.hh:1719
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