GridInterfaceEuler.hpp
1 //===========================================================================
2 //
3 // File: GridInterfaceEuler.hpp
4 //
5 // Created: Mon Jun 15 12:53:38 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 Statoil ASA.
19 
20  This file is part of The Open Reservoir Simulator Project (OpenRS).
21 
22  OpenRS is free software: you can redistribute it and/or modify
23  it under the terms of the GNU General Public License as published by
24  the Free Software Foundation, either version 3 of the License, or
25  (at your option) any later version.
26 
27  OpenRS is distributed in the hope that it will be useful,
28  but WITHOUT ANY WARRANTY; without even the implied warranty of
29  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30  GNU General Public License for more details.
31 
32  You should have received a copy of the GNU General Public License
33  along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
34 */
35 
36 #ifndef OPENRS_GRIDINTERFACEEULER_HEADER
37 #define OPENRS_GRIDINTERFACEEULER_HEADER
38 
39 #include <opm/core/utility/SparseTable.hpp>
40 #include <opm/core/utility/StopWatch.hpp>
41 #include <dune/grid/CpGrid.hpp> // How to avoid this? Needed for the explicit mapper specialization below.
42 
43 #include <opm/common/utility/platform_dependent/disable_warnings.h>
44 
45 #include <dune/common/fvector.hh>
46 #include <dune/grid/common/mcmgmapper.hh>
47 
48 #include <boost/iterator/iterator_facade.hpp>
49 #include <boost/scoped_ptr.hpp>
50 
51 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
52 
53 #include <climits>
54 #include <iostream>
55 
56 
57 
58 namespace Opm
59 {
60 
62  template<int dim>
63  struct AllCellsLayout {
64  bool contains(Dune::GeometryType gt) { return gt.dim() == dim; }
65  };
66 
68  template <class DuneGrid>
69  struct GICellMapper
70  {
71  typedef Dune::LeafMultipleCodimMultipleGeomTypeMapper<DuneGrid, AllCellsLayout> Type;
72  };
73 
76  {
77  explicit CpGridCellMapper(const Dune::CpGrid&)
78  {
79  }
80  template<class EntityType>
81  int map (const EntityType& e) const
82  {
83  return e.index();
84  }
85  };
86 
90  template <>
91  struct GICellMapper<Dune::CpGrid>
92  {
93  typedef CpGridCellMapper Type;
94  };
95 
96 
97  namespace GIE
98  {
99  template <class GridInterface, class EntityPointerType>
100  class Cell;
101 
102  template <class GI>
103  class Face {
104  public:
105  typedef typename GI::DuneIntersectionIterator DuneIntersectionIter;
106  typedef typename GI::GridType::Traits::template Codim<0>::EntityPointer CellPtr;
107  typedef typename GI::GridType::ctype Scalar;
108  typedef Dune::FieldVector<Scalar, GI::GridType::dimension> Vector;
109  typedef Dune::FieldVector<Scalar, GI::GridType::dimension-1> LocalVector;
110  typedef int Index;
112 
113  enum { BoundaryMarkerIndex = -999,
114  LocalEndIndex = INT_MAX };
115 
116  Face()
117  : pgrid_(0), iter_(), local_index_(-1)
118  {
119  }
120 
121  Face(const GI& grid,
122  const DuneIntersectionIter& it,
123  const Index loc_ind)
124  : pgrid_(&grid), iter_(it), local_index_(loc_ind)
125  {
126  }
127 
131  Scalar area() const
132  {
133  return iter_->geometry().volume();
134  }
135 
139  Vector centroid() const
140  {
141  //return iter_->geometry().global(localCentroid());
142  return iter_->geometry().center();
143  }
144 
148  Vector normal() const
149  {
150  //return iter_->unitOuterNormal(localCentroid());
151  return iter_->centerUnitOuterNormal();
152  }
153 
157  bool boundary() const
158  {
159  return iter_->boundary();
160  }
161 
165  int boundaryId() const
166  {
167  return iter_->boundaryId();
168  }
169 
173  Cell cell() const
174  {
175  return Cell(*pgrid_, iter_->inside());
176  }
177 
181  Index cellIndex() const
182  {
183  return pgrid_->mapper().map(*iter_->inside());
184  }
185 
190  {
191  return Cell(*pgrid_, iter_->outside());
192  }
193 
197  Index neighbourCellIndex() const
198  {
199  if (iter_->boundary()) {
200  return BoundaryMarkerIndex;
201  } else {
202  return pgrid_->mapper().map(*iter_->outside());
203  }
204  }
205 
209  Index index() const
210  {
211  return pgrid_->faceIndex(cellIndex(), localIndex());
212  }
213 
217  Index localIndex() const
218  {
219  return local_index_;
220  }
221 
225  Scalar neighbourCellVolume() const
226  {
227  return iter_->outside()->geometry().volume();
228  }
229 
230  protected:
231  const GI* pgrid_;
232  DuneIntersectionIter iter_;
233  Index local_index_;
234 
235  private:
236 // LocalVector localCentroid() const
237 // {
238 // typedef Dune::ReferenceElements<Scalar, GI::GridType::dimension-1> RefElems;
239 // return RefElems::general(iter_->type()).position(0,0);
240 // }
241  };
242 
243 
249  template <class GridInterface>
250  class FaceIterator :
251  public boost::iterator_facade<FaceIterator<GridInterface>,
252  const Face<GridInterface>,
253  boost::forward_traversal_tag>,
254  public Face<GridInterface>
255  {
256  private:
258 
259  public:
264 
267  : FaceType()
268  {
269  }
270 
282  FaceIterator(const GridInterface& grid,
283  const DuneIntersectionIter& it,
284  const int local_index)
285  : FaceType(grid, it, local_index)
286  {
287  }
288 
290  const FaceIterator& dereference() const
291  {
292  return *this;
293  }
294 
296  bool equal(const FaceIterator& other) const
297  {
298  // Note that we do not compare the local_index_ members,
299  // since they may or may not be equal for end iterators.
300  return FaceType::iter_ == other.FaceType::iter_;
301  }
302 
304  void increment()
305  {
306  ++FaceType::iter_;
307  ++FaceType::local_index_;
308  }
309 
311  bool operator<(const FaceIterator& other) const
312  {
313  if (FaceType::cellIndex() == other.FaceType::cellIndex()) {
314  return FaceType::localIndex() < other.FaceType::localIndex();
315  } else {
316  return FaceType::cellIndex() < other.FaceType::cellIndex();
317  }
318  }
319  };
320 
321 
322  template <class GridInterface, class EntityPointerType>
323  class Cell
324  {
325  public:
326  Cell()
327  : pgrid_(0), iter_()
328  {
329  }
330  Cell(const GridInterface& grid,
331  const EntityPointerType& it)
332  : pgrid_(&grid), iter_(it)
333  {
334  }
335  typedef GIE::FaceIterator<GridInterface> FaceIterator;
336  typedef typename FaceIterator::Vector Vector;
337  typedef typename FaceIterator::Scalar Scalar;
338  typedef typename FaceIterator::Index Index;
339 
340  FaceIterator facebegin() const
341  {
342  return FaceIterator(*pgrid_, iter_->ileafbegin(), 0);
343  }
344 
345  FaceIterator faceend() const
346  {
347  return FaceIterator(*pgrid_, iter_->ileafend(),
348  FaceIterator::LocalEndIndex);
349  }
350 
351  Scalar volume() const
352  {
353  return iter_->geometry().volume();
354  }
355 
356  Vector centroid() const
357  {
358 // typedef Dune::ReferenceElements<Scalar, GridInterface::GridType::dimension> RefElems;
359 // Vector localpt
360 // = RefElems::general(iter_->type()).position(0,0);
361 // return iter_->geometry().global(localpt);
362  return iter_->geometry().center();
363  }
364 
365  Index index() const
366  {
367  return pgrid_->mapper().map(*iter_);
368  }
369  protected:
370  const GridInterface* pgrid_;
371  EntityPointerType iter_;
372  };
373 
374 
375  template <class GridInterface>
377  : public boost::iterator_facade<CellIterator<GridInterface>,
378  const Cell<GridInterface,
379  typename GridInterface::GridType::template Codim<0>::LeafIterator>,
380  boost::forward_traversal_tag>,
381  public Cell<GridInterface, typename GridInterface::GridType::template Codim<0>::LeafIterator>
382  {
383  private:
384  typedef typename GridInterface::GridType::template Codim<0>::LeafIterator DuneCellIter;
386  public:
387  typedef typename CellType::Vector Vector;
388  typedef typename CellType::Scalar Scalar;
389  typedef typename CellType::Index Index;
390 
391  CellIterator(const GridInterface& grid, DuneCellIter it)
392  : CellType(grid, it)
393  {
394  }
396  const CellIterator& dereference() const
397  {
398  return *this;
399  }
401  bool equal(const CellIterator& other) const
402  {
403  return CellType::iter_ == other.CellType::iter_;
404  }
406  void increment()
407  {
408  ++CellType::iter_;
409  }
410  };
411 
412  } // namespace GIE
413 
414 
415  template <class DuneGrid>
417  {
418  public:
419  typedef typename DuneGrid::LeafIntersectionIterator DuneIntersectionIterator;
420  typedef DuneGrid GridType;
423  typedef typename CellIterator::Vector Vector;
424  typedef typename CellIterator::Scalar Scalar;
425  typedef typename CellIterator::Index Index;
426 
427  typedef typename GICellMapper<DuneGrid>::Type Mapper;
428 
429  enum { Dimension = DuneGrid::dimension };
430 
432  : pgrid_(0), num_faces_(0), max_faces_per_cell_(0)
433  {
434  }
435  explicit GridInterfaceEuler(const DuneGrid& grid, bool build_facemap = true)
436  : pgrid_(&grid), pmapper_(new Mapper(grid)), num_faces_(0), max_faces_per_cell_(0)
437  {
438  if (build_facemap) {
439  buildFaceIndices();
440  }
441  }
442  void init(const DuneGrid& grid, bool build_facemap = true)
443  {
444  pgrid_ = &grid;
445  pmapper_.reset(new Mapper(grid));
446  if (build_facemap) {
447  buildFaceIndices();
448  }
449  }
450  CellIterator cellbegin() const
451  {
452  return CellIterator(*this, grid().template leafbegin<0>());
453  }
454  CellIterator cellend() const
455  {
456  return CellIterator(*this, grid().template leafend<0>());
457  }
458  int numberOfCells() const
459  {
460  return grid().size(0);
461  }
462  int numberOfFaces() const
463  {
464  assert(num_faces_ != 0);
465  return num_faces_;
466  }
467  int maxFacesPerCell() const
468  {
469  assert(max_faces_per_cell_ != 0);
470  return max_faces_per_cell_;
471  }
472  const DuneGrid& grid() const
473  {
474  assert(pgrid_);
475  return *pgrid_;
476  }
477 
478  // The following are primarily helpers for the implementation,
479  // perhaps they should be private?
480  const Mapper& mapper() const
481  {
482  assert (pmapper_);
483  return *pmapper_;
484  }
485  Index faceIndex(int cell_index, int local_face_index) const
486  {
487  assert(num_faces_ != 0);
488  return face_indices_[cell_index][local_face_index];
489  }
490  typedef Opm::SparseTable<int>::row_type Indices;
491  Indices faceIndices(int cell_index) const
492  {
493  assert(num_faces_ != 0);
494  return face_indices_[cell_index];
495  }
496  private:
497  const DuneGrid* pgrid_;
498  boost::scoped_ptr<Mapper> pmapper_;
499  int num_faces_;
500  int max_faces_per_cell_;
501  Opm::SparseTable<int> face_indices_;
502 
503  void buildFaceIndices()
504  {
505 #ifdef VERBOSE
506  std::cout << "Building unique face indices... " << std::flush;
507  Opm::time::StopWatch clock;
508  clock.start();
509 #endif
510  typedef CellIterator CI;
511  typedef typename CI::FaceIterator FI;
512 
513  // We build the actual cell to face mapping in two passes.
514  // [code mostly lifted from IncompFlowSolverHybrid::enumerateGridDof(),
515  // but with a twist: This code builds a mapping from cells in index
516  // order to unique face numbers, while the mapping built in the
517  // enumerateGridDof() method was ordered by cell iterator order]
518 
519  // Allocate and reserve structures.
520  const int nc = numberOfCells();
521  std::vector<int> cell(nc, -1);
522  std::vector<int> num_faces(nc); // In index order.
523  std::vector<int> fpos; fpos .reserve(nc + 1);
524  std::vector<int> num_cf; num_cf.reserve(nc); // In iterator order.
525  std::vector<int> faces ;
526 
527  // First pass: enumerate internal faces.
528  int cellno = 0; fpos.push_back(0);
529  int tot_ncf = 0, tot_ncf2 = 0, max_ncf = 0;
530  for (CI c = cellbegin(); c != cellend(); ++c, ++cellno) {
531  const int c0 = c->index();
532  assert((0 <= c0) && (c0 < nc) && (cell[c0] == -1));
533  cell[c0] = cellno;
534  num_cf.push_back(0);
535  int& ncf = num_cf.back();
536  for (FI f = c->facebegin(); f != c-> faceend(); ++f) {
537  if (!f->boundary()) {
538  const int c1 = f->neighbourCellIndex();
539  assert((0 <= c1) && (c1 < nc) && (c1 != c0));
540 
541  if (cell[c1] == -1) {
542  // Previously undiscovered internal face.
543  faces.push_back(c1);
544  }
545  }
546  ++ncf;
547  }
548  num_faces[c0] = ncf;
549  fpos.push_back(int(faces.size()));
550  max_ncf = std::max(max_ncf, ncf);
551  tot_ncf += ncf;
552  tot_ncf2 += ncf * ncf;
553  }
554  assert(cellno == nc);
555 
556  // Build cumulative face sizes enabling direct insertion of
557  // face indices into cfdata later.
558  std::vector<int> cumul_num_faces(numberOfCells() + 1);
559  cumul_num_faces[0] = 0;
560  std::partial_sum(num_faces.begin(), num_faces.end(), cumul_num_faces.begin() + 1);
561 
562  // Avoid (most) allocation(s) inside 'c' loop.
563  std::vector<int> l2g;
564  l2g.reserve(max_ncf);
565  std::vector<double> cfdata(tot_ncf);
566  int total_num_faces = int(faces.size());
567 
568  // Second pass: build cell-to-face mapping, including boundary.
569  typedef std::vector<int>::iterator VII;
570  for (CI c = cellbegin(); c != cellend(); ++c) {
571  const int c0 = c->index();
572  assert ((0 <= c0 ) && ( c0 < nc) &&
573  (0 <= cell[c0]) && (cell[c0] < nc));
574  const int ncf = num_cf[cell[c0]];
575  l2g.resize(ncf, 0);
576  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
577  if (f->boundary()) {
578  // External, not counted before. Add new face...
579  l2g[f->localIndex()] = total_num_faces++;
580  } else {
581  // Internal face. Need to determine during
582  // traversal of which cell we discovered this
583  // face first, and extract the face number
584  // from the 'faces' table range of that cell.
585 
586  // Note: std::find() below is potentially
587  // *VERY* expensive (e.g., large number of
588  // seeks in moderately sized data in case of
589  // faulted cells).
590  const int c1 = f->neighbourCellIndex();
591  assert ((0 <= c1 ) && ( c1 < nc) &&
592  (0 <= cell[c1]) && (cell[c1] < nc));
593 
594  int t = c0, seek = c1;
595  if (cell[seek] < cell[t])
596  std::swap(t, seek);
597  int s = fpos[cell[t]], e = fpos[cell[t] + 1];
598  VII p = std::find(faces.begin() + s, faces.begin() + e, seek);
599  assert(p != faces.begin() + e);
600  l2g[f->localIndex()] = p - faces.begin();
601  }
602  }
603  assert(int(l2g.size()) == num_faces[c0]);
604  std::copy(l2g.begin(), l2g.end(), cfdata.begin() + cumul_num_faces[c0]);
605  }
606  num_faces_ = total_num_faces;
607  max_faces_per_cell_ = max_ncf;
608  face_indices_.assign(cfdata.begin(), cfdata.end(),
609  num_faces.begin(), num_faces.end());
610 
611 #ifdef VERBOSE
612  clock.stop();
613  double elapsed = clock.secsSinceStart();
614  std::cout << "done. Time elapsed: " << elapsed << std::endl;
615 #endif
616  }
617 
618  };
619 
620 
621 
622 } // namespace Opm
623 
624 
625 #endif // OPENRS_GRIDINTERFACEEULER_HEADER
Vector normal() const
Definition: GridInterfaceEuler.hpp:148
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition: matrixops.hpp:33
void increment()
Used by iterator facade.
Definition: GridInterfaceEuler.hpp:304
Definition: GridInterfaceEuler.hpp:376
Cell cell() const
Definition: GridInterfaceEuler.hpp:173
bool equal(const CellIterator &other) const
Used by iterator facade.
Definition: GridInterfaceEuler.hpp:401
Definition: GridInterfaceEuler.hpp:100
Class for immiscible dead oil and dry gas.
Definition: applier.hpp:18
Definition: GridInterfaceEuler.hpp:416
bool equal(const FaceIterator &other) const
Used by iterator facade.
Definition: GridInterfaceEuler.hpp:296
Vector centroid() const
Definition: GridInterfaceEuler.hpp:139
Index cellIndex() const
Definition: GridInterfaceEuler.hpp:181
Scalar neighbourCellVolume() const
Definition: GridInterfaceEuler.hpp:225
General cell layout.
Definition: GridInterfaceEuler.hpp:63
Mapper for general grids.
Definition: GridInterfaceEuler.hpp:69
FaceIterator()
Default constructor.
Definition: GridInterfaceEuler.hpp:266
FaceIterator(const GridInterface &grid, const DuneIntersectionIter &it, const int local_index)
Constructor.
Definition: GridInterfaceEuler.hpp:282
Scalar area() const
Definition: GridInterfaceEuler.hpp:131
void increment()
Used by iterator facade.
Definition: GridInterfaceEuler.hpp:406
A mapper for Dune::CpGrid cells only.
Definition: GridInterfaceEuler.hpp:75
bool boundary() const
Definition: GridInterfaceEuler.hpp:157
Index neighbourCellIndex() const
Definition: GridInterfaceEuler.hpp:197
Intersection (face) iterator for solver-near grid interface.
Definition: GridInterfaceEuler.hpp:250
Index localIndex() const
Definition: GridInterfaceEuler.hpp:217
Definition: GridInterfaceEuler.hpp:103
bool operator<(const FaceIterator &other) const
Gives an ordering of intersectionIterators.
Definition: GridInterfaceEuler.hpp:311
Face< GridInterface >::DuneIntersectionIter DuneIntersectionIter
Type of low-level intersection iterator.
Definition: GridInterfaceEuler.hpp:263
const CellIterator & dereference() const
Used by iterator facade.
Definition: GridInterfaceEuler.hpp:396
Index index() const
Definition: GridInterfaceEuler.hpp:209
const FaceIterator & dereference() const
Used by iterator facade.
Definition: GridInterfaceEuler.hpp:290
Cell neighbourCell() const
Definition: GridInterfaceEuler.hpp:189
int boundaryId() const
Definition: GridInterfaceEuler.hpp:165