vcfvstencil.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_VCFV_STENCIL_HH
#define EWOMS_VCFV_STENCIL_HH

#include <dune/common/version.hh>
#include <dune/geometry/referenceelements.hh>

#include <dune/grid/common/intersectioniterator.hh>
#include <dune/grid/common/mcmgmapper.hh>

#include <opm/models/utils/quadraturegeometries.hh>

#if HAVE_DUNE_LOCALFUNCTIONS
#include <dune/localfunctions/lagrange/pqkfactory.hh>
#include <vector>
#endif // HAVE_DUNE_LOCALFUNCTIONS

#include <array>
#include <algorithm>
#include <cassert>
#include <cmath>
#include <stdexcept>
#include <string>
#include <utility>

namespace Opm {

enum class ElementType {
    none,
    simplex,
    cube,
};

template <class Scalar, unsigned dim, ElementType basicGeomType>
class VcfvScvGeometries;

// local geometries for 1D elements
template <class Scalar>
class VcfvScvGeometries<Scalar, /*dim=*/1, ElementType::cube>
{
    enum { dim = 1 };
    enum { numScv = 2 };

public:
    using ScvLocalGeometry = QuadrialteralQuadratureGeometry<Scalar, dim>;

    static void init()
    {
        // 1D LINE SEGMENTS
        const Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] = {
            { // corners of the first sub control volume
                { 0.0 },
                { 0.5 }
            },

            { // corners of the second sub control volume
                { 0.5 },
                { 1.0 }
            }
        };
        for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx) {
            scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
        }
    }

    static const ScvLocalGeometry& get(unsigned scvIdx)
    { return scvGeoms_[scvIdx]; }

private:
    static std::array<ScvLocalGeometry, numScv> scvGeoms_;
};

template <class Scalar>
class VcfvScvGeometries<Scalar, /*dim=*/1, ElementType::simplex>
{
    enum { dim = 1 };
    enum { numScv = 2 };

public:
    using ScvLocalGeometry = QuadrialteralQuadratureGeometry<Scalar, dim>;

    static const ScvLocalGeometry& get(unsigned)
    {
        throw std::logic_error("Not implemented: VcfvScvGeometries<Scalar, 1, ElementType::simplex>");
    }
};

// local geometries for 2D elements
template <class Scalar>
class VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::simplex>
{
    enum { dim = 2 };
    enum { numScv = 3 };

public:
    using ScvLocalGeometry = QuadrialteralQuadratureGeometry<Scalar, dim>;

    static const ScvLocalGeometry& get(unsigned scvIdx)
    { return scvGeoms_[scvIdx]; }

    static void init()
    {
        // 2D SIMPLEX
        const Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
        {
            { // SCV 0 corners
                { 0.0,   0.0 },
                { 1.0/2, 0.0 },
                { 1.0/3, 1.0/3 },
                { 0.0,   1.0/2 },
            },

            { // SCV 1 corners
                { 1.0/2, 0.0 },
                { 1.0,   0.0 },
                { 1.0/3, 1.0/3 },
                { 1.0/2, 1.0/2 },
            },

            { // SCV 2 corners
                { 0.0,   1.0/2 },
                { 1.0/3, 1.0/3 },
                { 0.0,   1.0 },
                { 1.0/2, 1.0/2 },
            }
        };

        for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx) {
            scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
        }
    }

private:
    static std::array<ScvLocalGeometry, numScv> scvGeoms_;
};

template <class Scalar>
class VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::cube>
{
    enum { dim = 2 };
    enum { numScv = 4 };

public:
    using ScvLocalGeometry = QuadrialteralQuadratureGeometry<Scalar, dim>;

    static const ScvLocalGeometry& get(unsigned scvIdx)
    { return scvGeoms_[scvIdx]; }

    static void init()
    {
        // 2D CUBE
        const Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
        {
            { // SCV 0 corners
                { 0.0, 0.0 },
                { 0.5, 0.0 },
                { 0.0, 0.5 },
                { 0.5, 0.5 }
            },

            { // SCV 1 corners
                { 0.5, 0.0 },
                { 1.0, 0.0 },
                { 0.5, 0.5 },
                { 1.0, 0.5 }
            },

            { // SCV 2 corners
                { 0.0,  0.5 },
                { 0.5,  0.5 },
                { 0.0,  1.0 },
                { 0.5,  1.0 }
            },

            { // SCV 3 corners
                { 0.5,   0.5 },
                { 1.0,   0.5 },
                { 0.5,   1.0 },
                { 1.0,   1.0 }
            }
        };

        for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx) {
            scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
        }
    }

    static std::array<ScvLocalGeometry, numScv> scvGeoms_;
};

// local geometries for 3D elements
template <class Scalar>
class VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::simplex>
{
    enum { dim = 3 };
    enum { numScv = 4 };

public:
    using ScvLocalGeometry = QuadrialteralQuadratureGeometry<Scalar, dim>;

    static const ScvLocalGeometry& get(unsigned scvIdx)
    { return scvGeoms_[scvIdx]; }

