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
46namespace 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 >
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 >
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 >
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
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 >
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
892 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 >
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
Partition< All_Partition >::LevelGridView LevelGridView
View types for All_Partition.
Definition: grid.hh:249
Traits::template Codim< EntitySeed::codimension >::EntityPointer entityPointer(const EntitySeed &seed) const
obtain EntityPointer from EntitySeed.
Definition: grid.hh:818
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
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
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
void scatterData(DataHandle &handle) const
Moves data from the global (all data on process) view to the distributed view.
Definition: grid.hh:893
Traits::LocalIdSet LocalIdSet
type of local id set
Definition: grid.hh:304
const CommunicationType & comm() const
obtain CollectiveCommunication object
Definition: grid.hh:721
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 >::LeafIterator LeafIterator
type of level iterator
Definition: grid.hh:1766
Traits::template Codim< codim >::Entity Entity
type of entity
Definition: grid.hh:1710
Traits::template Codim< codim >::EntityPointer EntityPointer
type of entity pointer
Definition: grid.hh:1716
Traits::template Codim< codim >::LocalGeometry LocalGeometry
type of local geometry
Definition: grid.hh:1741
Partition< All_Partition >::LevelIterator LevelIterator
type of leaf iterator
Definition: grid.hh:1775
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