nav.hpp

Go to the documentation of this file.

#ifndef OPM_VERTEQ_NAV_HPP_INCLUDED
#define OPM_VERTEQ_NAV_HPP_INCLUDED

// Copyright (C) 2013 Uni Research AS
// This file is licensed under the GNU General Public License v3.0

#ifndef OPM_GRID_HEADER_INCLUDED
#include <opm/core/grid.h>
#endif

#include <cstdlib>  // div_t
#include <iosfwd>   // ostream

struct Coord2D {
    Coord2D (int i_, int j_) : m_i (i_), m_j (j_) { }

    int i() const { return m_i; }
    int j() const { return m_j; }

    bool operator == (const Coord2D& rhs) const {
        return (m_i == rhs.m_i) && (m_j == rhs.m_j);
    }

protected:
    const int m_i;
    const int m_j;

    friend std::ostream& operator << (std::ostream& s, const Coord2D& c);
};

struct Coord3D : public Coord2D {
    Coord3D (int i_, int j_, int k_)
        : Coord2D (i_, j_)
        , m_k (k_) {
    }

    int k() const { return m_k; }

protected:
    const int m_k;

    friend std::ostream& operator << (std::ostream& s, const Coord3D& c);
};

// forward declaration
template <typename Dim > struct Side;

struct Dir {
    static const Dir DEC; // = 0

    static const Dir INC; // = 1

    static const int COUNT = 2;

    const int val;

    Dir (const Dir& rhs) : val (rhs.val) {}
    bool operator == (const Dir& rhs) const { return val == rhs.val; }

    Dir opposite () const { return Dir (-(val - 1)); }

protected:
    Dir (int i) : val (i) { }

    template <typename Dim> friend struct Side;

    friend std::ostream& operator << (std::ostream& os, const Dir& d);
};

struct Dim1D {
    static const int COUNT = 1;
};

struct Dim2D : public Dim1D {
    // two spatial directions
    static const Dim2D X; // = 0
    static const Dim2D Y; // = 1

    // number of dimensions
    static const int COUNT = 2;

    const int val;

    Dim2D (const Dim2D& rhs) : val (rhs.val) { }
    bool operator == (const Dim2D& rhs) const { return val == rhs.val; }

    Dim2D orthogonal () const { return Dim2D (-(val - 1)); }

protected:
    Dim2D (int i) : val (i) { }

    friend struct Side <Dim2D>;

    friend std::ostream& operator << (std::ostream& os, const Dim2D& d);
};

struct Dim3D : public Dim2D {
    // added dimension in 3D
    static const Dim3D Z; // = 2

    // number of dimensions (shadows Dim2D)
    static const int COUNT = 3;

    Dim3D (const Dim3D& rhs) : Dim2D (rhs.val) { }
    bool operator == (const Dim2D& rhs) const { return val == rhs.val; }

    // allow X and Y to work in 3D too
    Dim3D (const Dim2D& rhs) : Dim2D (rhs) {}

protected:
    Dim3D (int i) : Dim2D (i) { }

    friend struct Side <Dim3D>;
};

template <typename Dim>
struct Side {
    Side (Dim dim_, Dir dir_) : m_dim (dim_), m_dir (dir_) { }
    Side (const Side& rhs) : m_dim (rhs.m_dim), m_dir (rhs.m_dir) { }

    int facetag () const {
        return dim().val * Dir::COUNT + dir().val;
    }

    static Side <Dim> from_tag (int tag);

    Dim dim() const { return m_dim; }
    Dir dir() const { return m_dir; }

    static const int COUNT = Dim::COUNT * Dir::COUNT;

    static const Side* begin () { return &ALL[0]; }
    static const Side* end () { return &ALL[COUNT]; }

    bool operator == (const Side <Dim>& rhs) const {
        return (m_dim == rhs.m_dim) && (m_dir == rhs.m_dir);
    }

protected:
    const Dim m_dim;
    const Dir m_dir;

    // fixed enumeration of all sides
    static const Side ALL [];

    template <typename T>
    friend std::ostream& operator << (std::ostream& os, const Side<T>& s);
};

// specializations for the dimensions we work with
typedef Side <Dim2D> Side2D;
typedef Side <Dim3D> Side3D;

// forward declaration of stream operator present in library
std::ostream& operator << (std::ostream& os, const Side2D& s);
std::ostream& operator << (std::ostream& os, const Side3D& s);

// standalone constants for sides that we use; we call them 'up' and
// 'down' so that U and D are mnemonics, in contrast to 'top' and 'bottom'
// where the 'b' would conflict with 'back'.
extern const Side3D UP;
extern const Side3D DOWN;