    static void init()
    {
        // 3D SIMPLEX
        const Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
        {
            { // SCV 0 corners
                { 0.0,   0.0,   0.0 },
                { 1.0/2, 0.0,   0.0 },
                { 0.0,   1.0/2, 0.0 },
                { 1.0/3, 1.0/3, 0.0 },

                { 0.0,   0.0,   0.5 },
                { 1.0/3, 0.0,   1.0/3 },
                { 0.0,   1.0/3, 1.0/3 },
                { 1.0/4, 1.0/4, 1.0/4 },
            },

            { // SCV 1 corners
                { 1.0/2, 0.0,   0.0 },
                { 1.0,   0.0,   0.0 },
                { 1.0/3, 1.0/3, 0.0 },
                { 1.0/2, 1.0/2, 0.0 },

                { 1.0/3, 0.0,   1.0/3 },
                { 1.0/2, 0.0,   1.0/2 },
                { 1.0/4, 1.0/4, 1.0/4 },
                { 1.0/3, 1.0/3, 1.0/3 },
            },

            { // SCV 2 corners
                { 0.0,   1.0/2, 0.0 },
                { 1.0/3, 1.0/3, 0.0 },
                { 0.0,   1.0,   0.0 },
                { 1.0/2, 1.0/2, 0.0 },

                { 0.0,   1.0/3, 1.0/3 },
                { 1.0/4, 1.0/4, 1.0/4 },
                { 0.0,   1.0/2, 1.0/2 },
                { 1.0/3, 1.0/3, 1.0/3 },
            },

            { // SCV 3 corners
                { 0.0,   0.0,   1.0/2 },
                { 1.0/3, 0.0,   1.0/3 },
                { 0.0,   1.0/3, 1.0/3 },
                { 1.0/4, 1.0/4, 1.0/4 },

                { 0.0,   0.0,   1.0   },
                { 1.0/2, 0.0,   1.0/2 },
                { 0.0,   1.0/2, 1.0/2 },
                { 1.0/3, 1.0/3, 1.0/3 },
            }
        };

        for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx) {
            scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx],  ScvLocalGeometry::numCorners);
        }
    }

private:
    static std::array<ScvLocalGeometry, numScv> scvGeoms_;
};

template <class Scalar>
class VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::cube>
{
    enum { dim = 3 };
    enum { numScv = 8 };

public:
    using ScvLocalGeometry = QuadrialteralQuadratureGeometry<Scalar, dim>;

    static const ScvLocalGeometry& get(unsigned scvIdx)
    { return scvGeoms_[scvIdx]; }

    static void init()
    {
        // 3D CUBE
        const Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
        {
            { // SCV 0 corners
                { 0.0,   0.0,   0.0 },
                { 1.0/2, 0.0,   0.0 },
                { 0.0,   1.0/2, 0.0 },
                { 1.0/2, 1.0/2, 0.0 },

                { 0.0,   0.0,   1.0/2 },
                { 1.0/2, 0.0,   1.0/2 },
                { 0.0,   1.0/2, 1.0/2 },
                { 1.0/2, 1.0/2, 1.0/2 },
            },

            { // SCV 1 corners
                { 1.0/2, 0.0,   0.0 },
                { 1.0,   0.0,   0.0 },
                { 1.0/2, 1.0/2, 0.0 },
                { 1.0,   1.0/2, 0.0 },

                { 1.0/2, 0.0,   1.0/2 },
                { 1.0,   0.0,   1.0/2 },
                { 1.0/2, 1.0/2, 1.0/2 },
                { 1.0,   1.0/2, 1.0/2 },
            },

            { // SCV 2 corners
                { 0.0,   1.0/2, 0.0 },
                { 1.0/2, 1.0/2, 0.0 },
                { 0.0,   1.0,   0.0 },
                { 1.0/2, 1.0,   0.0 },

                { 0.0,   1.0/2, 1.0/2 },
                { 1.0/2, 1.0/2, 1.0/2 },
                { 0.0,   1.0,   1.0/2 },
                { 1.0/2, 1.0,   1.0/2 },
            },

            { // SCV 3 corners
                { 1.0/2, 1.0/2, 0.0 },
                { 1.0,   1.0/2, 0.0 },
                { 1.0/2, 1.0,   0.0 },
                { 1.0,   1.0,   0.0 },

                { 1.0/2, 1.0/2, 1.0/2 },
                { 1.0,   1.0/2, 1.0/2 },
                { 1.0/2, 1.0,   1.0/2 },
                { 1.0,   1.0,   1.0/2 },
            },

            { // SCV 4 corners
                { 0.0,   0.0,   1.0/2 },
                { 1.0/2, 0.0,   1.0/2 },
                { 0.0,   1.0/2, 1.0/2 },
                { 1.0/2, 1.0/2, 1.0/2 },

                { 0.0,   0.0,   1.0 },
                { 1.0/2, 0.0,   1.0 },
                { 0.0,   1.0/2, 1.0 },
                { 1.0/2, 1.0/2, 1.0 },
            },

            { // SCV 5 corners
                { 1.0/2, 0.0,   1.0/2 },
                { 1.0,   0.0,   1.0/2 },
                { 1.0/2, 1.0/2, 1.0/2 },
                { 1.0,   1.0/2, 1.0/2 },

                { 1.0/2, 0.0,   1.0 },
                { 1.0,   0.0,   1.0 },
                { 1.0/2, 1.0/2, 1.0 },
                { 1.0,   1.0/2, 1.0 },
            },

            { // SCV 6 corners
                { 0.0,   1.0/2, 1.0/2 },
                { 1.0/2, 1.0/2, 1.0/2 },
                { 0.0,   1.0,   1.0/2 },
                { 1.0/2, 1.0,   1.0/2 },

                { 0.0,   1.0/2, 1.0 },
                { 1.0/2, 1.0/2, 1.0 },
                { 0.0,   1.0,   1.0 },
                { 1.0/2, 1.0,   1.0 },
            },

            { // SCV 7 corners
                { 1.0/2, 1.0/2, 1.0/2 },
                { 1.0,   1.0/2, 1.0/2 },
                { 1.0/2, 1.0,   1.0/2 },
                { 1.0,   1.0,   1.0/2 },

                { 1.0/2, 1.0/2, 1.0 },
                { 1.0,   1.0/2, 1.0 },
                { 1.0/2, 1.0,   1.0 },
                { 1.0,   1.0,   1.0 },
            },
        };

        for (unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx) {
            scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
        }
    }
private:
    static std::array<ScvLocalGeometry, numScv> scvGeoms_;
};

