My Project
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/grid/utility/SparseTable.hpp>
40 #include <opm/grid/utility/StopWatch.hpp>
41 #include <opm/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 #include <dune/grid/common/defaultgridview.hh>
48 
49 #include <boost/iterator/iterator_facade.hpp>
50 
51 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
52 
53 #include <climits>
54 #include <iostream>
55 #include <memory>
56 
57 
58 
59 namespace Opm
60 {
61 
63  template <class DuneGrid>
64  struct GICellMapper
65  {
66  using Type = Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<DuneGrid>>>;
67  };
68 
71  {
72  explicit CpGridCellMapper(const Dune::CpGrid&)
73  {
74  }
75  template<class EntityType>
76  int map (const EntityType& e) const
77  {
78  return e.index();
79  }
80  };
81 
82  namespace GIE
83  {
84  template <class GridInterface, class EntityPointerType>
85  class Cell;
86 
87  template <class GI>
88  class Face {
89  public:
90  typedef typename GI::DuneIntersectionIterator DuneIntersectionIter;
91  typedef typename GI::GridType::Traits::template Codim<0>::Entity CellPtr;
92  typedef typename GI::GridType::ctype Scalar;
93  typedef Dune::FieldVector<Scalar, GI::GridType::dimension> Vector;
94  typedef Dune::FieldVector<Scalar, GI::GridType::dimension-1> LocalVector;
95  typedef int Index;
97 
98  enum { BoundaryMarkerIndex = -999,
99  LocalEndIndex = INT_MAX };
100 
101  Face()
102  : pgrid_(0), iter_(), local_index_(-1)
103  {
104  }
105 
106  Face(const GI& grid,
107  const DuneIntersectionIter& it,
108  const Index loc_ind)
109  : pgrid_(&grid), iter_(it), local_index_(loc_ind)
110  {
111  }
112 
116  Scalar area() const
117  {
118  return iter_->geometry().volume();
119  }
120 
124  Vector centroid() const
125  {
126  //return iter_->geometry().global(localCentroid());
127  return iter_->geometry().center();
128  }
129 
133  Vector normal() const
134  {
135  //return iter_->unitOuterNormal(localCentroid());
136  return iter_->centerUnitOuterNormal();
137  }
138 
142  bool boundary() const
143  {
144  return iter_->boundary();
145  }
146 
150  int boundaryId() const
151  {
152  return iter_->boundaryId();
153  }
154 
158  Cell cell() const
159  {
160  return Cell(*pgrid_, iter_->inside());
161  }
162 
166  Index cellIndex() const
167  {
168  return pgrid_->mapper().index(iter_->inside());
169  }
170 
175  {
176  return Cell(*pgrid_, iter_->outside());
177  }
178 
182  Index neighbourCellIndex() const
183  {
184  if (iter_->boundary()) {
185  return BoundaryMarkerIndex;
186  } else {
187  return pgrid_->mapper().index(
188  iter_->outside());
189  }
190  }
191 
195  Index index() const
196  {
197  return pgrid_->faceIndex(cellIndex(), localIndex());
198  }
199 
203  Index localIndex() const
204  {
205  return local_index_;
206  }
207 
211  Scalar neighbourCellVolume() const
212  {
213  return iter_->outside()->geometry().volume();
214  }
215 
216  protected:
217  const GI* pgrid_;
218  DuneIntersectionIter iter_;
219  Index local_index_;
220 
221  private:
222 // LocalVector localCentroid() const
223 // {
224 // typedef Dune::ReferenceElements<Scalar, GI::GridType::dimension-1> RefElems;
225 // return RefElems::general(iter_->type()).position(0,0);
226 // }
227  };
228 
229 
235  template <class GridInterface>
236  class FaceIterator :
237  public boost::iterator_facade<FaceIterator<GridInterface>,
238  const Face<GridInterface>,
239  boost::forward_traversal_tag>,
240  public Face<GridInterface>
241  {
242  private:
244 
245  public:
250 
253  : FaceType()
254  {
255  }
256 
268  FaceIterator(const GridInterface& grid,
269  const DuneIntersectionIter& it,
270  const int local_index)
271  : FaceType(grid, it, local_index)
272  {
273  }
274 
276  const FaceIterator& dereference() const
277  {
278  return *this;
279  }
280 
282  bool equal(const FaceIterator& other) const
283  {
284  // Note that we do not compare the local_index_ members,
285  // since they may or may not be equal for end iterators.
286  return FaceType::iter_ == other.FaceType::iter_;
287  }
288 
290  void increment()
291  {
292  ++FaceType::iter_;
293  ++FaceType::local_index_;
294  }
295 
297  bool operator<(const FaceIterator& other) const
298  {
299  if (FaceType::cellIndex() == other.FaceType::cellIndex()) {
300  return FaceType::localIndex() < other.FaceType::localIndex();
301  } else {
302  return FaceType::cellIndex() < other.FaceType::cellIndex();
303  }
304  }
305  };
306 
307 
308  template <class GridInterface, class EntityType>
309  class Cell
310  {
311  public:
312  Cell()
313  : pgrid_(0), cell_()
314  {
315  }
316  Cell(const GridInterface& grid,
317  const EntityType& cell)
318  : pgrid_(&grid), cell_(cell)
319  {
320  }
322  typedef typename FaceIterator::Vector Vector;
323  typedef typename FaceIterator::Scalar Scalar;
324  typedef typename FaceIterator::Index Index;
325 
326  FaceIterator facebegin() const
327  {
328  return FaceIterator(*pgrid_, cell_.ileafbegin(), 0);
329  }
330 
331  FaceIterator faceend() const
332  {
333  return FaceIterator(*pgrid_, cell_.ileafend(),
334  FaceIterator::LocalEndIndex);
335  }
336 
337  Scalar volume() const
338  {
339  return cell_.geometry().volume();
340  }
341 
342  Vector centroid() const
343  {
344 // typedef Dune::ReferenceElements<Scalar, GridInterface::GridType::dimension> RefElems;
345 // Vector localpt
346 // = RefElems::general(cell_->type()).position(0,0);
347 // return cell_->geometry().global(localpt);
348  return cell_.geometry().center();
349  }
350 
351  Index index() const
352  {
353  return pgrid_->mapper().index(cell_);
354  }
355  protected:
356  const GridInterface* pgrid_;
357  EntityType cell_;
358  };
359 
360 
361  template <class GridInterface>
363  : public boost::iterator_facade<CellIterator<GridInterface>,
364  const Cell<GridInterface,
365  typename GridInterface::GridType::template Codim<0>::LeafIterator>,
366  boost::forward_traversal_tag>,
367  public Cell<GridInterface, typename GridInterface::GridType::template Codim<0>::LeafIterator>
368  {
369  private:
370  typedef typename GridInterface::GridType::template Codim<0>::LeafIterator DuneCellIter;
372  public:
373  typedef typename CellType::Vector Vector;
374  typedef typename CellType::Scalar Scalar;
375  typedef typename CellType::Index Index;
376 
377  CellIterator(const GridInterface& grid, DuneCellIter it)
378  : CellType(grid, it)
379  {
380  }
382  const CellIterator& dereference() const
383  {
384  return *this;
385  }
387  bool equal(const CellIterator& other) const
388  {
389  return CellType::cell_ == other.CellType::cell_;
390  }
392  void increment()
393  {
394  ++CellType::cell_;
395  }
396  };
397 
398  } // namespace GIE
399 
400 
401  template <class DuneGrid>
403  {
404  public:
405  typedef typename DuneGrid::LeafIntersectionIterator DuneIntersectionIterator;
406  typedef DuneGrid GridType;
409  typedef typename CellIterator::Vector Vector;
410  typedef typename CellIterator::Scalar Scalar;
411  typedef typename CellIterator::Index Index;
412 
413  typedef typename GICellMapper<DuneGrid>::Type Mapper;
414 
415  enum { Dimension = DuneGrid::dimension };
416 
418  : pgrid_(0), num_faces_(0), max_faces_per_cell_(0)
419  {
420  }
421  explicit GridInterfaceEuler(const DuneGrid& grid, bool build_facemap = true)
422  : pgrid_(&grid), pmapper_(new Mapper(grid.leafGridView(), Dune::mcmgElementLayout())), num_faces_(0), max_faces_per_cell_(0)
423  {
424  if (build_facemap) {
425  buildFaceIndices();
426  }
427  }
428  void init(const DuneGrid& grid, bool build_facemap = true)
429  {
430  pgrid_ = &grid;
431  pmapper_.reset(new Mapper(grid.leafGridView(), Dune::mcmgElementLayout()));
432  if (build_facemap) {
433  buildFaceIndices();
434  }
435  }
436  CellIterator cellbegin() const
437  {
438  return CellIterator(*this, grid().template leafbegin<0>());
439  }
440  CellIterator cellend() const
441  {
442  return CellIterator(*this, grid().template leafend<0>());
443  }
444  int numberOfCells() const
445  {
446  return grid().size(0);
447  }
448  int numberOfFaces() const
449  {
450  assert(num_faces_ != 0);
451  return num_faces_;
452  }
453  int maxFacesPerCell() const
454  {
455  assert(max_faces_per_cell_ != 0);
456  return max_faces_per_cell_;
457  }
458  const DuneGrid& grid() const
459  {
460  assert(pgrid_);
461  return *pgrid_;
462  }
463 
464  // The following are primarily helpers for the implementation,
465  // perhaps they should be private?
466  const Mapper& mapper() const
467  {
468  assert (pmapper_);
469  return *pmapper_;
470  }
471  Index faceIndex(int cell_index, int local_face_index) const
472  {
473  assert(num_faces_ != 0);
474  return face_indices_[cell_index][local_face_index];
475  }
476  typedef Opm::SparseTable<int>::row_type Indices;
477  Indices faceIndices(int cell_index) const
478  {
479  assert(num_faces_ != 0);
480  return face_indices_[cell_index];
481  }
482  private:
483  const DuneGrid* pgrid_;
484  std::unique_ptr<Mapper> pmapper_;
485  int num_faces_;
486  int max_faces_per_cell_;
487  Opm::SparseTable<int> face_indices_;
488 
489  void buildFaceIndices()
490  {
491 #ifdef VERBOSE
492  std::cout << "Building unique face indices... " << std::flush;
493  Opm::time::StopWatch clock;
494  clock.start();
495 #endif
496  typedef CellIterator CI;
497  typedef typename CI::FaceIterator FI;
498 
499  // We build the actual cell to face mapping in two passes.
500  // [code mostly lifted from IncompFlowSolverHybrid::enumerateGridDof(),
501  // but with a twist: This code builds a mapping from cells in index
502  // order to unique face numbers, while the mapping built in the
503  // enumerateGridDof() method was ordered by cell iterator order]
504 
505  // Allocate and reserve structures.
506  const int nc = numberOfCells();
507  std::vector<int> cell(nc, -1);
508  std::vector<int> num_faces(nc); // In index order.
509  std::vector<int> fpos; fpos .reserve(nc + 1);
510  std::vector<int> num_cf; num_cf.reserve(nc); // In iterator order.
511  std::vector<int> faces ;
512 
513  // First pass: enumerate internal faces.
514  int cellno = 0; fpos.push_back(0);
515  int tot_ncf = 0, max_ncf = 0;
516  for (CI c = cellbegin(); c != cellend(); ++c, ++cellno) {
517  const int c0 = c->index();
518  assert((0 <= c0) && (c0 < nc) && (cell[c0] == -1));
519  cell[c0] = cellno;
520  num_cf.push_back(0);
521  int& ncf = num_cf.back();
522  for (FI f = c->facebegin(); f != c-> faceend(); ++f) {
523  if (!f->boundary()) {
524  const int c1 = f->neighbourCellIndex();
525  assert((0 <= c1) && (c1 < nc) && (c1 != c0));
526 
527  if (cell[c1] == -1) {
528  // Previously undiscovered internal face.
529  faces.push_back(c1);
530  }
531  }
532  ++ncf;
533  }
534  num_faces[c0] = ncf;
535  fpos.push_back(int(faces.size()));
536  max_ncf = std::max(max_ncf, ncf);
537  tot_ncf += ncf;
538  }
539  assert(cellno == nc);
540 
541  // Build cumulative face sizes enabling direct insertion of
542  // face indices into cfdata later.
543  std::vector<int> cumul_num_faces(numberOfCells() + 1);
544  cumul_num_faces[0] = 0;
545  std::partial_sum(num_faces.begin(), num_faces.end(), cumul_num_faces.begin() + 1);
546 
547  // Avoid (most) allocation(s) inside 'c' loop.
548  std::vector<int> l2g;
549  l2g.reserve(max_ncf);
550  std::vector<double> cfdata(tot_ncf);
551  int total_num_faces = int(faces.size());
552 
553  // Second pass: build cell-to-face mapping, including boundary.
554  typedef std::vector<int>::iterator VII;
555  for (CI c = cellbegin(); c != cellend(); ++c) {
556  const int c0 = c->index();
557  assert ((0 <= c0 ) && ( c0 < nc) &&
558  (0 <= cell[c0]) && (cell[c0] < nc));
559  const int ncf = num_cf[cell[c0]];
560  l2g.resize(ncf, 0);
561  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
562  if (f->boundary()) {
563  // External, not counted before. Add new face...
564  l2g[f->localIndex()] = total_num_faces++;
565  } else {
566  // Internal face. Need to determine during
567  // traversal of which cell we discovered this
568  // face first, and extract the face number
569  // from the 'faces' table range of that cell.
570 
571  // Note: std::find() below is potentially
572  // *VERY* expensive (e.g., large number of
573  // seeks in moderately sized data in case of
574  // faulted cells).
575  const int c1 = f->neighbourCellIndex();
576  assert ((0 <= c1 ) && ( c1 < nc) &&
577  (0 <= cell[c1]) && (cell[c1] < nc));
578 
579  int t = c0, seek = c1;
580  if (cell[seek] < cell[t])
581  std::swap(t, seek);
582  int s = fpos[cell[t]], e = fpos[cell[t] + 1];
583  VII p = std::find(faces.begin() + s, faces.begin() + e, seek);
584  assert(p != faces.begin() + e);
585  l2g[f->localIndex()] = p - faces.begin();
586  }
587  }
588  assert(int(l2g.size()) == num_faces[c0]);
589  std::copy(l2g.begin(), l2g.end(), cfdata.begin() + cumul_num_faces[c0]);
590  }
591  num_faces_ = total_num_faces;
592  max_faces_per_cell_ = max_ncf;
593  face_indices_.assign(cfdata.begin(), cfdata.end(),
594  num_faces.begin(), num_faces.end());
595 
596 #ifdef VERBOSE
597  clock.stop();
598  double elapsed = clock.secsSinceStart();
599  std::cout << "done. Time elapsed: " << elapsed << std::endl;
600 #endif
601  }
602 
603  };
604 
605 
606 
607 } // namespace Opm
608 
609 
610 #endif // OPENRS_GRIDINTERFACEEULER_HEADER
Definition: GridInterfaceEuler.hpp:368
const CellIterator & dereference() const
Used by iterator facade.
Definition: GridInterfaceEuler.hpp:382
bool equal(const CellIterator &other) const
Used by iterator facade.
Definition: GridInterfaceEuler.hpp:387
void increment()
Used by iterator facade.
Definition: GridInterfaceEuler.hpp:392
Definition: GridInterfaceEuler.hpp:310
Intersection (face) iterator for solver-near grid interface.
Definition: GridInterfaceEuler.hpp:241
const FaceIterator & dereference() const
Used by iterator facade.
Definition: GridInterfaceEuler.hpp:276
bool equal(const FaceIterator &other) const
Used by iterator facade.
Definition: GridInterfaceEuler.hpp:282
void increment()
Used by iterator facade.
Definition: GridInterfaceEuler.hpp:290
Face< GridInterface >::DuneIntersectionIter DuneIntersectionIter
Type of low-level intersection iterator.
Definition: GridInterfaceEuler.hpp:249
bool operator<(const FaceIterator &other) const
Gives an ordering of intersectionIterators.
Definition: GridInterfaceEuler.hpp:297
FaceIterator(const GridInterface &grid, const DuneIntersectionIter &it, const int local_index)
Constructor.
Definition: GridInterfaceEuler.hpp:268
FaceIterator()
Default constructor.
Definition: GridInterfaceEuler.hpp:252
Definition: GridInterfaceEuler.hpp:88
Vector normal() const
Definition: GridInterfaceEuler.hpp:133
Vector centroid() const
Definition: GridInterfaceEuler.hpp:124
Index neighbourCellIndex() const
Definition: GridInterfaceEuler.hpp:182
bool boundary() const
Definition: GridInterfaceEuler.hpp:142
Cell cell() const
Definition: GridInterfaceEuler.hpp:158
Index cellIndex() const
Definition: GridInterfaceEuler.hpp:166
Cell neighbourCell() const
Definition: GridInterfaceEuler.hpp:174
int boundaryId() const
Definition: GridInterfaceEuler.hpp:150
Index localIndex() const
Definition: GridInterfaceEuler.hpp:203
Scalar neighbourCellVolume() const
Definition: GridInterfaceEuler.hpp:211
Scalar area() const
Definition: GridInterfaceEuler.hpp:116
Index index() const
Definition: GridInterfaceEuler.hpp:195
Definition: GridInterfaceEuler.hpp:403
Inverting small matrices.
Definition: ImplicitAssembly.hpp:43
A mapper for Dune::CpGrid cells only.
Definition: GridInterfaceEuler.hpp:71
Mapper for general grids.
Definition: GridInterfaceEuler.hpp:65