struct Corn2D {
    Corn2D (Dir i_, Dir j_) : m_i (i_), m_j (j_) { }
    Corn2D (const Corn2D& rhs) : m_i (rhs.m_i), m_j (rhs.m_j) { }

    Dir i() const { return m_i; }
    Dir j() const { return m_j; }

protected:
    const Dir m_i;
    const Dir m_j;
};

struct Corn3D : public Corn2D {
    Corn3D (Dir i_, Dir j_, Dir k_) : Corn2D (i_, j_), m_k (k_) { }
    Corn3D (const Corn3D& rhs) : Corn2D (rhs), m_k (rhs.m_k) { }

    Corn3D pivot (Dim3D dim, Dir dir) {
        return Corn3D (dim == Dim3D::X ? dir : m_i,
                                     dim == Dim3D::Y ? dir : m_j,
                                     dim == Dim3D::Z ? dir : m_k);
    }

    Dir k() const { return m_k; }

    bool operator == (const Corn3D& rhs) const {
        return (m_i == rhs.m_i) && (m_j == rhs.m_j) && (m_k == rhs.m_k);
    }

protected:
    const Dir m_k;

    friend std::ostream& operator << (std::ostream& os, const Corn3D& c);
};

struct Cart2D {
    // number of cells in each direction
    const int ni;
    const int nj;

    // initialize from the size of the grid
    Cart2D (int ni_, int nj_) : ni (ni_), nj (nj_) { }

    // use these two value types for structured coordinates and
    // flattened indices, respectively
    typedef Coord2D coord_t;
    typedef int elem_t;
    typedef int node_t;
    typedef int face_t;

    static const int NO_ELEM; // = -1
    static const int NO_FACE; // = -1
    static const int NO_NODE; // = -1

    int num_elems () const {
        return ni * nj;
    }

    elem_t cart_ndx (const coord_t& coord) const {
        return coord.j() * ni + coord.i();
    }

    coord_t coord (const elem_t& cart_ndx) const {
        const std::div_t strip = std::div (cart_ndx, ni);
        const int i = strip.rem;
        const int j = strip.quot;
        return coord_t (i, j);
    }

    int num_nodes () const {
        return (ni + 1) * (nj + 1);
    }

    node_t node_ndx (const coord_t& coord, const Corn2D& corn) {
        return (coord.j() + corn.j().val) * (ni + 1) + (coord.i() + corn.i().val);
    }

    int num_faces () const {
        // there are ni+1 faces oriented in j-direction in each column,
        // for a total of (ni+1)*nj, and nj+1 faces in i-direction in
        // each row, for a total of (nj+1)*ni.
        return (ni + 1) * nj + ni * (nj + 1);
    }

    face_t face_ndx (const coord_t& coord, const Side2D& side) {
        // flags that code 1 (include) or 0 (don't) for each of the tests
        const int dirv = side.dir().val;
        const int idim = side.dim() == Dim2D::X ? 1 : 0;
        const int idir = idim ? dirv : 0;
        const int jdir = idim ? 0 : dirv;
        // there is a left side and a lower side face for each element in a
        // column, plus one face at the top of each column; 2*ni+1. the j-
        // coordinate plus one extra if we are selecting the right face, tells
        // us which column which is to the left of or "above" the face.
        // if the side is in the i-direction to the center of the element, then
        // the face is aligned with the j axis, so we skip idim*ni i-aligned
        // faces before it.
        // finally, determine the row by using the i-coordinate, and i-direction
        // to determine whether it is the lower or upper face (if applicable).
        return (coord.j() + jdir) * (2 * ni + 1) + (idim * ni) + (coord.i() + idir);
    }
};

struct Cart3D {
    // number of cells in each direction
    const int ni;
    const int nj;
    const int nk;

    Cart3D (const UnstructuredGrid& g)
        : ni (g.cartdims [0])
        , nj (g.cartdims [1])
        , nk (g.cartdims [2]) { }

    Cart2D project () const {
        return Cart2D (ni, nj);
    }

    // use these two value types for structured coordinates and
    // flattened indices, respectively
    typedef Coord3D coord_t;
    typedef int elem_t;

    coord_t coord (const elem_t& cart_ndx) const {
        // the i-index moves fastest, as this is Fortran-indexing
        const div_t strip = div (cart_ndx, ni);
        const int i = strip.rem;
        const div_t plane = div (strip.quot, nj);
        const int j = plane.rem;
        const int k = plane.quot;
        return coord_t (i, j, k);
    }
};

#endif // OPM_VERTEQ_NAV_HPP_INCLUDED