template <class Scalar, class GridView>
class VcfvStencil
{
    enum { dim = GridView::dimension };
    enum { dimWorld = GridView::dimensionworld };
    enum { maxNC = dim < 3 ? 4 : 8 };
    enum { maxNE = dim < 3 ? 4 : 12 };
    enum { maxNF = dim < 3 ? 1 : 6 };
    enum { maxBF = dim < 3 ? 8 : 24 };
    using CoordScalar = typename GridView::ctype;
    using Element = typename GridView::Traits::template Codim<0>::Entity;

public:
    using Entity = typename GridView::Traits::template Codim<dim>::Entity;

private:
    using Geometry = typename Element::Geometry;
    using DimVector = Dune::FieldVector<Scalar, dimWorld>;
    using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
    using LocalPosition = Dune::FieldVector<CoordScalar, dim>;
    using IntersectionIterator = typename GridView::IntersectionIterator;

    using ScvLocalGeometry = QuadrialteralQuadratureGeometry<Scalar, dim>;

#if HAVE_DUNE_LOCALFUNCTIONS
    using LocalFiniteElementCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
    using LocalFiniteElement = typename LocalFiniteElementCache::FiniteElementType;
    using LocalBasisTraits = typename LocalFiniteElement::Traits::LocalBasisType::Traits;
    using ShapeJacobian = typename LocalBasisTraits::JacobianType;
#endif // HAVE_DUNE_LOCALFUNCTIONS

    Scalar quadrilateralArea(const GlobalPosition& p0,
                             const GlobalPosition& p1,
                             const GlobalPosition& p2,
                             const GlobalPosition& p3)
    {
        return 0.5 * std::abs((p3[0] - p1[0]) * (p2[1] - p0[1]) -
                              (p3[1] - p1[1]) * (p2[0] - p0[0]));
    }

    void crossProduct(DimVector& c, const DimVector& a, const DimVector& b)
    {
        c[0] = a[1] * b[2] - a[2] * b[1];
        c[1] = a[2] * b[0] - a[0] * b[2];
        c[2] = a[0] * b[1] - a[1] * b[0];
    }

    Scalar pyramidVolume(const GlobalPosition& p0,
                         const GlobalPosition& p1,
                         const GlobalPosition& p2,
                         const GlobalPosition& p3,
                         const GlobalPosition& p4)
    {
        DimVector a(p2); a -= p0;
        DimVector b(p3); b -= p1;

        DimVector n;
        crossProduct(n, a, b);

        a = p4; a -= p0;

        return 1.0 / 6.0 * (n * a);
    }

    Scalar prismVolume(const GlobalPosition& p0,
                       const GlobalPosition& p1,
                       const GlobalPosition& p2,
                       const GlobalPosition& p3,
                       const GlobalPosition& p4,
                       const GlobalPosition& p5)
    {
        DimVector a(p4);
        for (unsigned k = 0; k < dimWorld; ++k) {
            a[k] -= p0[k];
        }
        DimVector b(p1);
        for (unsigned k = 0; k < dimWorld; ++k) {
            b[k] -= p3[k];
        }
        DimVector m;
        crossProduct(m, a, b);

        for (unsigned k = 0; k < dimWorld; ++k) {
            a[k] = p1[k] - p0[k];
        }
        for (unsigned k = 0; k < dimWorld; ++k) {
            b[k] = p2[k] - p0[k];
        }
        DimVector n;
        crossProduct(n, a, b);
        n += m;

        for (unsigned k = 0; k < dimWorld; ++k) {
            a[k] = p5[k] - p0[k];
        }

        return std::abs(1.0 / 6.0 * (n * a));
    }

    Scalar hexahedronVolume(const GlobalPosition& p0,
                            const GlobalPosition& p1,
                            const GlobalPosition& p2,
                            const GlobalPosition& p3,
                            const GlobalPosition& p4,
                            const GlobalPosition& p5,
                            const GlobalPosition& p6,
                            const GlobalPosition& p7)
    {
        return prismVolume(p0, p1, p2, p4, p5, p6) +
               prismVolume(p0, p2, p3, p4, p6, p7);
    }

    void normalOfQuadrilateral3D(DimVector& normal,
                                 const GlobalPosition& p0,
                                 const GlobalPosition& p1,
                                 const GlobalPosition& p2,
                                 const GlobalPosition& p3)
    {
        DimVector a(p2);
        for (unsigned k = 0; k < dimWorld; ++k) {
            a[k] -= p0[k];
        }
        DimVector b(p3);
        for (unsigned k = 0; k < dimWorld; ++k) {
            b[k] -= p1[k];
        }

        crossProduct(normal, a, b);
        normal *= 0.5;
    }

    Scalar quadrilateralArea3D(const GlobalPosition& p0,
                               const GlobalPosition& p1,
                               const GlobalPosition& p2,
                               const GlobalPosition& p3)
    {
        DimVector normal;
        normalOfQuadrilateral3D(normal, p0, p1, p2, p3);
        return normal.two_norm();
    }

    void getFaceIndices(unsigned numElemVertices, unsigned k, unsigned& leftFace, unsigned& rightFace)
    {
        static const unsigned edgeToFaceTet[2][6] = {
            {1, 0, 3, 2, 1, 3},
            {0, 2, 0, 1, 3, 2}
        };
        static const unsigned edgeToFacePyramid[2][8] = {
            {1, 2, 3, 4, 1, 3, 4, 2},
            {0, 0, 0, 0, 3, 2, 1, 4}
        };
        static const unsigned edgeToFacePrism[2][9] = {
            {1, 0, 2, 0, 1, 2, 4, 4, 4},
            {0, 2, 1, 3, 3, 3, 0, 1, 2}
        };
        static const unsigned edgeToFaceHex[2][12] = {
            {0, 2, 3, 1, 4, 1, 2, 4, 0, 5, 5, 3},
            {2, 1, 0, 3, 0, 4, 4, 3, 5, 1, 2, 5}
        };

        switch (numElemVertices) {
        case 4:
            leftFace = edgeToFaceTet[0][k];
            rightFace = edgeToFaceTet[1][k];
            break;
        case 5:
            leftFace = edgeToFacePyramid[0][k];
            rightFace = edgeToFacePyramid[1][k];
            break;
        case 6:
            leftFace = edgeToFacePrism[0][k];
            rightFace = edgeToFacePrism[1][k];
            break;
        case 8:
            leftFace = edgeToFaceHex[0][k];
            rightFace = edgeToFaceHex[1][k];
            break;
        default:
            throw std::logic_error("Not implemented: VcfvStencil::getFaceIndices for "
                                   + std::to_string(numElemVertices) + " vertices");
            break;
        }
    }

