Loading...
Searching...
No Matches
ecfvstencil.hh
Go to the documentation of this file.
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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 2 of the License, or
(at your option) any later version.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#ifndef EWOMS_ECFV_STENCIL_HH
#define EWOMS_ECFV_STENCIL_HH
#include <dune/common/fvector.hh>
#include <dune/common/version.hh>
#include <dune/geometry/type.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/intersectioniterator.hh>
#include <opm/common/ErrorMacros.hpp>
#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
#include <opm/material/common/ConditionalStorage.hpp>
#include <opm/models/utils/quadraturegeometries.hh>
#include <cassert>
#include <cstddef>
#include <stdexcept>
#include <vector>
namespace Opm {
template <class Scalar,
class GridView,
bool needFaceIntegrationPos = true,
bool needFaceNormal = true>
class EcfvStencil
{
enum { dimWorld = GridView::dimensionworld };
using CoordScalar = typename GridView::ctype;
using Intersection = typename GridView::Intersection;
using Element = typename GridView::template Codim<0>::Entity;
using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
using WorldVector = Dune::FieldVector<Scalar, dimWorld>;
public:
using Entity = Element;
using Mapper = ElementMapper;
using LocalGeometry = typename Element::Geometry;
class SubControlVolume
{
public:
// default construct an uninitialized object.
// this is only here because std::vector needs it...
SubControlVolume()
{}
explicit SubControlVolume(const Element& element)
: element_(element)
{ update(); }
void update(const Element& element)
{ element_ = element; }
void update()
{}
decltype(auto) globalPos() const
{ return element_.geometry().center(); }
decltype(auto) center() const
{ return element_.geometry().center(); }
Scalar volume() const
{ return element_.geometry().volume(); }
const LocalGeometry geometry() const
{ return element_.geometry(); }
const LocalGeometry localGeometry() const
{ return element_.geometryInFather(); }
private:
Element element_;
};
template <bool needIntegrationPos, bool needNormal>
class EcfvSubControlVolumeFace
{
public:
EcfvSubControlVolumeFace()
{}
EcfvSubControlVolumeFace(const Intersection& intersection, unsigned localNeighborIdx)
{
exteriorIdx_ = static_cast<unsigned short>(localNeighborIdx);
if constexpr (needNormal) {
(*normal_) = intersection.centerUnitOuterNormal();
}
const auto& geometry = intersection.geometry();
if constexpr (needIntegrationPos) {
(*integrationPos_) = geometry.center();
}
area_ = geometry.volume();
// TODO: facedir_ = intersection.boundaryId(); // This did not compile for some reason
// using indexInInside() as in ecltransmissibility.cc instead
// Note that indexInInside() returns -1 for NNC faces, that's the reason we
// cannot convert directly to a FaceDir::DirEnum (see FaceDir.hpp)
dirId_ = intersection.indexInInside(); // Legal values: -1, 0, 1, 2, 3, 4, 5
}
unsigned short interiorIndex() const
{
// The local index of the control volume in the interior
// of a face of the stencil in the element centered finite
// volume discretization is always the "central"
// element. In this implementation this element always has
// index 0....
return 0;
}
unsigned short exteriorIndex() const
{ return exteriorIdx_; }
const GlobalPosition& integrationPos() const
{ return *integrationPos_; }
const WorldVector& normal() const
{ return *normal_; }
Scalar area() const
{ return area_; }
int dirId() const
{ return dirId_; }
FaceDir::DirEnum faceDirFromDirId() const
{
using Dir = FaceDir::DirEnum;
if (dirId_ == -1) {
OPM_THROW(std::runtime_error, "NNC faces does not have a face id");
}
switch(dirId_) {
case 0:
case 1:
return Dir::XPlus;
case 2:
case 3:
return Dir::YPlus;
case 4:
case 5:
return Dir::ZPlus;
default:
OPM_THROW(std::runtime_error,
"Unexpected face id" + std::to_string(dirId_));
}
}
private:
ConditionalStorage<needIntegrationPos, GlobalPosition> integrationPos_;
ConditionalStorage<needNormal, WorldVector> normal_;
Scalar area_;
int dirId_;
unsigned short exteriorIdx_;
};
using SubControlVolumeFace = EcfvSubControlVolumeFace<needFaceIntegrationPos, needFaceNormal>;
using BoundaryFace = EcfvSubControlVolumeFace</*needFaceIntegrationPos=*/true, needFaceNormal>;
EcfvStencil(const GridView& gridView, const Mapper& mapper)
: gridView_(gridView)
, elementMapper_(mapper)
{
// try to ensure that the mapper passed indeed maps elements
assert(int(gridView.size(/*codim=*/0)) == int(elementMapper_.size()));
}
void updateTopology(const Element& element)
{
// add the "center" element of the stencil
subControlVolumes_.clear();
subControlVolumes_.emplace_back(/*SubControlVolume(*/element/*)*/);
elements_.clear();
elements_.emplace_back(element);
interiorFaces_.clear();
boundaryFaces_.clear();
for (const auto& intersection : intersections(gridView_, element)) {
// if the current intersection has a neighbor, add a
// degree of freedom and an internal face, else add a
// boundary face
if (intersection.neighbor()) {
elements_.emplace_back(intersection.outside());
subControlVolumes_.emplace_back(/*SubControlVolume(*/elements_.back()/*)*/);
interiorFaces_.emplace_back(/*SubControlVolumeFace(*/intersection, subControlVolumes_.size() - 1/*)*/);
}
else {
boundaryFaces_.emplace_back(/*SubControlVolumeFace(*/intersection, - 10000/*)*/);
}
}
}
void updatePrimaryTopology(const Element& element)
{
// add the "center" element of the stencil
subControlVolumes_.clear();
subControlVolumes_.emplace_back(/*SubControlVolume(*/element/*)*/);
elements_.clear();
elements_.emplace_back(element);
}
void update(const Element& element)
{
updateTopology(element);
}
void updateCenterGradients()
{
assert(false); // not yet implemented
}
const Element& element() const
{ return element(0); }
const GridView& gridView() const
{ return *gridView_; }
std::size_t numDof() const
{ return subControlVolumes_.size(); }
std::size_t numPrimaryDof() const
{ return 1; }
unsigned globalSpaceIndex(unsigned dofIdx) const
{
assert(dofIdx < numDof());
return static_cast<unsigned>(elementMapper_.index(element(dofIdx)));
}
Dune::PartitionType partitionType(unsigned dofIdx) const
{ return elements_[dofIdx]->partitionType(); }
const Element& element(unsigned dofIdx) const
{
assert(dofIdx < numDof());
return elements_[dofIdx];
}
const Entity& entity(unsigned dofIdx) const
{ return element( dofIdx ); }
const SubControlVolume& subControlVolume(unsigned dofIdx) const
{ return subControlVolumes_[dofIdx]; }
std::size_t numInteriorFaces() const
{ return interiorFaces_.size(); }
const SubControlVolumeFace& interiorFace(unsigned faceIdx) const
{ return interiorFaces_[faceIdx]; }
std::size_t numBoundaryFaces() const
{ return boundaryFaces_.size(); }
const BoundaryFace& boundaryFace(unsigned bfIdx) const
{ return boundaryFaces_[bfIdx]; }
protected:
const GridView& gridView_;
const ElementMapper& elementMapper_;
std::vector<Element> elements_;
std::vector<SubControlVolume> subControlVolumes_;
std::vector<SubControlVolumeFace> interiorFaces_;
std::vector<BoundaryFace> boundaryFaces_;
};
} // namespace Opm
#endif
Opm::EcfvStencil::EcfvSubControlVolumeFace::integrationPos
const GlobalPosition & integrationPos() const
Returns the global position of the face's integration point.
Definition ecfvstencil.hh:202
Opm::EcfvStencil::EcfvSubControlVolumeFace::interiorIndex
unsigned short interiorIndex() const
Returns the local index of the degree of freedom to the face's interior.
Definition ecfvstencil.hh:181
Opm::EcfvStencil::EcfvSubControlVolumeFace::exteriorIndex
unsigned short exteriorIndex() const
Returns the local index of the degree of freedom to the face's outside.
Definition ecfvstencil.hh:195
Opm::EcfvStencil::EcfvSubControlVolumeFace::normal
const WorldVector & normal() const
Returns the outer unit normal at the face's integration point.
Definition ecfvstencil.hh:209
Opm::EcfvStencil::EcfvSubControlVolumeFace::faceDirFromDirId
FaceDir::DirEnum faceDirFromDirId() const
Returns the direction of the face.
Definition ecfvstencil.hh:230
Opm::EcfvStencil::EcfvSubControlVolumeFace::area
Scalar area() const
Returns the area [m^2] of the face.
Definition ecfvstencil.hh:215
Opm::EcfvStencil::EcfvSubControlVolumeFace::dirId
int dirId() const
Returns the direction id of the face w.r.t the cell.
Definition ecfvstencil.hh:224
Opm::EcfvStencil::SubControlVolume::geometry
const LocalGeometry geometry() const
The geometry of the sub-control volume.
Definition ecfvstencil.hh:135
Opm::EcfvStencil::SubControlVolume::volume
Scalar volume() const
The volume [m^3] occupied by the sub-control volume.
Definition ecfvstencil.hh:129
Opm::EcfvStencil::SubControlVolume::localGeometry
const LocalGeometry localGeometry() const
Geometry of the sub-control volume relative to parent.
Definition ecfvstencil.hh:141
Opm::EcfvStencil::SubControlVolume::center
decltype(auto) center() const
The center of the sub-control volume.
Definition ecfvstencil.hh:123
Opm::EcfvStencil::SubControlVolume::globalPos
decltype(auto) globalPos() const
The global position associated with the sub-control volume.
Definition ecfvstencil.hh:117
Opm::EcfvStencil
Represents the stencil (finite volume geometry) of a single element in the ECFV discretization.
Definition ecfvstencil.hh:69
Opm::EcfvStencil::interiorFace
const SubControlVolumeFace & interiorFace(unsigned faceIdx) const
Returns the face object belonging to a given face index in the interior of the domain.
Definition ecfvstencil.hh:404
Opm::EcfvStencil::boundaryFace
const BoundaryFace & boundaryFace(unsigned bfIdx) const
Returns the boundary face object belonging to a given boundary face index.
Definition ecfvstencil.hh:417
Opm::EcfvStencil::element
const Element & element(unsigned dofIdx) const
Return the element given the index of a degree of freedom.
Definition ecfvstencil.hh:373
Opm::EcfvStencil::element
const Element & element() const
Return the element to which the stencil refers.
Definition ecfvstencil.hh:320
Opm::EcfvStencil::entity
const Entity & entity(unsigned dofIdx) const
Return the entity given the index of a degree of freedom.
Definition ecfvstencil.hh:384
Opm::EcfvStencil::numDof
std::size_t numDof() const
Returns the number of degrees of freedom which the current element interacts with.
Definition ecfvstencil.hh:334
Opm::EcfvStencil::globalSpaceIndex
unsigned globalSpaceIndex(unsigned dofIdx) const
Return the global space index given the index of a degree of freedom.
Definition ecfvstencil.hh:353
Opm::EcfvStencil::numPrimaryDof
std::size_t numPrimaryDof() const
Returns the number of degrees of freedom which are contained by within the current element.
Definition ecfvstencil.hh:346
Opm::EcfvStencil::gridView
const GridView & gridView() const
Return the grid view of the element to which the stencil refers.
Definition ecfvstencil.hh:327
Opm::EcfvStencil::partitionType
Dune::PartitionType partitionType(unsigned dofIdx) const
Return partition type of a given degree of freedom.
Definition ecfvstencil.hh:363
Opm::EcfvStencil::numInteriorFaces
std::size_t numInteriorFaces() const
Returns the number of interior faces of the stencil.
Definition ecfvstencil.hh:397
Opm::EcfvStencil::numBoundaryFaces
std::size_t numBoundaryFaces() const
Returns the number of boundary faces of the stencil.
Definition ecfvstencil.hh:410
Opm::EcfvStencil::subControlVolume
const SubControlVolume & subControlVolume(unsigned dofIdx) const
Returns the sub-control volume object belonging to a given degree of freedom.
Definition ecfvstencil.hh:391
Opm
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilbioeffectsmodules.hh:43
quadraturegeometries.hh
Quadrature geometry for quadrilaterals.