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