    void getEdgeIndices(unsigned numElemVertices, unsigned face, unsigned vert,
                        unsigned& leftEdge, unsigned& rightEdge)
    {
        static const int faceAndVertexToLeftEdgeTet[4][4] = {
                { 0, 0, 2, -1},
                { 0, 0, -1, 3},
                { 1, -1, 1, 3},
                {-1, 2, 2, 4}
        };
        static const int faceAndVertexToRightEdgeTet[4][4] = {
                { 1, 2, 1, -1},
                { 3, 4, -1, 4},
                { 3, -1, 5, 5},
                {-1, 4, 5, 5}
        };
        static const int faceAndVertexToLeftEdgePyramid[5][5] = {
            { 0, 2, 3, 1, -1},
            { 0, -1, 0, -1, 4},
            {-1, 1, -1, 1, 5},
            { 2, 2, -1, -1, 4},
            {-1, -1, 3, 3, 7}
        };
        static const int faceAndVertexToRightEdgePyramid[5][5] = {
            { 2, 1, 0, 3, -1},
            { 4, -1, 6, -1, 6},
            {-1, 5, -1, 7, 7},
            { 4, 5, -1, -1, 5},
            {-1, -1, 6, 7, 6}
        };
        static const int faceAndVertexToLeftEdgePrism[5][6] = {
            { 3, 3, -1, 0, 1, -1},
            { 4, -1, 4, 0, -1, 2},
            {-1, 5, 5, -1, 1, 2},
            { 3, 3, 5, -1, -1, -1},
            {-1, -1, -1, 6, 6, 8}
        };
        static const int faceAndVertexToRightEdgePrism[5][6] = {
            { 0, 1, -1, 6, 6, -1},
            { 0, -1, 2, 7, -1, 7},
            {-1, 1, 2, -1, 8, 8},
            { 4, 5, 4, -1, -1, -1},
            {-1, -1, -1, 7, 8, 7}
        };
        static const int faceAndVertexToLeftEdgeHex[6][8] = {
            { 0, -1, 4, -1, 8, -1, 2, -1},
            {-1, 5, -1, 3, -1, 1, -1, 9},
            { 6, 1, -1, -1, 0, 10, -1, -1},
            {-1, -1, 2, 7, -1, -1, 11, 3},
            { 4, 6, 7, 5, -1, -1, -1, -1},
            {-1, -1, -1, -1, 10, 9, 8, 11}
        };
        static const int faceAndVertexToRightEdgeHex[6][8] = {
            { 4, -1, 2, -1, 0, -1, 8, -1},
            {-1, 1, -1, 5, -1, 9, -1, 3},
            { 0, 6, -1, -1, 10, 1, -1, -1},
            {-1, -1, 7, 3, -1, -1, 2, 11},
            { 6, 5, 4, 7, -1, -1, -1, -1},
            {-1, -1, -1, -1, 8, 10, 11, 9}
        };

        switch (numElemVertices) {
        case 4:
            leftEdge = static_cast<unsigned>(faceAndVertexToLeftEdgeTet[face][vert]);
            rightEdge = static_cast<unsigned>(faceAndVertexToRightEdgeTet[face][vert]);
            break;
        case 5:
            leftEdge = static_cast<unsigned>(faceAndVertexToLeftEdgePyramid[face][vert]);
            rightEdge = static_cast<unsigned>(faceAndVertexToRightEdgePyramid[face][vert]);
            break;
        case 6:
            leftEdge = static_cast<unsigned>(faceAndVertexToLeftEdgePrism[face][vert]);
            rightEdge = static_cast<unsigned>(faceAndVertexToRightEdgePrism[face][vert]);
            break;
        case 8:
            leftEdge = static_cast<unsigned>(faceAndVertexToLeftEdgeHex[face][vert]);
            rightEdge = static_cast<unsigned>(faceAndVertexToRightEdgeHex[face][vert]);
            break;
        default:
            throw std::logic_error("Not implemented: VcfvStencil::getFaceIndices for "
                                   +std::to_string(numElemVertices)+" vertices");
            break;
        }
    }

public:
    using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;

    class ScvGeometry
    {
    public:
        const GlobalPosition center() const
        { return global(localGeometry_->center()); }

        const GlobalPosition corner(unsigned cornerIdx) const
        { return global(localGeometry_->corner(cornerIdx)); }

        const GlobalPosition global(const LocalPosition& localPos) const
        { return element_->geometry().global(localPos); }

        const ScvLocalGeometry& localGeometry() const
        { return *localGeometry_; }

        const ScvLocalGeometry* localGeometry_;
        const Element* element_;
    };

    struct SubControlVolume 
    {
        const GlobalPosition& globalPos() const
        { return global; }

        const GlobalPosition center() const
        { return geometry_.center(); }

        Scalar volume() const
        { return volume_; }

        const ScvLocalGeometry& localGeometry() const
        { return geometry_.localGeometry(); }

        const ScvGeometry& geometry() const
        { return geometry_; }

        LocalPosition local;

        GlobalPosition global;

