36 #ifndef OPENRS_GRIDINTERFACEEULER_HEADER 37 #define OPENRS_GRIDINTERFACEEULER_HEADER 39 #include <opm/core/utility/SparseTable.hpp> 40 #include <opm/core/utility/StopWatch.hpp> 41 #include <dune/grid/CpGrid.hpp> 43 #include <opm/common/utility/platform_dependent/disable_warnings.h> 45 #include <dune/common/fvector.hh> 46 #include <dune/grid/common/mcmgmapper.hh> 48 #include <boost/iterator/iterator_facade.hpp> 49 #include <boost/scoped_ptr.hpp> 51 #include <opm/common/utility/platform_dependent/reenable_warnings.h> 64 bool contains(Dune::GeometryType gt) {
return gt.dim() == dim; }
68 template <
class DuneGr
id>
71 typedef Dune::LeafMultipleCodimMultipleGeomTypeMapper<DuneGrid, AllCellsLayout> Type;
80 template<
class EntityType>
81 int map (
const EntityType& e)
const 99 template <
class Gr
idInterface,
class EntityPo
interType>
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;
113 enum { BoundaryMarkerIndex = -999,
114 LocalEndIndex = INT_MAX };
117 : pgrid_(0), iter_(), local_index_(-1)
122 const DuneIntersectionIter& it,
124 : pgrid_(&grid), iter_(it), local_index_(loc_ind)
133 return iter_->geometry().volume();
142 return iter_->geometry().center();
151 return iter_->centerUnitOuterNormal();
159 return iter_->boundary();
167 return iter_->boundaryId();
175 return Cell(*pgrid_, iter_->inside());
183 return pgrid_->mapper().map(*iter_->inside());
191 return Cell(*pgrid_, iter_->outside());
199 if (iter_->boundary()) {
200 return BoundaryMarkerIndex;
202 return pgrid_->mapper().map(*iter_->outside());
227 return iter_->outside()->geometry().volume();
232 DuneIntersectionIter iter_;
249 template <
class Gr
idInterface>
251 public boost::iterator_facade<FaceIterator<GridInterface>,
252 const Face<GridInterface>,
253 boost::forward_traversal_tag>,
254 public Face<GridInterface>
284 const int local_index)
300 return FaceType::iter_ == other.FaceType::iter_;
307 ++FaceType::local_index_;
322 template <
class Gr
idInterface,
class EntityPo
interType>
330 Cell(
const GridInterface& grid,
331 const EntityPointerType& it)
332 : pgrid_(&grid), iter_(it)
335 typedef GIE::FaceIterator<GridInterface> FaceIterator;
336 typedef typename FaceIterator::Vector Vector;
337 typedef typename FaceIterator::Scalar Scalar;
338 typedef typename FaceIterator::Index Index;
340 FaceIterator facebegin()
const 342 return FaceIterator(*pgrid_, iter_->ileafbegin(), 0);
345 FaceIterator faceend()
const 347 return FaceIterator(*pgrid_, iter_->ileafend(),
348 FaceIterator::LocalEndIndex);
351 Scalar volume()
const 353 return iter_->geometry().volume();
362 return iter_->geometry().center();
367 return pgrid_->mapper().map(*iter_);
370 const GridInterface* pgrid_;
371 EntityPointerType iter_;
375 template <
class Gr
idInterface>
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>
384 typedef typename GridInterface::GridType::template Codim<0>::LeafIterator DuneCellIter;
387 typedef typename CellType::Vector Vector;
388 typedef typename CellType::Scalar Scalar;
389 typedef typename CellType::Index Index;
391 CellIterator(
const GridInterface& grid, DuneCellIter it)
403 return CellType::iter_ == other.CellType::iter_;
415 template <
class DuneGr
id>
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;
427 typedef typename GICellMapper<DuneGrid>::Type Mapper;
429 enum { Dimension = DuneGrid::dimension };
432 : pgrid_(0), num_faces_(0), max_faces_per_cell_(0)
436 : pgrid_(&grid), pmapper_(
new Mapper(grid)), num_faces_(0), max_faces_per_cell_(0)
442 void init(
const DuneGrid& grid,
bool build_facemap =
true)
445 pmapper_.reset(
new Mapper(grid));
452 return CellIterator(*
this, grid().
template leafbegin<0>());
456 return CellIterator(*
this, grid().
template leafend<0>());
458 int numberOfCells()
const 460 return grid().size(0);
462 int numberOfFaces()
const 464 assert(num_faces_ != 0);
467 int maxFacesPerCell()
const 469 assert(max_faces_per_cell_ != 0);
470 return max_faces_per_cell_;
472 const DuneGrid& grid()
const 480 const Mapper& mapper()
const 485 Index faceIndex(
int cell_index,
int local_face_index)
const 487 assert(num_faces_ != 0);
488 return face_indices_[cell_index][local_face_index];
490 typedef Opm::SparseTable<int>::row_type Indices;
491 Indices faceIndices(
int cell_index)
const 493 assert(num_faces_ != 0);
494 return face_indices_[cell_index];
497 const DuneGrid* pgrid_;
498 boost::scoped_ptr<Mapper> pmapper_;
500 int max_faces_per_cell_;
501 Opm::SparseTable<int> face_indices_;
503 void buildFaceIndices()
506 std::cout <<
"Building unique face indices... " << std::flush;
507 Opm::time::StopWatch clock;
511 typedef typename CI::FaceIterator FI;
520 const int nc = numberOfCells();
521 std::vector<int> cell(nc, -1);
522 std::vector<int> num_faces(nc);
523 std::vector<int> fpos; fpos .reserve(nc + 1);
524 std::vector<int> num_cf; num_cf.reserve(nc);
525 std::vector<int> 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));
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));
541 if (cell[c1] == -1) {
549 fpos.push_back(
int(faces.size()));
550 max_ncf = std::max(max_ncf, ncf);
552 tot_ncf2 += ncf * ncf;
554 assert(cellno == nc);
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);
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());
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]];
576 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
579 l2g[f->localIndex()] = total_num_faces++;
590 const int c1 = f->neighbourCellIndex();
591 assert ((0 <= c1 ) && ( c1 < nc) &&
592 (0 <= cell[c1]) && (cell[c1] < nc));
594 int t = c0, seek = c1;
595 if (cell[seek] < cell[t])
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();
603 assert(
int(l2g.size()) == num_faces[c0]);
604 std::copy(l2g.begin(), l2g.end(), cfdata.begin() + cumul_num_faces[c0]);
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());
613 double elapsed = clock.secsSinceStart();
614 std::cout <<
"done. Time elapsed: " << elapsed << std::endl;
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