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