        Scalar volume_;

        ScvGeometry geometry_;

        Dune::FieldVector<DimVector, maxNC> gradCenter;
    };

    struct SubControlVolumeFace 
    {
        const DimVector& normal() const
        { return normal_; }

        unsigned short interiorIndex() const
        { return i; }

        unsigned short exteriorIndex() const
        { return j; }

        Scalar area() const
        { return area_; }

        const LocalPosition& localPos() const
        { return ipLocal_; }

        const GlobalPosition& integrationPos() const
        { return ipGlobal_; }

        unsigned short i,j;

        LocalPosition ipLocal_;

        GlobalPosition ipGlobal_;

        DimVector normal_;

        Scalar area_;
    };

    using BoundaryFace = SubControlVolumeFace;

    VcfvStencil(const GridView& gridView, const Mapper& mapper)
        : gridView_(gridView)
        , vertexMapper_(mapper )
        , element_(*gridView.template begin</*codim=*/0>())
    {
        // try to check if the mapper really maps the vertices
        assert(static_cast<int>(gridView.size(/*codim=*/dimWorld)) == static_cast<int>(mapper.size()));

        static bool localGeometriesInitialized = false;
        if (!localGeometriesInitialized) {
            localGeometriesInitialized = true;

            VcfvScvGeometries<Scalar, /*dim=*/1, ElementType::cube>::init();
            VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::cube>::init();
            VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::simplex>::init();
            VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::cube>::init();
            VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::simplex>::init();
        }
    }

    void updateTopology(const Element& e)
    {
        element_ = e;

        numVertices = e.subEntities(/*codim=*/dim);
        numEdges = e.subEntities(/*codim=*/dim-1);
        if constexpr (dim == 3) {
            numFaces = e.subEntities(/*codim=*/1);
        }
        else {
            numFaces = 0;
        }

        numBoundarySegments_ = 0; // TODO: really required here(?)

        // compute the local and global coordinates of the element
        const Geometry& geometry = e.geometry();
        geometryType_ = geometry.type();
        const auto& referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
        // should be the same, but let's let the compiler know to quell (false) diagnostics
        assert(numVertices == static_cast<unsigned>(referenceElement.size(dim)));
        numVertices = std::min(numVertices, static_cast<unsigned>(referenceElement.size(dim)));
        for (unsigned vertexIdx = 0; vertexIdx < numVertices; ++vertexIdx) {
            subContVol[vertexIdx].local = referenceElement.position(static_cast<int>(vertexIdx), dim);
            subContVol[vertexIdx].global = geometry.corner(static_cast<int>(vertexIdx));
        }
    }

    void updatePrimaryTopology(const Element& element)
    {
        // since all degrees of freedom in a stencil are "primary" DOFs for the
        // vertex-centered finite volume method, there's no difference to
        // updateTopology()
        updateTopology(element);
    }

    void update(const Element& e)
    {
        updateTopology(e);

        const Geometry& geometry = e.geometry();
        geometryType_ = geometry.type();

        const auto& referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);

        elementVolume = geometry.volume();
        elementLocal = referenceElement.position(0, 0);
        elementGlobal = geometry.global(elementLocal);

        // corners:
        for (unsigned vert = 0; vert < numVertices; ++vert) {
            subContVol[vert].local = referenceElement.position(static_cast<int>(vert), dim);
            subContVol[vert].global = geometry.global(subContVol[vert].local);
        }

        // edges:
        for (unsigned edge = 0; edge < numEdges; ++edge) {
            edgeCoord[edge] = geometry.global(referenceElement.position(static_cast<int>(edge), dim - 1));
        }

        // faces:
        for (unsigned face = 0; face < numFaces; ++face) {
            faceCoord[face] = geometry.global(referenceElement.position(static_cast<int>(face), 1));
        }

        // fill sub control volume data use specialization for this
        // \todo maybe it would be a good idea to migrate everything
        // which is dependend of the grid's dimension to
        // _VcfvFVElemGeomHelper in order to benefit from more aggressive
        // compiler optimizations...
        fillSubContVolData_();

