Loading...
Searching...
No Matches
GridInterfaceEuler.hpp
//===========================================================================
//
// File: GridInterfaceEuler.hpp
//
// Created: Mon Jun 15 12:53:38 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <bard.skaflestad@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
OpenRS is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenRS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPENRS_GRIDINTERFACEEULER_HEADER
#define OPENRS_GRIDINTERFACEEULER_HEADER
#include <opm/grid/utility/SparseTable.hpp>
#include <opm/grid/utility/StopWatch.hpp>
#include <opm/grid/CpGrid.hpp> // How to avoid this? Needed for the explicit mapper specialization below.
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <dune/common/fvector.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/defaultgridview.hh>
#include <boost/iterator/iterator_facade.hpp>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <climits>
#include <iostream>
#include <memory>
namespace Opm
{
template <class DuneGrid>
struct GICellMapper
{
using Type = Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<DuneGrid>>>;
};
struct CpGridCellMapper
{
explicit CpGridCellMapper(const Dune::CpGrid&)
{
}
template<class EntityType>
int map (const EntityType& e) const
{
return e.index();
}
};
namespace GIE
{
template <class GridInterface, class EntityPointerType>
class Cell;
template <class GI>
class Face {
public:
typedef typename GI::DuneIntersectionIterator DuneIntersectionIter;
typedef typename GI::GridType::Traits::template Codim<0>::Entity CellPtr;
typedef typename GI::GridType::ctype Scalar;
typedef Dune::FieldVector<Scalar, GI::GridType::dimension> Vector;
typedef Dune::FieldVector<Scalar, GI::GridType::dimension-1> LocalVector;
typedef int Index;
typedef GIE::Cell<GI, CellPtr> Cell;
enum { BoundaryMarkerIndex = -999,
LocalEndIndex = INT_MAX };
Face()
: pgrid_(0), iter_(), local_index_(-1)
{
}
Face(const GI& grid,
const DuneIntersectionIter& it,
const Index loc_ind)
: pgrid_(&grid), iter_(it), local_index_(loc_ind)
{
}
Scalar area() const
{
return iter_->geometry().volume();
}
Vector centroid() const
{
//return iter_->geometry().global(localCentroid());
return iter_->geometry().center();
}
Vector normal() const
{
//return iter_->unitOuterNormal(localCentroid());
return iter_->centerUnitOuterNormal();
}
bool boundary() const
{
return iter_->boundary();
}
int boundaryId() const
{
return iter_->boundaryId();
}
Cell cell() const
{
return Cell(*pgrid_, iter_->inside());
}
Index cellIndex() const
{
return pgrid_->mapper().index(iter_->inside());
}
Cell neighbourCell() const
{
return Cell(*pgrid_, iter_->outside());
}
Index neighbourCellIndex() const
{
if (iter_->boundary()) {
return BoundaryMarkerIndex;
} else {
return pgrid_->mapper().index(
iter_->outside());
}
}
Index index() const
{
return pgrid_->faceIndex(cellIndex(), localIndex());
}
Index localIndex() const
{
return local_index_;
}
Scalar neighbourCellVolume() const
{
return iter_->outside()->geometry().volume();
}
protected:
const GI* pgrid_;
DuneIntersectionIter iter_;
Index local_index_;
private:
// LocalVector localCentroid() const
// {
// typedef Dune::ReferenceElements<Scalar, GI::GridType::dimension-1> RefElems;
// return RefElems::general(iter_->type()).position(0,0);
// }
};
template <class GridInterface>
class FaceIterator :
public boost::iterator_facade<FaceIterator<GridInterface>,
const Face<GridInterface>,
boost::forward_traversal_tag>,
public Face<GridInterface>
{
private:
typedef Face<GridInterface> FaceType;
public:
typedef typename Face<GridInterface>::DuneIntersectionIter DuneIntersectionIter;
FaceIterator()
: FaceType()
{
}
FaceIterator(const GridInterface& grid,
const DuneIntersectionIter& it,
const int local_index)
: FaceType(grid, it, local_index)
{
}
const FaceIterator& dereference() const
{
return *this;
}
bool equal(const FaceIterator& other) const
{
// Note that we do not compare the local_index_ members,
// since they may or may not be equal for end iterators.
return FaceType::iter_ == other.FaceType::iter_;
}
void increment()
{
++FaceType::iter_;
++FaceType::local_index_;
}
bool operator<(const FaceIterator& other) const
{
if (FaceType::cellIndex() == other.FaceType::cellIndex()) {
return FaceType::localIndex() < other.FaceType::localIndex();
} else {
return FaceType::cellIndex() < other.FaceType::cellIndex();
}
}
};
template <class GridInterface, class EntityType>
class Cell
{
public:
Cell()
: pgrid_(0), cell_()
{
}
Cell(const GridInterface& grid,
const EntityType& cell)
: pgrid_(&grid), cell_(cell)
{
}
typedef GIE::FaceIterator<GridInterface> FaceIterator;
typedef typename FaceIterator::Vector Vector;
typedef typename FaceIterator::Scalar Scalar;
typedef typename FaceIterator::Index Index;
FaceIterator facebegin() const
{
return FaceIterator(*pgrid_, cell_.ileafbegin(), 0);
}
FaceIterator faceend() const
{
return FaceIterator(*pgrid_, cell_.ileafend(),
FaceIterator::LocalEndIndex);
}
Scalar volume() const
{
return cell_.geometry().volume();
}
Vector centroid() const
{
// typedef Dune::ReferenceElements<Scalar, GridInterface::GridType::dimension> RefElems;
// Vector localpt
// = RefElems::general(cell_->type()).position(0,0);
// return cell_->geometry().global(localpt);
return cell_.geometry().center();
}
Index index() const
{
return pgrid_->mapper().index(cell_);
}
protected:
const GridInterface* pgrid_;
EntityType cell_;
};
template <class GridInterface>
class CellIterator
: public boost::iterator_facade<CellIterator<GridInterface>,
const Cell<GridInterface,
typename GridInterface::GridType::template Codim<0>::LeafIterator>,
boost::forward_traversal_tag>,
public Cell<GridInterface, typename GridInterface::GridType::template Codim<0>::LeafIterator>
{
private:
typedef typename GridInterface::GridType::template Codim<0>::LeafIterator DuneCellIter;
typedef Cell<GridInterface, DuneCellIter> CellType;
public:
typedef typename CellType::Vector Vector;
typedef typename CellType::Scalar Scalar;
typedef typename CellType::Index Index;
CellIterator(const GridInterface& grid, DuneCellIter it)
: CellType(grid, it)
{
}
const CellIterator& dereference() const
{
return *this;
}
bool equal(const CellIterator& other) const
{
return CellType::cell_ == other.CellType::cell_;
}
void increment()
{
++CellType::cell_;
}
};
} // namespace GIE
template <class DuneGrid>
class GridInterfaceEuler
{
public:
typedef typename DuneGrid::LeafIntersectionIterator DuneIntersectionIterator;
typedef DuneGrid GridType;
typedef GridInterfaceEuler<DuneGrid> InterfaceType;
typedef GIE::CellIterator<InterfaceType> CellIterator;
typedef typename CellIterator::Vector Vector;
typedef typename CellIterator::Scalar Scalar;
typedef typename CellIterator::Index Index;
typedef typename GICellMapper<DuneGrid>::Type Mapper;
enum { Dimension = DuneGrid::dimension };
GridInterfaceEuler()
: pgrid_(0), num_faces_(0), max_faces_per_cell_(0)
{
}
explicit GridInterfaceEuler(const DuneGrid& grid, bool build_facemap = true)
: pgrid_(&grid), pmapper_(new Mapper(grid.leafGridView(), Dune::mcmgElementLayout())), num_faces_(0), max_faces_per_cell_(0)
{
if (build_facemap) {
buildFaceIndices();
}
}
void init(const DuneGrid& grid, bool build_facemap = true)
{
pgrid_ = &grid;
pmapper_.reset(new Mapper(grid.leafGridView(), Dune::mcmgElementLayout()));
if (build_facemap) {
buildFaceIndices();
}
}
CellIterator cellbegin() const
{
return CellIterator(*this, grid().template leafbegin<0>());
}
CellIterator cellend() const
{
return CellIterator(*this, grid().template leafend<0>());
}
int numberOfCells() const
{
return grid().size(0);
}
int numberOfFaces() const
{
assert(num_faces_ != 0);
return num_faces_;
}
int maxFacesPerCell() const
{
assert(max_faces_per_cell_ != 0);
return max_faces_per_cell_;
}
const DuneGrid& grid() const
{
assert(pgrid_);
return *pgrid_;
}
// The following are primarily helpers for the implementation,
// perhaps they should be private?
const Mapper& mapper() const
{
assert (pmapper_);
return *pmapper_;
}
Index faceIndex(int cell_index, int local_face_index) const
{
assert(num_faces_ != 0);
return face_indices_[cell_index][local_face_index];
}
typedef Opm::SparseTable<int>::row_type Indices;
Indices faceIndices(int cell_index) const
{
assert(num_faces_ != 0);
return face_indices_[cell_index];
}
private:
const DuneGrid* pgrid_;
std::unique_ptr<Mapper> pmapper_;
int num_faces_;
int max_faces_per_cell_;
Opm::SparseTable<int> face_indices_;
void buildFaceIndices()
{
#ifdef VERBOSE
std::cout << "Building unique face indices... " << std::flush;
Opm::time::StopWatch clock;
clock.start();
#endif
typedef CellIterator CI;
typedef typename CI::FaceIterator FI;
// We build the actual cell to face mapping in two passes.
// [code mostly lifted from IncompFlowSolverHybrid::enumerateGridDof(),
// but with a twist: This code builds a mapping from cells in index
// order to unique face numbers, while the mapping built in the
// enumerateGridDof() method was ordered by cell iterator order]
// Allocate and reserve structures.
const int nc = numberOfCells();
std::vector<int> cell(nc, -1);
std::vector<int> num_faces(nc); // In index order.
std::vector<int> fpos; fpos .reserve(nc + 1);
std::vector<int> num_cf; num_cf.reserve(nc); // In iterator order.
std::vector<int> faces ;
// First pass: enumerate internal faces.
int cellno = 0; fpos.push_back(0);
int tot_ncf = 0, max_ncf = 0;
for (CI c = cellbegin(); c != cellend(); ++c, ++cellno) {
const int c0 = c->index();
assert((0 <= c0) && (c0 < nc) && (cell[c0] == -1));
cell[c0] = cellno;
num_cf.push_back(0);
int& ncf = num_cf.back();
for (FI f = c->facebegin(); f != c-> faceend(); ++f) {
if (!f->boundary()) {
const int c1 = f->neighbourCellIndex();
assert((0 <= c1) && (c1 < nc) && (c1 != c0));
if (cell[c1] == -1) {
// Previously undiscovered internal face.
faces.push_back(c1);
}
}
++ncf;
}
num_faces[c0] = ncf;
fpos.push_back(int(faces.size()));
max_ncf = std::max(max_ncf, ncf);
tot_ncf += ncf;
}
assert(cellno == nc);
// Build cumulative face sizes enabling direct insertion of
// face indices into cfdata later.
std::vector<int> cumul_num_faces(numberOfCells() + 1);
cumul_num_faces[0] = 0;
std::partial_sum(num_faces.begin(), num_faces.end(), cumul_num_faces.begin() + 1);
// Avoid (most) allocation(s) inside 'c' loop.
std::vector<int> l2g;
l2g.reserve(max_ncf);
std::vector<double> cfdata(tot_ncf);
int total_num_faces = int(faces.size());
// Second pass: build cell-to-face mapping, including boundary.
typedef std::vector<int>::iterator VII;
for (CI c = cellbegin(); c != cellend(); ++c) {
const int c0 = c->index();
assert ((0 <= c0 ) && ( c0 < nc) &&
(0 <= cell[c0]) && (cell[c0] < nc));
const int ncf = num_cf[cell[c0]];
l2g.resize(ncf, 0);
for (FI f = c->facebegin(); f != c->faceend(); ++f) {
if (f->boundary()) {
// External, not counted before. Add new face...
l2g[f->localIndex()] = total_num_faces++;
} else {
// Internal face. Need to determine during
// traversal of which cell we discovered this
// face first, and extract the face number
// from the 'faces' table range of that cell.
// Note: std::find() below is potentially
// *VERY* expensive (e.g., large number of
// seeks in moderately sized data in case of
// faulted cells).
const int c1 = f->neighbourCellIndex();
assert ((0 <= c1 ) && ( c1 < nc) &&
(0 <= cell[c1]) && (cell[c1] < nc));
int t = c0, seek = c1;
if (cell[seek] < cell[t])
std::swap(t, seek);
int s = fpos[cell[t]], e = fpos[cell[t] + 1];
VII p = std::find(faces.begin() + s, faces.begin() + e, seek);
assert(p != faces.begin() + e);
l2g[f->localIndex()] = p - faces.begin();
}
}
assert(int(l2g.size()) == num_faces[c0]);
std::copy(l2g.begin(), l2g.end(), cfdata.begin() + cumul_num_faces[c0]);
}
num_faces_ = total_num_faces;
max_faces_per_cell_ = max_ncf;
face_indices_.assign(cfdata.begin(), cfdata.end(),
num_faces.begin(), num_faces.end());
#ifdef VERBOSE
clock.stop();
double elapsed = clock.secsSinceStart();
std::cout << "done. Time elapsed: " << elapsed << std::endl;
#endif
}
};
} // namespace Opm
#endif // OPENRS_GRIDINTERFACEEULER_HEADER
Opm::GIE::CellIterator
Definition GridInterfaceEuler.hpp:368
Opm::GIE::CellIterator::dereference
const CellIterator & dereference() const
Used by iterator facade.
Definition GridInterfaceEuler.hpp:382
Opm::GIE::CellIterator::equal
bool equal(const CellIterator &other) const
Used by iterator facade.
Definition GridInterfaceEuler.hpp:387
Opm::GIE::CellIterator::increment
void increment()
Used by iterator facade.
Definition GridInterfaceEuler.hpp:392
Opm::GIE::Cell
Definition GridInterfaceEuler.hpp:310
Opm::GIE::FaceIterator
Intersection (face) iterator for solver-near grid interface.
Definition GridInterfaceEuler.hpp:241
Opm::GIE::FaceIterator::equal
bool equal(const FaceIterator &other) const
Used by iterator facade.
Definition GridInterfaceEuler.hpp:282
Opm::GIE::FaceIterator::increment
void increment()
Used by iterator facade.
Definition GridInterfaceEuler.hpp:290
Opm::GIE::FaceIterator::DuneIntersectionIter
Face< GridInterface >::DuneIntersectionIter DuneIntersectionIter
Type of low-level intersection iterator.
Definition GridInterfaceEuler.hpp:249
Opm::GIE::FaceIterator::operator<
bool operator<(const FaceIterator &other) const
Gives an ordering of intersectionIterators.
Definition GridInterfaceEuler.hpp:297
Opm::GIE::FaceIterator::dereference
const FaceIterator & dereference() const
Used by iterator facade.
Definition GridInterfaceEuler.hpp:276
Opm::GIE::FaceIterator::FaceIterator
FaceIterator(const GridInterface &grid, const DuneIntersectionIter &it, const int local_index)
Constructor.
Definition GridInterfaceEuler.hpp:268
Opm::GIE::FaceIterator::FaceIterator
FaceIterator()
Default constructor.
Definition GridInterfaceEuler.hpp:252
Opm::GIE::DuneIntersectionIter
Opm::GIE::Face::normal
Vector normal() const
Definition GridInterfaceEuler.hpp:133
Opm::GIE::Face::centroid
Vector centroid() const
Definition GridInterfaceEuler.hpp:124
Opm::GIE::Face::neighbourCellIndex
Index neighbourCellIndex() const
Definition GridInterfaceEuler.hpp:182
Opm::GIE::Face::boundary
bool boundary() const
Definition GridInterfaceEuler.hpp:142
Opm::GIE::Face::cell
Cell cell() const
Definition GridInterfaceEuler.hpp:158
Opm::GIE::Face::cellIndex
Index cellIndex() const
Definition GridInterfaceEuler.hpp:166
Opm::GIE::Face::neighbourCell
Cell neighbourCell() const
Definition GridInterfaceEuler.hpp:174
Opm::GIE::Face::boundaryId
int boundaryId() const
Definition GridInterfaceEuler.hpp:150
Opm::GIE::Face::localIndex
Index localIndex() const
Definition GridInterfaceEuler.hpp:203
Opm::GIE::Face::neighbourCellVolume
Scalar neighbourCellVolume() const
Definition GridInterfaceEuler.hpp:211
Opm::GIE::Face::area
Scalar area() const
Definition GridInterfaceEuler.hpp:116
Opm::GIE::Face::index
Index index() const
Definition GridInterfaceEuler.hpp:195
Opm::GridInterface
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Opm::GICellMapper
Mapper for general grids.
Definition GridInterfaceEuler.hpp:65