CpGrid.hpp
1 //===========================================================================
2 //
3 // File: CpGrid.hpp
4 //
5 // Created: Fri May 29 20:26:36 2009
6 //
7 // Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8 // B�rd Skaflestad <bard.skaflestad@sintef.no>
9 //
10 // $Date$
11 //
12 // $Revision$
13 //
14 //===========================================================================
15 
16 /*
17  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
18  Copyright 2009, 2010, 2014 Statoil ASA.
19  Copyright 2014, 2015 Dr. Blatt - HPC-Simulartion-Software & Services
20  Copyright 2015 NTNU
21 
22  This file is part of The Open Porous Media project (OPM).
23 
24  OPM is free software: you can redistribute it and/or modify
25  it under the terms of the GNU General Public License as published by
26  the Free Software Foundation, either version 3 of the License, or
27  (at your option) any later version.
28 
29  OPM is distributed in the hope that it will be useful,
30  but WITHOUT ANY WARRANTY; without even the implied warranty of
31  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32  GNU General Public License for more details.
33 
34  You should have received a copy of the GNU General Public License
35  along with OPM. If not, see <http://www.gnu.org/licenses/>.
36 */
37 
38 #ifndef OPM_CPGRID_HEADER
39 #define OPM_CPGRID_HEADER
40 
41 #include <string>
42 #include <map>
43 #include <array>
44 #include <unordered_set>
45 #include <opm/grid/utility/ErrorMacros.hpp>
46 
47 // Warning suppression for Dune includes.
48 #include <opm/grid/utility/platform_dependent/disable_warnings.h>
49 
50 #include <dune/common/version.hh>
51 
52 #if HAVE_MPI
53 #include <dune/common/parallel/variablesizecommunicator.hh>
54 #endif
55 
56 #include <dune/grid/common/capabilities.hh>
57 #include <dune/grid/common/grid.hh>
58 #include <dune/grid/common/gridenums.hh>
59 
60 #include <opm/grid/utility/platform_dependent/reenable_warnings.h>
61 
62 #include "cpgrid/Intersection.hpp"
63 #include "cpgrid/Entity.hpp"
64 #include "cpgrid/Geometry.hpp"
65 #include "cpgrid/Iterators.hpp"
66 #include "cpgrid/Indexsets.hpp"
67 #include "cpgrid/DefaultGeometryPolicy.hpp"
68 #include "common/Volumes.hpp"
70 
71 #include <opm/grid/utility/OpmParserIncludes.hpp>
72 
73 #include <iostream>
74 
75 namespace Dune
76 {
77 
78  class CpGrid;
79 
80  namespace cpgrid
81  {
82  class CpGridData;
83  }
84 
86  //
87  // CpGridTraits
88  //
90 
91  struct CpGridTraits
92  {
94  typedef CpGrid Grid;
95 
104 
107 
110  template <int cd>
111  struct Codim
112  {
115  typedef cpgrid::Geometry<3-cd, 3> Geometry;
116  //typedef Dune::Geometry<3-cd, 3, CpGrid, cpgrid::Geometry> Geometry;
118  typedef cpgrid::Geometry<3-cd, 3> LocalGeometry;
119  //typedef Dune::Geometry<3-cd, 3, CpGrid, cpgrid::Geometry> LocalGeometry;
122 
125 
128 
131 
134 
137  template <PartitionIteratorType pitype>
138  struct Partition
139  {
144  };
145  };
146 
149  template <PartitionIteratorType pitype>
150  struct Partition
151  {
152 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 5)
153  typedef Dune::GridView<DefaultLevelGridViewTraits<CpGrid> > LevelGridView;
156  typedef Dune::GridView<DefaultLeafGridViewTraits<CpGrid> > LeafGridView;
157 #else
158  typedef Dune::GridView<DefaultLevelGridViewTraits<CpGrid, pitype> > LevelGridView;
161  typedef Dune::GridView<DefaultLeafGridViewTraits<CpGrid, pitype> > LeafGridView;
162 #endif
163 
164  };
165 
166 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 5)
167  typedef Dune::GridView<DefaultLevelGridViewTraits<CpGrid> > LevelGridView;
170  typedef Dune::GridView<DefaultLeafGridViewTraits<CpGrid> > LeafGridView;
171 #else
172  typedef Dune::GridView<DefaultLevelGridViewTraits<CpGrid, Dune::All_Partition> > LevelGridView;
175  typedef Dune::GridView<DefaultLeafGridViewTraits<CpGrid, Dune::All_Partition> > LeafGridView;
176 #endif
177 
186 
188 
189  typedef Dune::MPIHelper::MPICommunicator MPICommunicator;
190  typedef Dune::CollectiveCommunication<MPICommunicator> CollectiveCommunication;
191  };
192 
193 
194 
196  //
197  // CpGridFamily
198  //
200 
202  {
203  typedef CpGridTraits Traits;
204  };
205 
207  //
208  // CpGrid
209  //
211 
213  class CpGrid
214  : public GridDefaultImplementation<3, 3, double, CpGridFamily >
215  {
216 
217  // defined in OpmParserIncludes in case opm-parser is not available
218  typedef cpgrid::OpmEclipseStateType OpmEclipseStateType;
219 
220  public:
221 
222  // --- Typedefs ---
223 
224 
227 
228 
229  // --- Methods ---
230 
231 
233  CpGrid();
234 
236 
237  void readSintefLegacyFormat(const std::string& grid_prefix);
241 
242 
246  void writeSintefLegacyFormat(const std::string& grid_prefix) const;
247 
248 
249 #if HAVE_OPM_PARSER
250  void processEclipseFormat(const Opm::EclipseGrid& ecl_grid, bool periodic_extension, bool turn_normals = false, bool clip_z = false,
259  const std::vector<double>& poreVolume = std::vector<double>());
260 #endif
261 
267  void processEclipseFormat(const grdecl& input_data, double z_tolerance, bool remove_ij_boundary, bool turn_normals = false);
268 
270 
276 
277  void createCartesian(const array<int, 3>& dims,
281  const array<double, 3>& cellsize);
282 
286  const std::array<int, 3>& logicalCartesianSize() const
287  {
288  return current_view_data_->logical_cartesian_size_;
289  }
290 
298  const std::vector<int>& globalCell() const
299  {
300  return current_view_data_->global_cell_;
301  }
302 
310  void getIJK(const int c, std::array<int,3>& ijk) const
311  {
312  current_view_data_->getIJK(c, ijk);
313  }
315 
319  bool uniqueBoundaryIds() const
320  {
321  return current_view_data_->uniqueBoundaryIds();
322  }
323 
326  void setUniqueBoundaryIds(bool uids)
327  {
328  current_view_data_->setUniqueBoundaryIds(uids);
329  }
330 
331  // --- Dune interface below ---
332 
334  // \@{
339  std::string name() const
340  {
341  return "CpGrid";
342  }
343 
344 
347  int maxLevel() const
348  {
349  return 0;
350  }
351 
352 
354  template<int codim>
355  typename Traits::template Codim<codim>::LevelIterator lbegin (int level) const
356  {
357  if (level<0 || level>maxLevel())
358  DUNE_THROW(GridError, "levelIndexSet of nonexisting level " << level << " requested!");
359  return cpgrid::Iterator<codim, All_Partition>(*current_view_data_, 0, true);
360  }
361 
362 
364  template<int codim>
365  typename Traits::template Codim<codim>::LevelIterator lend (int level) const
366  {
367  if (level<0 || level>maxLevel())
368  DUNE_THROW(GridError, "levelIndexSet of nonexisting level " << level << " requested!");
369  return cpgrid::Iterator<codim,All_Partition>(*current_view_data_, size(codim), true );
370 
371  }
372 
373 
375  template<int codim, PartitionIteratorType PiType>
376  typename Traits::template Codim<codim>::template Partition<PiType>::LevelIterator lbegin (int level) const
377  {
378  if (level<0 || level>maxLevel())
379  DUNE_THROW(GridError, "levelIndexSet of nonexisting level " << level << " requested!");
380  return cpgrid::Iterator<codim,PiType>(*current_view_data_, 0, true );
381  }
382 
383 
385  template<int codim, PartitionIteratorType PiType>
386  typename Traits::template Codim<codim>::template Partition<PiType>::LevelIterator lend (int level) const
387  {
388  if (level<0 || level>maxLevel())
389  DUNE_THROW(GridError, "levelIndexSet of nonexisting level " << level << " requested!");
390  return cpgrid::Iterator<codim,PiType>(*current_view_data_, size(codim), true);
391  }
392 
393 
395  template<int codim>
396  typename Traits::template Codim<codim>::LeafIterator leafbegin() const
397  {
398  return cpgrid::Iterator<codim, All_Partition>(*current_view_data_, 0, true);
399  }
400 
401 
403  template<int codim>
404  typename Traits::template Codim<codim>::LeafIterator leafend() const
405  {
406  return cpgrid::Iterator<codim, All_Partition>(*current_view_data_, size(codim), true);
407  }
408 
409 
411  template<int codim, PartitionIteratorType PiType>
412  typename Traits::template Codim<codim>::template Partition<PiType>::LeafIterator leafbegin() const
413  {
414  return cpgrid::Iterator<codim, PiType>(*current_view_data_, 0, true);
415  }
416 
417 
419  template<int codim, PartitionIteratorType PiType>
420  typename Traits::template Codim<codim>::template Partition<PiType>::LeafIterator leafend() const
421  {
422  return cpgrid::Iterator<codim, PiType>(*current_view_data_, size(codim), true);
423  }
424 
425 
427  int size (int level, int codim) const
428  {
429  if (level<0 || level>maxLevel())
430  DUNE_THROW(GridError, "levelIndexSet of nonexisting level " << level << " requested!");
431  return size(codim);
432  }
433 
434 
436  int size (int codim) const
437  {
438  return current_view_data_->size(codim);
439  }
440 
441 
443  int size (int level, GeometryType type) const
444  {
445  if (level<0 || level>maxLevel())
446  DUNE_THROW(GridError, "levelIndexSet of nonexisting level " << level << " requested!");
447  return size(type);
448  }
449 
450 
452  int size (GeometryType type) const
453  {
454  return current_view_data_->size(type);
455  }
456 
457 
459  const Traits::GlobalIdSet& globalIdSet() const
460  {
461  return *current_view_data_->global_id_set_;
462  }
463 
464 
466  const Traits::LocalIdSet& localIdSet() const
467  {
468  return *current_view_data_->local_id_set_;
469  }
470 
471 
473  const Traits::LevelIndexSet& levelIndexSet(int level) const
474  {
475  if (level<0 || level>maxLevel())
476  DUNE_THROW(GridError, "levelIndexSet of nonexisting level " << level << " requested!");
477  return *current_view_data_->index_set_;
478  }
479 
480 
482  const Traits::LeafIndexSet& leafIndexSet() const
483  {
484  return *current_view_data_->index_set_;
485  }
486 
487 
489  void globalRefine (int)
490  {
491  std::cout << "Warning: Global refinement not implemented, yet." << std::endl;
492  }
493 
494  const std::vector< Dune :: GeometryType >& geomTypes( const int codim ) const
495  {
496  return leafIndexSet().geomTypes( codim );
497  }
498 
500  template <int codim>
502  {
503  return cpgrid::Entity<codim>( *seed );
504  }
505 
506  /* No refinement implemented. GridDefaultImplementation's methods will be used.
507 
517 
518  bool mark(int refCount, const typename Traits::template Codim<0>::EntityPointer & e)
519  {
520  return hostgrid_->mark(refCount, getHostEntity<0>(*e));
521  }
522 
526 
527  int getMark(const typename Traits::template Codim<0>::EntityPointer & e) const
528  {
529  return hostgrid_->getMark(getHostEntity<0>(*e));
530  }
531 
533  bool preAdapt() {
534  return hostgrid_->preAdapt();
535  }
536 
537 
539  bool adapt()
540  {
541  return hostgrid_->adapt();
542  }
543 
545  void postAdapt() {
546  return hostgrid_->postAdapt();
547  }
548 
549  end of refinement section */
550 
552  unsigned int overlapSize(int) const {
553  return 1;
554  }
555 
556 
558  unsigned int ghostSize(int) const {
559  return 0;
560  }
561 
562 
564  unsigned int overlapSize(int, int) const {
565  return 1;
566  }
567 
568 
570  unsigned int ghostSize(int, int) const {
571  return 0;
572  }
573 
575  unsigned int numBoundarySegments() const
576  {
577  if( uniqueBoundaryIds() )
578  {
579  return current_view_data_->unique_boundary_ids_.size();
580  }
581  else
582  {
583  unsigned int numBndSegs = 0;
584  const int num_faces = numFaces();
585  for (int i = 0; i < num_faces; ++i) {
586  cpgrid::EntityRep<1> face(i, true);
587  if (current_view_data_->face_to_cell_[face].size() == 1) {
588  ++numBndSegs;
589  }
590  }
591  return numBndSegs;
592  }
593  }
594 
595 
596  // loadbalance is not part of the grid interface therefore we skip it.
597 
601  bool loadBalance(int overlapLayers=1)
602  {
603  using std::get;
604  return get<0>(scatterGrid(nullptr, nullptr, overlapLayers ));
605  }
606 
607  // loadbalance is not part of the grid interface therefore we skip it.
608 
619  std::pair<bool, std::unordered_set<std::string> >
620  loadBalance(const OpmEclipseStateType* ecl,
621  const double* transmissibilities = nullptr,
622  int overlapLayers=1)
623  {
624  return scatterGrid(ecl, transmissibilities, overlapLayers);
625  }
626 
640  template<class DataHandle>
641  std::pair<bool, std::unordered_set<std::string> >
642  loadBalance(DataHandle& data,
643  const OpmEclipseStateType* ecl,
644  const double* transmissibilities = nullptr,
645  int overlapLayers=1)
646  {
647  auto ret = loadBalance(ecl, transmissibilities, overlapLayers);
648  scatterData(data);
649  return ret;
650  }
651 
658  template<class DataHandle>
659  bool loadBalance(DataHandle& data,
660  int overlapLayers=1)
661  {
662  bool ret = loadBalance(overlapLayers);
663  scatterData(data);
664  return ret;
665  }
666 
674  template<class DataHandle>
675  void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int /*level*/) const
676  {
677  communicate(data, iftype, dir);
678  }
679 
687  template<class DataHandle>
688  void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir) const
689  {
690  current_view_data_->communicate(data, iftype, dir);
691  }
692 
694  const CollectiveCommunication& comm () const
695  {
696  return current_view_data_->ccobj_;
697  }
699 
700  // ------------ End of Dune interface, start of simplified interface --------------
701 
707 
708  // enum { dimension = 3 }; // already defined
709 
710  typedef Dune::FieldVector<double, 3> Vector;
711 
712 
713  const std::vector<double>& zcornData() const {
714  return data_->zcornData();
715  }
716 
717 
718  // Topology
720  int numCells() const
721  {
722  return current_view_data_->cell_to_face_.size();
723  }
725  int numFaces() const
726  {
727  return current_view_data_->face_to_cell_.size();
728  }
730  int numVertices() const
731  {
732  return current_view_data_->geomVector<3>().size();
733  }
740  int numCellFaces(int cell) const
741  {
742  return current_view_data_->cell_to_face_[cpgrid::EntityRep<0>(cell, true)].size();
743  }
748  int cellFace(int cell, int local_index) const
749  {
750  return current_view_data_->cell_to_face_[cpgrid::EntityRep<0>(cell, true)][local_index].index();
751  }
752 
756  {
757  return current_view_data_->cell_to_face_[cpgrid::EntityRep<0>(cell, true)];
758  }
769  int faceCell(int face, int local_index) const
770  {
771  // In the parallel case we store non-existent cells for faces along
772  // the front region. Theses marked with index std::numeric_limits<int>::max(),
773  // orientation might be arbitrary, though.
775  = current_view_data_->face_to_cell_[cpgrid::EntityRep<1>(face, true)];
776  bool a = (local_index == 0);
777  bool b = r[0].orientation();
778  bool use_first = a ? b : !b;
779  // The number of valid cells.
780  int r_size = r.size();
781  // In the case of only one valid cell, this is the index of it.
782  int index = 0;
783  if(r[0].index()==std::numeric_limits<int>::max()){
784  assert(r_size==2);
785  --r_size;
786  index=1;
787  }
788  if(r.size()>1 && r[1].index()==std::numeric_limits<int>::max())
789  {
790  assert(r_size==2);
791  --r_size;
792  }
793  if (r_size == 2) {
794  return use_first ? r[0].index() : r[1].index();
795  } else {
796  return use_first ? r[index].index() : -1;
797  }
798  }
805  int numCellFaces() const
806  {
807  return current_view_data_->cell_to_face_.dataSize();
808  }
809  int numFaceVertices(int face) const
810  {
811  return current_view_data_->face_to_point_[face].size();
812  }
817  int faceVertex(int face, int local_index) const
818  {
819  return current_view_data_->face_to_point_[face][local_index];
820  }
823  double cellCenterDepth(int cell_index) const
824  {
825  // Here cell center depth is computed as a raw average of cell corner depths.
826  // This generally gives slightly different results than using the cell centroid.
827  double zz = 0.0;
828  const int nv = current_view_data_->cell_to_point_[cell_index].size();
829  const int nd = 3;
830  for (int i=0; i<nv; ++i) {
831  zz += vertexPosition(current_view_data_->cell_to_point_[cell_index][i])[nd-1];
832  }
833  return zz/nv;
834  }
835 
836  const Vector faceCenterEcl(int cell_index, int faceTag) const
837  {
838  // This method is an alternative to the method faceCentroid(...).
839  // The face center is computed as a raw average of cell corners.
840  // For faulted cells this gives different results then average of face nodes
841  // that seems to agree more with eclipse.
842  // This assumes the cell nodes are ordered
843  // 6---7
844  // | T |
845  // 4---5
846  // 2---3
847  // | B |
848  // 0---1
849 
850  // this follows the DUNE reference cube
851  static const int faceVxMap[ 6 ][ 4 ] = { {0, 2, 4, 6}, // face 0
852  {1, 3, 5, 7}, // face 1
853  {0, 1, 4, 5}, // face 2
854  {2, 3, 6, 7}, // face 3
855  {0, 1, 2, 3}, // face 4
856  {4, 5, 6, 7} // face 5
857  };
858 
859 
860  assert (current_view_data_->cell_to_point_[cell_index].size() == 8);
861  Vector center(0.0);
862  for( int i=0; i<4; ++i )
863  {
864  center += vertexPosition(current_view_data_->cell_to_point_[cell_index][ faceVxMap[ faceTag ][ i ] ]);
865  }
866 
867  for (int i=0; i<3; ++i) {
868  center[i] /= 4;
869  }
870  return center;
871 
872  }
873 
874  const Vector faceAreaNormalEcl(int face) const
875  {
876  // same implementation as ResInsight
877  const int nd = Vector::dimension;
878  const int nv = numFaceVertices(face);
879  switch (nv)
880  {
881  case 0:
882  case 1:
883  case 2:
884  {
885  return Vector(0.0);
886  }
887  break;
888  case 3:
889  {
890  Vector a = vertexPosition(current_view_data_->face_to_point_[face][0]) - vertexPosition(current_view_data_->face_to_point_[face][2]);
891  Vector b = vertexPosition(current_view_data_->face_to_point_[face][1]) - vertexPosition(current_view_data_->face_to_point_[face][2]);
892  Vector areaNormal = cross(a,b);
893  for (int i=0; i<nd; ++i) {
894  areaNormal[i] /= 2;
895  }
896  return areaNormal;
897  }
898  break;
899  case 4:
900  {
901  Vector a = vertexPosition(current_view_data_->face_to_point_[face][0]) - vertexPosition(current_view_data_->face_to_point_[face][2]);
902  Vector b = vertexPosition(current_view_data_->face_to_point_[face][1]) - vertexPosition(current_view_data_->face_to_point_[face][3]);
903  Vector areaNormal = cross(a,b);
904  areaNormal *= 0.5;
905  return areaNormal;
906  }
907  break;
908  default:
909  {
910  int h = (nv - 1)/2;
911  int k = (nv % 2) ? 0 : nv - 1;
912 
913  Vector areaNormal(0.0);
914  // First quads
915  for (int i = 1; i < h; ++i)
916  {
917  Vector a = vertexPosition(current_view_data_->face_to_point_[face][2*i]) - vertexPosition(current_view_data_->face_to_point_[face][0]);
918  Vector b = vertexPosition(current_view_data_->face_to_point_[face][2*i+1]) - vertexPosition(current_view_data_->face_to_point_[face][2*i-1]);
919  areaNormal += cross(a,b);
920  }
921 
922  // Last triangle or quad
923  Vector a = vertexPosition(current_view_data_->face_to_point_[face][2*h]) - vertexPosition(current_view_data_->face_to_point_[face][0]);
924  Vector b = vertexPosition(current_view_data_->face_to_point_[face][k]) - vertexPosition(current_view_data_->face_to_point_[face][2*h-1]);
925  areaNormal += cross(a,b);
926 
927  areaNormal *= 0.5;
928 
929  return areaNormal;
930  }
931 
932  }
933  }
934 
935  // Geometry
939  const Vector& vertexPosition(int vertex) const
940  {
941  return current_view_data_->geomVector<3>()[cpgrid::EntityRep<3>(vertex, true)].center();
942  }
945  double faceArea(int face) const
946  {
947  return current_view_data_->geomVector<1>()[cpgrid::EntityRep<1>(face, true)].volume();
948  }
951  const Vector& faceCentroid(int face) const
952  {
953  return current_view_data_->geomVector<1>()[cpgrid::EntityRep<1>(face, true)].center();
954  }
958  const Vector& faceNormal(int face) const
959  {
960  return current_view_data_->face_normals_.get(face);
961  }
964  double cellVolume(int cell) const
965  {
966  return current_view_data_->geomVector<0>()[cpgrid::EntityRep<0>(cell, true)].volume();
967  }
970  const Vector& cellCentroid(int cell) const
971  {
972  return current_view_data_->geomVector<0>()[cpgrid::EntityRep<0>(cell, true)].center();
973  }
974 
977  template<int codim>
979  : public RandomAccessIteratorFacade<CentroidIterator<codim>,
980  FieldVector<double, 3>,
981  const FieldVector<double, 3>&, int>
982  {
983  public:
985  typedef typename std::vector<cpgrid::Geometry<3-codim, 3> >::const_iterator
990  : iter_(iter)
991  {}
992 
993  const FieldVector<double, 3>& dereference() const
994  {
995  return iter_->center();
996  }
997  void increment()
998  {
999  ++iter_;
1000  }
1001  const FieldVector<double, 3>& elementAt(int n)
1002  {
1003  return iter_[n]->center();
1004  }
1005  void advance(int n)
1006  {
1007  iter_+=n;
1008  }
1009  void decrement()
1010  {
1011  --iter_;
1012  }
1013  int distanceTo(const CentroidIterator& o)
1014  {
1015  return o-iter_;
1016  }
1017  bool equals(const CentroidIterator& o) const
1018  {
1019  return o==iter_;
1020  }
1021  private:
1023  GeometryIterator iter_;
1024  };
1025 
1028  {
1029  return CentroidIterator<0>(current_view_data_->geomVector<0>().begin());
1030  }
1031 
1034  {
1035  return CentroidIterator<1>(current_view_data_->geomVector<1>().begin());
1036  }
1037 
1038  // Extra
1039  int boundaryId(int face) const
1040  {
1041  // Note that this relies on the following implementation detail:
1042  // The grid is always construct such that the faces where
1043  // orientation() returns true are oriented along the positive IJK
1044  // direction. Oriented means that the first cell attached to face
1045  // has the lower index.
1046  int ret = 0;
1047  cpgrid::EntityRep<1> f(face, true);
1048  if (current_view_data_->face_to_cell_[f].size() == 1) {
1049  if (current_view_data_->uniqueBoundaryIds()) {
1050  // Use the unique boundary ids.
1051  ret = current_view_data_->unique_boundary_ids_[f];
1052  } else {
1053  // Use the face tag based ids, i.e. 1-6 for i-, i+, j-, j+, k-, k+.
1054  const bool normal_is_in =
1055  !(current_view_data_->face_to_cell_[f][0].orientation());
1056  enum face_tag tag = current_view_data_->face_tag_[f];
1057  switch (tag) {
1058  case LEFT:
1059  // LEFT : RIGHT
1060  ret = normal_is_in ? 1 : 2; // min(I) : max(I)
1061  break;
1062  case BACK:
1063  // BACK : FRONT
1064  ret = normal_is_in ? 3 : 4; // min(J) : max(J)
1065  break;
1066  case TOP:
1067  // Note: TOP at min(K) as 'z' measures *depth*.
1068  // TOP : BOTTOM
1069  ret = normal_is_in ? 5 : 6; // min(K) : max(K)
1070  break;
1071  }
1072  }
1073  }
1074  return ret;
1075  }
1076 
1083  template<class Cell2FacesRowIterator>
1084  int
1085  faceTag(const Cell2FacesRowIterator& cell_face) const
1086  {
1087  // Note that this relies on the following implementation detail:
1088  // The grid is always construct such that the interior faces constructed
1089  // with orientation set to true are
1090  // oriented along the positive IJK direction. Oriented means that
1091  // the first cell attached to face has the lower index.
1092  // For faces along the boundary (only one cell, always attached at index 0)
1093  // the orientation has to be determined by the orientation of the cell.
1094  // If it is true then in UnstructuredGrid it would be stored at index 0,
1095  // otherwise at index 1.
1096  const int cell = cell_face.getCellIndex();
1097  const int face = *cell_face;
1098  assert (0 <= cell); assert (cell < numCells());
1099  assert (0 <= face); assert (face < numFaces());
1100 
1102 
1103  const cpgrid::EntityRep<1> f(face, true);
1104  const F2C& f2c = current_view_data_->face_to_cell_[f];
1105  const face_tag tag = current_view_data_->face_tag_[f];
1106 
1107  assert ((f2c.size() == 1) || (f2c.size() == 2));
1108 
1109  int inside_cell = 0;
1110 
1111  if ( f2c.size() == 2 ) // Two cells => interior
1112  {
1113  if ( f2c[1].index() == cell )
1114  {
1115  inside_cell = 1;
1116  }
1117  }
1118  const bool normal_is_in = ! f2c[inside_cell].orientation();
1119 
1120  switch (tag) {
1121  case LEFT:
1122  // LEFT : RIGHT
1123  return normal_is_in ? 0 : 1; // min(I) : max(I)
1124 
1125  case BACK:
1126  // BACK : FRONT
1127  return normal_is_in ? 2 : 3; // min(J) : max(J)
1128 
1129  case TOP:
1130  // Note: TOP at min(K) as 'z' measures *depth*.
1131  // TOP : BOTTOM
1132  return normal_is_in ? 4 : 5; // min(K) : max(K)
1133  default:
1134  OPM_THROW(std::logic_error, "Unhandled face tag. This should never happen!");
1135  }
1136  }
1137 
1139 
1140  // ------------ End of simplified interface --------------
1141 
1142  //------------- methods not in the DUNE grid interface.
1143 
1148 
1149  template<class DataHandle>
1159  void scatterData(DataHandle& handle)
1160  {
1161 #if HAVE_MPI
1162  if(!distributed_data_)
1163  OPM_THROW(std::runtime_error, "Moving Data only allowed with a load balanced grid!");
1164  distributed_data_->scatterData(handle, data_.get(), distributed_data_.get());
1165 #else
1166  // Suppress warnings for unused argument.
1167  (void) handle;
1168 #endif
1169  }
1170 
1177  template<class DataHandle>
1178  void gatherData(DataHandle& handle)
1179  {
1180 #if HAVE_MPI
1181  if(!distributed_data_)
1182  OPM_THROW(std::runtime_error, "Moving Data only allowed with a load balance grid!");
1183  distributed_data_->gatherData(handle, data_.get(), distributed_data_.get());
1184 #else
1185  // Suppress warnings for unused argument.
1186  (void) handle;
1187 #endif
1188  }
1189 #if HAVE_MPI
1190  typedef VariableSizeCommunicator<>::InterfaceMap InterfaceMap;
1192 #else
1193  // bogus definition for the non parallel type. VariableSizeCommunicator not
1194  // availabe
1195 
1197  typedef std::map<int, std::list<int> > InterfaceMap;
1198 #endif
1199 
1229  {
1230  return *cell_scatter_gather_interfaces_;
1231  }
1232 
1235  {
1236  current_view_data_=data_.get();
1237  }
1238 
1241  {
1242  current_view_data_=distributed_data_.get();
1243  }
1245 
1246 #if HAVE_MPI
1247  typedef cpgrid::CpGridData::ParallelIndexSet ParallelIndexSet;
1250  typedef cpgrid::CpGridData::RemoteIndices RemoteIndices;
1251 
1252  ParallelIndexSet& getCellIndexSet()
1253  {
1254  return current_view_data_->cell_indexset_;
1255  }
1256 
1257  RemoteIndices& getCellRemoteIndices()
1258  {
1259  return current_view_data_->cell_remote_indices_;
1260  }
1261 
1262  const ParallelIndexSet& getCellIndexSet() const
1263  {
1264  return current_view_data_->cell_indexset_;
1265  }
1266 
1267  const RemoteIndices& getCellRemoteIndices() const
1268  {
1269  return current_view_data_->cell_remote_indices_;
1270  }
1271 
1272 #endif
1273 
1274  private:
1283  std::pair<bool, std::unordered_set<std::string> >
1284  scatterGrid(const OpmEclipseStateType* ecl, const double* transmissibilities,
1285  int overlapLayers);
1286 
1291  std::shared_ptr<cpgrid::CpGridData> data_;
1293  cpgrid::CpGridData* current_view_data_;
1295  std::shared_ptr<cpgrid::CpGridData> distributed_data_;
1301  std::shared_ptr<InterfaceMap> cell_scatter_gather_interfaces_;
1302  }; // end Class CpGrid
1303 
1304 
1305 
1306  namespace Capabilities
1307  {
1309  template <>
1310  struct hasEntity<CpGrid, 0>
1311  {
1312  static const bool v = true;
1313  };
1314 
1316  template <>
1317  struct hasEntity<CpGrid, 3>
1318  {
1319  static const bool v = true;
1320  };
1321 
1322 #if ! DUNE_VERSION_NEWER(DUNE_GRID, 2, 5)
1323  template <>
1325  struct isParallel<CpGrid>
1326  {
1327  static const bool v = true;
1328  };
1329 #endif
1330 
1331  template<>
1332  struct canCommunicate<CpGrid,0>
1333  {
1334  static const bool v = true;
1335  };
1336 
1337  template<>
1338  struct canCommunicate<CpGrid,3>
1339  {
1340  static const bool v = true;
1341  };
1342 
1344  template <>
1345  struct hasBackupRestoreFacilities<CpGrid>
1346  {
1347  static const bool v = false;
1348  };
1349 
1350  }
1351 
1352 } // namespace Dune
1353 
1354 #include <dune/grid/cpgrid/PersistentContainer.hpp>
1355 #include <dune/grid/cpgrid/CartesianIndexMapper.hpp>
1356 #endif // OPM_CPGRID_HEADER
const InterfaceMap & cellScatterGatherInterface() const
Get an interface for gathering/scattering data with communication.
Definition: CpGrid.hpp:1228
Traits::template Codim< codim >::LevelIterator lbegin(int level) const
Iterator to first entity of given codim on level.
Definition: CpGrid.hpp:355
const Traits::GlobalIdSet & globalIdSet() const
Access to the GlobalIdSet.
Definition: CpGrid.hpp:459
cpgrid::IdSet LocalIdSet
The type of the local id set.
Definition: CpGrid.hpp:185
Definition: CpGrid.hpp:201
Definition: CpGridData.hpp:93
Dune::MPIHelper::MPICommunicator MPICommunicator
The type of the collective communication.
Definition: CpGrid.hpp:189
void processEclipseFormat(const grdecl &input_data, double z_tolerance, bool remove_ij_boundary, bool turn_normals=false)
Read the Eclipse grid format (&#39;grdecl&#39;).
Definition: CpGrid.cpp:276
int cellFace(int cell, int local_index) const
Get a specific face of a cell.
Definition: CpGrid.hpp:748
std::vector< cpgrid::Geometry< 3-codim, 3 > >::const_iterator GeometryIterator
The type of the iterator over the codim geometries.
Definition: CpGrid.hpp:986
Traits associated with a specific grid partition type.
Definition: CpGrid.hpp:138
const Vector & cellCentroid(int cell) const
Get the coordinates of the center of a cell.
Definition: CpGrid.hpp:970
void switchToGlobalView()
Switch to the global view.
Definition: CpGrid.hpp:1234
unsigned int ghostSize(int, int) const
Size of the ghost cell layer on a given level.
Definition: CpGrid.hpp:570
double faceArea(int face) const
Get the area of a face.
Definition: CpGrid.hpp:945
unsigned int numBoundarySegments() const
returns the number of boundary segments within the macro grid
Definition: CpGrid.hpp:575
Dune::GridView< DefaultLevelGridViewTraits< CpGrid, Dune::All_Partition > > LevelGridView
The type of the level grid view associated with this partition type.
Definition: CpGrid.hpp:173
cpgrid::HierarchicIterator HierarchicIterator
The type of the hierarchic iterator.
Definition: CpGrid.hpp:106
Raw corner-point specification of a particular geological model.
Definition: preprocess.h:56
cpgrid::IndexSet LeafIndexSet
The type of the leaf index set.
Definition: CpGrid.hpp:181
cpgrid::GlobalIdSet GlobalIdSet
The type of the global id set.
Definition: CpGrid.hpp:183
Definition: CpGrid.hpp:91
Dune::GridView< DefaultLeafGridViewTraits< CpGrid, Dune::All_Partition > > LeafGridView
The type of the leaf grid view associated with this partition type.
Definition: CpGrid.hpp:175
int faceTag(const Cell2FacesRowIterator &cell_face) const
Get the cartesian tag associated with a face tag.
Definition: CpGrid.hpp:1085
const Vector & faceCentroid(int face) const
Get the coordinates of the center of a face.
Definition: CpGrid.hpp:951
Holds the implementation of the CpGrid as a pimple.
Definition: OpmParserIncludes.hpp:42
std::pair< bool, std::unordered_set< std::string > > loadBalance(const OpmEclipseStateType *ecl, const double *transmissibilities=nullptr, int overlapLayers=1)
Distributes this grid over the available nodes in a distributed machine.
Definition: CpGrid.hpp:620
void createCartesian(const array< int, 3 > &dims, const array< double, 3 > &cellsize)
Create a cartesian grid.
Definition: CpGrid.cpp:212
int size(int codim) const
number of leaf entities per codim in this process
Definition: CpGrid.hpp:436
void readSintefLegacyFormat(const std::string &grid_prefix)
Read the Sintef legacy grid format (&#39;topogeom&#39;).
Definition: CpGrid.cpp:254
Definition: Indexsets.hpp:192
void writeSintefLegacyFormat(const std::string &grid_prefix) const
Write the Sintef legacy grid format (&#39;topogeom&#39;).
Definition: CpGrid.cpp:258
Definition: Intersection.hpp:67
cpgrid::Intersection LeafIntersection
The type of the intersection at the leafs of the grid.
Definition: CpGrid.hpp:97
int size(int level, int codim) const
Number of grid entities per level and codim.
Definition: CpGrid.hpp:427
void setUniqueBoundaryIds(bool uids)
Set whether we want to have unique boundary ids.
Definition: CpGrid.hpp:326
int size() const
Returns the number of rows in the table.
Definition: SparseTable.hpp:124
const cpgrid::OrientedEntityTable< 0, 1 >::row_type cellFaceRow(int cell) const
Get a list of indices identifying all faces of a cell.
Definition: CpGrid.hpp:755
CentroidIterator< 1 > beginFaceCentroids() const
Get an iterator over the face centroids positioned at the first one.
Definition: CpGrid.hpp:1033
bool uniqueBoundaryIds() const
Is the grid currently using unique boundary ids?
Definition: CpGrid.hpp:319
[ provides Dune::Grid ]
Definition: CpGrid.hpp:213
Traits::template Codim< codim >::template Partition< PiType >::LevelIterator lend(int level) const
one past the end on this level
Definition: CpGrid.hpp:386
Traits associated with a specific grid partition type.
Definition: CpGrid.hpp:150
cpgrid::IntersectionIterator LevelIntersectionIterator
The type of the intersection iterator at the levels of the grid.
Definition: CpGrid.hpp:103
This class encapsulates geometry for both vertices, intersections and cells.
Definition: CpGridData.hpp:92
cpgrid::Entity< codim > entity(const cpgrid::EntityPointer< codim > &seed) const
given an EntitySeed (or EntityPointer) return an entity object
Definition: CpGrid.hpp:501
const std::vector< int > & globalCell() const
Retrieve mapping from internal ("compressed") active grid cells to external ("uncompressed") cells...
Definition: CpGrid.hpp:298
Iterator intended to be used as LeafIterator and LevelIterator (no difference due to no adaptivity) f...
Definition: Entity.hpp:52
std::map< int, std::list< int > > InterfaceMap
The type of the map describing communication interfaces.
Definition: CpGrid.hpp:1197
void setUniqueBoundaryIds(bool uids)
Set whether we want to have unique boundary ids.
Definition: CpGridData.hpp:215
cpgrid::Entity< cd > Entity
The type of the entity.
Definition: CpGrid.hpp:121
void communicate(DataHandle &data, InterfaceType iftype, CommunicationDirection dir)
communicate objects for all codims on a given level
Definition: CpGridData.hpp:512
Traits::template Codim< codim >::LevelIterator lend(int level) const
one past the end on this level
Definition: CpGrid.hpp:365
CentroidIterator< 0 > beginCellCentroids() const
Get an iterator over the cell centroids positioned at the first one.
Definition: CpGrid.hpp:1027
bool loadBalance(DataHandle &data, int overlapLayers=1)
Distributes this grid and data over the available nodes in a distributed machine. ...
Definition: CpGrid.hpp:659
int numFaces() const
Get the number of faces.
Definition: CpGrid.hpp:725
const Vector & faceNormal(int face) const
Get the unit normal of a face.
Definition: CpGrid.hpp:958
int size(int level, GeometryType type) const
number of entities per level and geometry type in this process
Definition: CpGrid.hpp:443
unsigned int overlapSize(int) const
Size of the overlap on the leaf level.
Definition: CpGrid.hpp:552
int numVertices() const
Get The number of vertices.
Definition: CpGrid.hpp:730
double cellCenterDepth(int cell_index) const
Get vertical position of cell center ("zcorn" average).
Definition: CpGrid.hpp:823
const Traits::LeafIndexSet & leafIndexSet() const
Access to the LeafIndexSet.
Definition: CpGrid.hpp:482
CpGrid Grid
The type that implements the grid.
Definition: CpGrid.hpp:94
Traits::template Codim< codim >::template Partition< PiType >::LevelIterator lbegin(int level) const
Iterator to first entity of given codim on level.
Definition: CpGrid.hpp:376
Traits::template Codim< codim >::LeafIterator leafbegin() const
Iterator to first leaf entity of given codim.
Definition: CpGrid.hpp:396
cpgrid::Iterator< cd, pitype > LevelIterator
The type of the iterator over the level entities of this codim on this partition. ...
Definition: CpGrid.hpp:141
Traits::template Codim< codim >::template Partition< PiType >::LeafIterator leafend() const
one past the end of the sequence of leaf entities
Definition: CpGrid.hpp:420
bool uniqueBoundaryIds() const
Is the grid currently using unique boundary ids?
Definition: CpGridData.hpp:208
void switchToDistributedView()
Switch to the distributed view.
Definition: CpGrid.hpp:1240
cpgrid::IntersectionIterator LeafIntersectionIterator
The type of the intersection iterator at the leafs of the grid.
Definition: CpGrid.hpp:101
CpGrid()
Default constructor.
Definition: CpGrid.cpp:58
double cellVolume(int cell) const
Get the volume of the cell.
Definition: CpGrid.hpp:964
unsigned int overlapSize(int, int) const
Size of the overlap on a given level.
Definition: CpGrid.hpp:564
Traits associated with a specific codim.
Definition: CpGrid.hpp:111
void communicate(DataHandle &data, InterfaceType iftype, CommunicationDirection dir, int) const
The new communication interface.
Definition: CpGrid.hpp:675
void communicate(DataHandle &data, InterfaceType iftype, CommunicationDirection dir) const
The new communication interface.
Definition: CpGrid.hpp:688
unsigned int ghostSize(int) const
Size of the ghost cell layer on the leaf level.
Definition: CpGrid.hpp:558
int maxLevel() const
Return maximum level defined in this grid.
Definition: CpGrid.hpp:347
int size(int codim) const
number of leaf entities per codim in this process
Definition: CpGridData.cpp:132
cpgrid::Iterator< cd, All_Partition > LeafIterator
The type of the iterator over all leaf entities of this codim.
Definition: CpGrid.hpp:127
void getIJK(int c, std::array< int, 3 > &ijk) const
Extract Cartesian index triplet (i,j,k) of an active cell.
Definition: CpGridData.hpp:194
Connection topologically parallel to I-K plane.
Definition: preprocess.h:69
An iterator over the centroids of the geometry of the entities.
Definition: CpGrid.hpp:978
std::pair< bool, std::unordered_set< std::string > > loadBalance(DataHandle &data, const OpmEclipseStateType *ecl, const double *transmissibilities=nullptr, int overlapLayers=1)
Distributes this grid and data over the available nodes in a distributed machine. ...
Definition: CpGrid.hpp:642
Only needs to provide interface for doing nothing.
Definition: Iterators.hpp:88
cpgrid::EntityPointer< cd > EntitySeed
The type of the entity pointer for entities of this codim.
Definition: CpGrid.hpp:133
Definition: Indexsets.hpp:250
Connection topologically parallel to I-J plane.
Definition: preprocess.h:70
cpgrid::Iterator< cd, pitype > LeafIterator
The type of the iterator over the leaf entities of this codim on this partition.
Definition: CpGrid.hpp:143
int faceVertex(int face, int local_index) const
Get the index identifying a vertex of a face.
Definition: CpGrid.hpp:817
cpgrid::Geometry< 3-cd, 3 > LocalGeometry
The type of the local geometry associated with the entity.
Definition: CpGrid.hpp:118
face_tag
Connection taxonomy.
Definition: preprocess.h:67
cpgrid::EntityPointer< cd > EntityPointer
The type of the entity pointer for entities of this codim.
Definition: CpGrid.hpp:130
CpGridFamily GridFamily
Family typedef, why is this not defined by Grid<>?
Definition: CpGrid.hpp:226
Dune::GridView< DefaultLeafGridViewTraits< CpGrid, pitype > > LeafGridView
The type of the leaf grid view associated with this partition type.
Definition: CpGrid.hpp:161
cpgrid::IndexSet LevelIndexSet
The type of the level index set.
Definition: CpGrid.hpp:179
const CollectiveCommunication & comm() const
Get the collective communication object.
Definition: CpGrid.hpp:694
Traits::template Codim< codim >::LeafIterator leafend() const
one past the end of the sequence of leaf entities
Definition: CpGrid.hpp:404
int dataSize() const
Returns the number of data elements.
Definition: SparseTable.hpp:144
CentroidIterator(GeometryIterator iter)
Constructs a new iterator from an iterator over the geometries.
Definition: CpGrid.hpp:989
int numCellFaces(int cell) const
Get the number of faces of a cell.
Definition: CpGrid.hpp:740
void globalRefine(int)
global refinement
Definition: CpGrid.hpp:489
cpgrid::Geometry< 3-cd, 3 > Geometry
The type of the geometry associated with the entity.
Definition: CpGrid.hpp:115
Low-level corner-point processing routines and supporting data structures.
const Traits::LevelIndexSet & levelIndexSet(int level) const
Access to the LevelIndexSets.
Definition: CpGrid.hpp:473
Represents an entity of a given codim, with positive or negative orientation.
Definition: CpGridData.hpp:94
const std::array< int, 3 > & logicalCartesianSize() const
The logical cartesian size of the global grid.
Definition: CpGrid.hpp:286
Definition: Indexsets.hpp:52
cpgrid::Intersection LevelIntersection
The type of the intersection at the levels of the grid.
Definition: CpGrid.hpp:99
int numCellFaces() const
Get the sum of all faces attached to all cells.
Definition: CpGrid.hpp:805
Definition: Intersection.hpp:286
void scatterData(DataHandle &handle)
Moves data from the global (all data on process) view to the distributed view.
Definition: CpGrid.hpp:1159
int faceCell(int face, int local_index) const
Get the index identifying a cell attached to a face.
Definition: CpGrid.hpp:769
std::string name() const
Get the grid name.
Definition: CpGrid.hpp:339
int numCells() const
Get the number of cells.
Definition: CpGrid.hpp:720
A class used as a row type for OrientedEntityTable.
Definition: OrientedEntityTable.hpp:55
const Vector & vertexPosition(int vertex) const
Get the Position of a vertex.
Definition: CpGrid.hpp:939
int size(GeometryType type) const
number of leaf entities per geometry type in this process
Definition: CpGrid.hpp:452
void gatherData(DataHandle &handle)
Moves data from the distributed view to the global (all data on process) view.
Definition: CpGrid.hpp:1178
void getIJK(const int c, std::array< int, 3 > &ijk) const
Extract Cartesian index triplet (i,j,k) of an active cell.
Definition: CpGrid.hpp:310
Connection topologically parallel to J-K plane.
Definition: preprocess.h:68
FieldVector< T, 3 > cross(const FieldVector< T, 3 > &a, const FieldVector< T, 3 > &b)
Definition: Volumes.hpp:58
int size() const
Returns the number of rows in the table.
Definition: SparseTable.hpp:124
Dune::GridView< DefaultLevelGridViewTraits< CpGrid, pitype > > LevelGridView
The type of the level grid view associated with this partition type.
Definition: CpGrid.hpp:159
cpgrid::Iterator< cd, All_Partition > LevelIterator
The type of the iterator over all level entities of this codim.
Definition: CpGrid.hpp:124
Traits::template Codim< codim >::template Partition< PiType >::LeafIterator leafbegin() const
Iterator to first leaf entity of given codim.
Definition: CpGrid.hpp:412
bool loadBalance(int overlapLayers=1)
Distributes this grid over the available nodes in a distributed machine.
Definition: CpGrid.hpp:601
const Traits::LocalIdSet & localIdSet() const
Access to the LocalIdSet.
Definition: CpGrid.hpp:466