        // fill sub control volume face data:
        for (unsigned k = 0; k < numEdges; ++k) { // begin loop over edges / sub control volume faces
            unsigned short i =
                static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(k), dim - 1, 0, dim));
            unsigned short j =
                static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(k), dim - 1, 1, dim));
            if (numEdges == 4 && (i == 2 || j == 2)) {
                std::swap(i, j);
            }
            subContVolFace[k].i = i;
            subContVolFace[k].j = j;

            // calculate the local integration point and
            // the face normal. note that since dim is a
            // constant which is known at compile time
            // the compiler can optimize away all if
            // cases which don't apply.
            LocalPosition ipLocal;
            DimVector diffVec;
            if constexpr (dim == 1) {
                subContVolFace[k].ipLocal_ = 0.5;
                subContVolFace[k].normal_ = 1.0;
                subContVolFace[k].area_ = 1.0;
                ipLocal = subContVolFace[k].ipLocal_;
            }
            else if constexpr (dim == 2) {
                ipLocal = referenceElement.position(static_cast<int>(k), dim - 1) + elementLocal;
                ipLocal *= 0.5;
                subContVolFace[k].ipLocal_ = ipLocal;
                for (unsigned m = 0; m < dimWorld; ++m) {
                    diffVec[m] = elementGlobal[m] - edgeCoord[k][m];
                }
                subContVolFace[k].normal_[0] = diffVec[1];
                subContVolFace[k].normal_[1] = -diffVec[0];

                for (unsigned m = 0; m < dimWorld; ++m) {
                    diffVec[m] = subContVol[j].global[m] - subContVol[i].global[m];
                }
                // make sure the normal points to the right direction
                if (subContVolFace[k].normal_ * diffVec < 0) {
                    subContVolFace[k].normal_ *= -1;
                }

                subContVolFace[k].area_ = subContVolFace[k].normal_.two_norm();
                subContVolFace[k].normal_ /= subContVolFace[k].area_;
            }
            else if constexpr (dim == 3) {
                unsigned leftFace;
                unsigned rightFace;
                getFaceIndices(numVertices, k, leftFace, rightFace);
                ipLocal = referenceElement.position(static_cast<int>(k), dim - 1) +
                          elementLocal +
                          referenceElement.position(static_cast<int>(leftFace), 1) +
                          referenceElement.position(static_cast<int>(rightFace), 1);
                ipLocal *= 0.25;
                subContVolFace[k].ipLocal_ = ipLocal;
                normalOfQuadrilateral3D(subContVolFace[k].normal_,
                                        edgeCoord[k], faceCoord[rightFace],
                                        elementGlobal, faceCoord[leftFace]);
                subContVolFace[k].area_ = subContVolFace[k].normal_.two_norm();
                subContVolFace[k].normal_ /= subContVolFace[k].area_;
            }

            // get the global integration point and the Jacobian inverse
            subContVolFace[k].ipGlobal_ = geometry.global(ipLocal);
        } // end loop over edges / sub control volume faces

        // fill boundary face data:
        for (const auto& intersection : intersections(gridView_, e)) {
            if (!intersection.boundary()) {
                continue;
            }

            const unsigned face = static_cast<unsigned>(intersection.indexInInside());
            const unsigned numVerticesOfFace = static_cast<unsigned>(referenceElement.size(static_cast<int>(face), 1, dim));
            for (unsigned vertInFace = 0; vertInFace < numVerticesOfFace; ++vertInFace) {
                const unsigned short vertInElement = static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(face), 1, static_cast<int>(vertInFace), dim));
                const unsigned bfIdx = numBoundarySegments_;
                ++numBoundarySegments_;

                if constexpr (dim == 1) {
                    boundaryFace_[bfIdx].ipLocal_ = referenceElement.position(static_cast<int>(vertInElement), dim);
                    boundaryFace_[bfIdx].area_ = 1.0;
                }
                else if constexpr (dim == 2) {
                    boundaryFace_[bfIdx].ipLocal_ =
                        referenceElement.position(static_cast<int>(vertInElement), dim) +
                        referenceElement.position(static_cast<int>(face), 1);
                    boundaryFace_[bfIdx].ipLocal_ *= 0.5;
                    boundaryFace_[bfIdx].area_ = 0.5 * intersection.geometry().volume();
                }
                else if constexpr (dim == 3) {
                    unsigned leftEdge;
                    unsigned rightEdge;
                    getEdgeIndices(numVertices, face, vertInElement, leftEdge, rightEdge);
                    boundaryFace_[bfIdx].ipLocal_ =
                        referenceElement.position(static_cast<int>(vertInElement), dim) +
                        referenceElement.position(static_cast<int>(face), 1) +
                        referenceElement.position(static_cast<int>(leftEdge), dim - 1) +
                        referenceElement.position(static_cast<int>(rightEdge), dim - 1);
                    boundaryFace_[bfIdx].ipLocal_ *= 0.25;
                    boundaryFace_[bfIdx].area_ =
                        quadrilateralArea3D(subContVol[vertInElement].global,
                                            edgeCoord[rightEdge],
                                            faceCoord[face],
                                            edgeCoord[leftEdge]);
                }
                else {
                    throw std::logic_error("Not implemented:VcfvStencil for dim = " + std::to_string(dim));
                }

                boundaryFace_[bfIdx].ipGlobal_ = geometry.global(boundaryFace_[bfIdx].ipLocal_);
                boundaryFace_[bfIdx].i = vertInElement;
                boundaryFace_[bfIdx].j = vertInElement;

                // ASSUME constant normal on the segment of the boundary face
                boundaryFace_[bfIdx].normal_ = intersection.centerUnitOuterNormal();
            }
        }

        updateScvGeometry(e);
    }

    void updateScvGeometry(const Element& element)
    {
        auto geomType = element.geometry().type();

        // get the local geometries of the sub control volumes
        if (geomType.isTriangle() || geomType.isTetrahedron()) {
            for (unsigned vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
                subContVol[vertIdx].geometry_.element_ = &element;
                subContVol[vertIdx].geometry_.localGeometry_ =
                    &VcfvScvGeometries<Scalar, dim, ElementType::simplex>::get(vertIdx);
            }
        }
        else if (geomType.isLine() || geomType.isQuadrilateral() || geomType.isHexahedron()) {
            for (unsigned vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
                subContVol[vertIdx].geometry_.element_ = &element;
                subContVol[vertIdx].geometry_.localGeometry_ =
                    &VcfvScvGeometries<Scalar, dim, ElementType::cube>::get(vertIdx);
            }
        }
        else {
            throw std::logic_error("Not implemented: SCV geometries for non hexahedron elements");
        }
    }

#if HAVE_DUNE_LOCALFUNCTIONS
    void updateCenterGradients()
    {
        const auto& localFiniteElement = feCache_.get(element_.type());
        const auto& geom = element_.geometry();

        std::vector<ShapeJacobian> localJac;

        for (unsigned scvIdx = 0; scvIdx < numVertices; ++scvIdx) {
            const auto& localCenter = subContVol[scvIdx].localGeometry().center();
            localFiniteElement.localBasis().evaluateJacobian(localCenter, localJac);
            const auto& globalPos = subContVol[scvIdx].geometry().center();

            const auto jacInvT = geom.jacobianInverseTransposed(globalPos);
            for (unsigned vert = 0; vert < numVertices; ++vert) {
                jacInvT.mv(localJac[vert][0], subContVol[scvIdx].gradCenter[vert]);
            }
        }
    }
