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;
98 enum { BoundaryMarkerIndex = -999,
99 LocalEndIndex = INT_MAX };
102 : pgrid_(0), iter_(), local_index_(-1)
107 const DuneIntersectionIter& it,
109 : pgrid_(&grid), iter_(it), local_index_(loc_ind)
118 return iter_->geometry().volume();
127 return iter_->geometry().center();
136 return iter_->centerUnitOuterNormal();
144 return iter_->boundary();
152 return iter_->boundaryId();
160 return Cell(*pgrid_, iter_->inside());
168 return pgrid_->mapper().index(iter_->inside());
176 return Cell(*pgrid_, iter_->outside());
184 if (iter_->boundary()) {
185 return BoundaryMarkerIndex;
187 return pgrid_->mapper().index(
213 return iter_->outside()->geometry().volume();
218 DuneIntersectionIter iter_;
316 Cell(
const GridInterface& grid,
317 const EntityType& cell)
318 : pgrid_(&grid), cell_(cell)
322 typedef typename FaceIterator::Vector Vector;
323 typedef typename FaceIterator::Scalar Scalar;
324 typedef typename FaceIterator::Index Index;
334 FaceIterator::LocalEndIndex);
337 Scalar volume()
const
339 return cell_.geometry().volume();
342 Vector centroid()
const
348 return cell_.geometry().center();
353 return pgrid_->mapper().index(cell_);
356 const GridInterface* pgrid_;
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;
413 typedef typename GICellMapper<DuneGrid>::Type Mapper;
415 enum { Dimension = DuneGrid::dimension };
418 : pgrid_(0), num_faces_(0), max_faces_per_cell_(0)
422 : pgrid_(&grid), pmapper_(
new Mapper(grid.leafGridView(), Dune::mcmgElementLayout())), num_faces_(0), max_faces_per_cell_(0)
428 void init(
const DuneGrid& grid,
bool build_facemap =
true)
431 pmapper_.reset(
new Mapper(grid.leafGridView(), Dune::mcmgElementLayout()));
438 return CellIterator(*
this, grid().
template leafbegin<0>());
442 return CellIterator(*
this, grid().
template leafend<0>());
444 int numberOfCells()
const
446 return grid().size(0);
448 int numberOfFaces()
const
450 assert(num_faces_ != 0);
453 int maxFacesPerCell()
const
455 assert(max_faces_per_cell_ != 0);
456 return max_faces_per_cell_;
458 const DuneGrid& grid()
const
466 const Mapper& mapper()
const
471 Index faceIndex(
int cell_index,
int local_face_index)
const
473 assert(num_faces_ != 0);
474 return face_indices_[cell_index][local_face_index];
476 typedef Opm::SparseTable<int>::row_type Indices;
477 Indices faceIndices(
int cell_index)
const
479 assert(num_faces_ != 0);
480 return face_indices_[cell_index];
483 const DuneGrid* pgrid_;
484 std::unique_ptr<Mapper> pmapper_;
486 int max_faces_per_cell_;
487 Opm::SparseTable<int> face_indices_;
489 void buildFaceIndices()
492 std::cout <<
"Building unique face indices... " << std::flush;
493 Opm::time::StopWatch clock;
497 typedef typename CI::FaceIterator FI;
506 const int nc = numberOfCells();
507 std::vector<int> cell(nc, -1);
508 std::vector<int> num_faces(nc);
509 std::vector<int> fpos; fpos .reserve(nc + 1);
510 std::vector<int> num_cf; num_cf.reserve(nc);
511 std::vector<int> 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));
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));
527 if (cell[c1] == -1) {
535 fpos.push_back(
int(faces.size()));
536 max_ncf = std::max(max_ncf, ncf);
539 assert(cellno == nc);
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);
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());
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]];
561 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
564 l2g[f->localIndex()] = total_num_faces++;
575 const int c1 = f->neighbourCellIndex();
576 assert ((0 <= c1 ) && ( c1 < nc) &&
577 (0 <= cell[c1]) && (cell[c1] < nc));
579 int t = c0, seek = c1;
580 if (cell[seek] < cell[t])
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();
588 assert(
int(l2g.size()) == num_faces[c0]);
589 std::copy(l2g.begin(), l2g.end(), cfdata.begin() + cumul_num_faces[c0]);
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());
598 double elapsed = clock.secsSinceStart();
599 std::cout <<
"done. Time elapsed: " << elapsed << std::endl;