#endif

    unsigned numDof() const
    { return numVertices; }

    unsigned numPrimaryDof() const
    { return numDof(); }

    Dune::PartitionType partitionType(unsigned scvIdx) const
    { return element_.template subEntity</*codim=*/dim>(scvIdx)->partitionType(); }

    const SubControlVolume& subControlVolume(unsigned scvIdx) const
    {
        assert(scvIdx < numDof());
        return subContVol[scvIdx];
    }

    unsigned numInteriorFaces() const
    { return numEdges; }

    unsigned numBoundaryFaces() const
    { return numBoundarySegments_; }

    const SubControlVolumeFace& interiorFace(unsigned faceIdx) const
    { return subContVolFace[faceIdx]; }

    const BoundaryFace& boundaryFace(unsigned bfIdx) const
    { return boundaryFace_[bfIdx]; }

    unsigned globalSpaceIndex(unsigned dofIdx) const
    {
        assert(dofIdx < numDof());

        return static_cast<unsigned>(vertexMapper_.subIndex(element_, static_cast<int>(dofIdx), /*codim=*/dim));
    }

    Entity entity(unsigned dofIdx) const
    {
        assert(dofIdx < numDof());
        return element_.template subEntity<dim>(static_cast<int>(dofIdx));
    }

private:
#if __GNUC__ || __clang__
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wpragmas"
#pragma GCC diagnostic ignored "-Warray-bounds"
#endif
    void fillSubContVolData_()
    {
        if constexpr (dim == 1) {
            // 1D
            subContVol[0].volume_ = 0.5 * elementVolume;
            subContVol[1].volume_ = 0.5 * elementVolume;
        }
        else if constexpr (dim == 2) {
            switch (numVertices) {
            case 3: // 2D, triangle
                subContVol[0].volume_ = elementVolume / 3;
                subContVol[1].volume_ = elementVolume / 3;
                subContVol[2].volume_ = elementVolume / 3;
                break;
            case 4: // 2D, quadrilinear
                subContVol[0].volume_ =
                    quadrilateralArea(subContVol[0].global,
                                      edgeCoord[2],
                                      elementGlobal,
                                      edgeCoord[0]);
                subContVol[1].volume_ =
                    quadrilateralArea(subContVol[1].global,
                                      edgeCoord[1],
                                      elementGlobal,
                                      edgeCoord[2]);
                subContVol[2].volume_ =
                    quadrilateralArea(subContVol[2].global,
                                      edgeCoord[0],
                                      elementGlobal,
                                      edgeCoord[3]);
                subContVol[3].volume_ =
                    quadrilateralArea(subContVol[3].global,
                                      edgeCoord[3],
                                      elementGlobal,
                                      edgeCoord[1]);
                break;
            default:
                throw std::logic_error("Not implemented:VcfvStencil dim = " + std::to_string(dim) +
                                       ", numVertices = " + std::to_string(numVertices));
            }
        }
        else if constexpr (dim == 3) {
            switch (numVertices) {
            case 4: // 3D, tetrahedron
                for (unsigned k = 0; k < numVertices; k++) {
                    subContVol[k].volume_ = elementVolume / 4.0;
                }
                break;
            case 5: // 3D, pyramid
                subContVol[0].volume_ =
                    hexahedronVolume(subContVol[0].global,
                                     edgeCoord[2],
                                     faceCoord[0],
                                     edgeCoord[0],
                                     edgeCoord[4],
                                     faceCoord[3],
                                     elementGlobal,
                                     faceCoord[1]);
                subContVol[1].volume_ =
                    hexahedronVolume(subContVol[1].global,
                                     edgeCoord[1],
                                     faceCoord[0],
                                     edgeCoord[2],
                                     edgeCoord[5],
                                     faceCoord[2],
                                     elementGlobal,
                                     faceCoord[3]);
                subContVol[2].volume_ =
                    hexahedronVolume(subContVol[2].global,
                                     edgeCoord[0],
                                     faceCoord[0],
                                     edgeCoord[3],
                                     edgeCoord[6],
                                     faceCoord[1],
                                     elementGlobal,
                                     faceCoord[4]);
                subContVol[3].volume_ =
                    hexahedronVolume(subContVol[3].global,
                                     edgeCoord[3],
                                     faceCoord[0],
                                     edgeCoord[1],
                                     edgeCoord[7],
                                     faceCoord[4],
                                     elementGlobal,
                                     faceCoord[2]);
                subContVol[4].volume_ = elementVolume -
                                        subContVol[0].volume_ -
                                        subContVol[1].volume_ -
                                        subContVol[2].volume_ -
                                        subContVol[3].volume_;
                break;
            case 6: // 3D, prism
                subContVol[0].volume_ =
                    hexahedronVolume(subContVol[0].global,
                                     edgeCoord[3],
                                     faceCoord[3],
                                     edgeCoord[4],
                                     edgeCoord[0],
                                     faceCoord[0],
                                     elementGlobal,
                                     faceCoord[1]);
                subContVol[1].volume_ =
                    hexahedronVolume(subContVol[1].global,
                                     edgeCoord[5],
                                     faceCoord[3],
                                     edgeCoord[3],
                                     edgeCoord[1],
                                     faceCoord[2],
                                     elementGlobal,
                                     faceCoord[0]);
                subContVol[2].volume_ =
                    hexahedronVolume(subContVol[2].global,
                                     edgeCoord[4],
                                     faceCoord[3],
                                     edgeCoord[5],
                                     edgeCoord[2],
                                     faceCoord[1],
                                     elementGlobal,
                                     faceCoord[2]);
                subContVol[3].volume_ =
                    hexahedronVolume(edgeCoord[0],
                                     faceCoord[0],
                                     elementGlobal,
                                     faceCoord[1],
                                     subContVol[3].global,
                                     edgeCoord[6],
                                     faceCoord[4],
                                     edgeCoord[7]);
                subContVol[4].volume_ =
                    hexahedronVolume(edgeCoord[1],
                                     faceCoord[2],
                                     elementGlobal,
                                     faceCoord[0],
                                     subContVol[4].global,
                                     edgeCoord[8],
                                     faceCoord[4],
                                     edgeCoord[6]);
                subContVol[5].volume_ =
                    hexahedronVolume(edgeCoord[2],
                                     faceCoord[1],
                                     elementGlobal,
                                     faceCoord[2],
                                     subContVol[5].global,
                                     edgeCoord[7],
                                     faceCoord[4],
                                     edgeCoord[8]);
                break;
            case 8: // 3D, hexahedron
                subContVol[0].volume_ =
                    hexahedronVolume(subContVol[0].global,
                                     edgeCoord[6],
                                     faceCoord[4],
                                     edgeCoord[4],
                                     edgeCoord[0],
                                     faceCoord[2],
                                     elementGlobal,
                                     faceCoord[0]);
                subContVol[1].volume_ =
                    hexahedronVolume(subContVol[1].global,
                                     edgeCoord[5],
                                     faceCoord[4],
                                     edgeCoord[6],
                                     edgeCoord[1],
                                     faceCoord[1],
                                     elementGlobal,
                                     faceCoord[2]);
                subContVol[2].volume_ =
                    hexahedronVolume(subContVol[2].global,
                                     edgeCoord[4],
                                     faceCoord[4],
                                     edgeCoord[7],
                                     edgeCoord[2],
                                     faceCoord[0],
                                     elementGlobal,
                                     faceCoord[3]);
                subContVol[3].volume_ =
                    hexahedronVolume(subContVol[3].global,
                                     edgeCoord[7],
                                     faceCoord[4],
                                     edgeCoord[5],
                                     edgeCoord[3],
                                     faceCoord[3],
                                     elementGlobal,
                                     faceCoord[1]);
                subContVol[4].volume_ =
                    hexahedronVolume(edgeCoord[0],
                                     faceCoord[2],
                                     elementGlobal,
                                     faceCoord[0],
                                     subContVol[4].global,
                                     edgeCoord[10],
                                     faceCoord[5],
                                     edgeCoord[8]);
                subContVol[5].volume_ =
                    hexahedronVolume(edgeCoord[1],
                                     faceCoord[1],
                                     elementGlobal,
                                     faceCoord[2],
                                     subContVol[5].global,
                                     edgeCoord[9],
                                     faceCoord[5],
                                     edgeCoord[10]);
                subContVol[6].volume_ =
                    hexahedronVolume(edgeCoord[2],
                                     faceCoord[0],
                                     elementGlobal,
                                     faceCoord[3],
                                     subContVol[6].global,
                                     edgeCoord[8],
                                     faceCoord[5],
                                     edgeCoord[11]);
                subContVol[7].volume_ =
                    hexahedronVolume(edgeCoord[3],
                                     faceCoord[3],
                                     elementGlobal,
                                     faceCoord[1],
                                     subContVol[7].global,
                                     edgeCoord[11],
                                     faceCoord[5],
                                     edgeCoord[9]);
                break;
            default:
                throw std::logic_error("Not implemented:VcfvStencil for dim = " + std::to_string(dim) +
                                       ", numVertices = " + std::to_string(numVertices));
            }
        }
        else {
            throw std::logic_error("Not implemented:VcfvStencil for dim = " + std::to_string(dim));
        }
    }
#if __GNUC__ || __clang__
#pragma GCC diagnostic pop
#endif

    const GridView&     gridView_;
    const Mapper& vertexMapper_;

    Element element_;

#if HAVE_DUNE_LOCALFUNCTIONS
    static LocalFiniteElementCache feCache_;
#endif // HAVE_DUNE_LOCALFUNCTIONS

    LocalPosition elementLocal;

    GlobalPosition elementGlobal;

    Scalar elementVolume;

    std::array<SubControlVolume, maxNC> subContVol{};

    std::array<SubControlVolumeFace, maxNE> subContVolFace{};

    std::array<BoundaryFace, maxBF> boundaryFace_{};

    unsigned numBoundarySegments_{};

    std::array<GlobalPosition, maxNE> edgeCoord{};

    std::array<GlobalPosition, maxNF> faceCoord{};

    unsigned numVertices{};

    unsigned numEdges{};

    unsigned numFaces{};

    Dune::GeometryType geometryType_;
};

#if HAVE_DUNE_LOCALFUNCTIONS
template<class Scalar, class GridView>
typename VcfvStencil<Scalar, GridView>::LocalFiniteElementCache
VcfvStencil<Scalar, GridView>::feCache_;
#endif // HAVE_DUNE_LOCALFUNCTIONS

template <class Scalar>
std::array<typename VcfvScvGeometries<Scalar, /*dim=*/1, ElementType::cube>::ScvLocalGeometry,
           VcfvScvGeometries<Scalar, /*dim=*/1, ElementType::cube>::numScv>
VcfvScvGeometries<Scalar, /*dim=*/1, ElementType::cube>::scvGeoms_{};

template <class Scalar>
std::array<typename VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::simplex>::ScvLocalGeometry,
           VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::simplex>::numScv>
VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::simplex>::scvGeoms_{};

template <class Scalar>
std::array<typename VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::cube>::ScvLocalGeometry,
           VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::cube>::numScv>
VcfvScvGeometries<Scalar, /*dim=*/2, ElementType::cube>::scvGeoms_{};

template <class Scalar>
std::array<typename VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::simplex>::ScvLocalGeometry,
           VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::simplex>::numScv>
VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::simplex>::scvGeoms_{};

template <class Scalar>
std::array<typename VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::cube>::ScvLocalGeometry,
           VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::cube>::numScv>
VcfvScvGeometries<Scalar, /*dim=*/3, ElementType::cube>::scvGeoms_{};

} // namespace Opm

#endif