volcartesian.cpp

/*---------------------------------------------------------------------------*\
 *
 *  bitpit
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of bitpit.
 *
 *  bitpit is free software: you can redistribute it and/or modify it
 *  under the terms of the GNU Lesser General Public License v3 (LGPL)
 *  as published by the Free Software Foundation.
 *
 *  bitpit 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 Lesser General Public
 *  License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with bitpit. If not, see <http://www.gnu.org/licenses/>.
 *
\*---------------------------------------------------------------------------*/
 
#include <cmath>
#include <bitset>
#if BITPIT_ENABLE_MPI==1
#   include <mpi.h>
#endif
 
#include "bitpit_common.hpp"
 
#include "volcartesian.hpp"
 
namespace bitpit {
 
VolCartesian::VolCartesian()
#if BITPIT_ENABLE_MPI==1
    : VolumeKernel(MPI_COMM_NULL, 0, ADAPTION_DISABLED, PARTITIONING_DISABLED)
#else
    : VolumeKernel(ADAPTION_DISABLED)
#endif
{
    initialize();
}
 
VolCartesian::VolCartesian(int dimension,
                               const std::array<double, 3> &origin,
                               const std::array<double, 3> &lengths,
                               const std::array<int, 3> &nCells)
    : VolCartesian(PatchManager::AUTOMATIC_ID, dimension, origin, lengths, nCells)
{
}
 
VolCartesian::VolCartesian(int id, int dimension,
                               const std::array<double, 3> &origin,
                               const std::array<double, 3> &lengths,
                               const std::array<int, 3> &nCells)
#if BITPIT_ENABLE_MPI==1
    : VolumeKernel(id, dimension, MPI_COMM_NULL, 0, ADAPTION_DISABLED, PARTITIONING_DISABLED)
#else
    : VolumeKernel(id, dimension, ADAPTION_DISABLED)
#endif
{
    initialize();
 
    setOrigin(origin);
    setLengths(lengths);
 
    setDiscretization(nCells);
}
 
VolCartesian::VolCartesian(int dimension,
                               const std::array<double, 3> &origin,
                               double length, int nCells)
    : VolCartesian(PatchManager::AUTOMATIC_ID, dimension, origin, length, nCells)
{
}
 
VolCartesian::VolCartesian(int id, int dimension,
                               const std::array<double, 3> &origin,
                               double length, int nCells)
#if BITPIT_ENABLE_MPI==1
    : VolumeKernel(id, dimension, MPI_COMM_NULL, 0, ADAPTION_DISABLED, PARTITIONING_DISABLED)
#else
    : VolumeKernel(id, dimension, ADAPTION_DISABLED)
#endif
{
    initialize();
 
    setOrigin(origin);
    setLengths({{length, length, length}});
 
    setDiscretization({{nCells, nCells, nCells}});
}
 
VolCartesian::VolCartesian(int dimension,
                               const std::array<double, 3> &origin,
                               double length, double dh)
    : VolCartesian(PatchManager::AUTOMATIC_ID, dimension, origin, length, dh)
{
}
 
VolCartesian::VolCartesian(int id, int dimension,
                               const std::array<double, 3> &origin,
                               double length, double dh)
#if BITPIT_ENABLE_MPI==1
    : VolumeKernel(id, dimension, MPI_COMM_NULL, 0, ADAPTION_DISABLED, PARTITIONING_DISABLED)
#else
    : VolumeKernel(id, dimension, ADAPTION_DISABLED)
#endif
{
    initialize();
 
    setOrigin(origin);
    setLengths({{length, length, length}});
 
    int nCells = (int) std::ceil(length / dh);
    setDiscretization({{nCells, nCells, nCells}});
}
 
VolCartesian::VolCartesian(std::istream &stream)
#if BITPIT_ENABLE_MPI==1
    : VolumeKernel(MPI_COMM_NULL, 0, ADAPTION_DISABLED, PARTITIONING_DISABLED)
#else
    : VolumeKernel(ADAPTION_DISABLED)
#endif
{
    initialize();
    restore(stream);
}
 
std::unique_ptr<PatchKernel> VolCartesian::clone()  const
{
    return std::unique_ptr<VolCartesian>(new VolCartesian(*this));
}
 
void VolCartesian::reset()
{
    // Switch to light memeory mode
    switchMemoryMode(MEMORY_LIGHT);
 
    // Reset geometry and discretization
    m_minCoords = {{0., 0., 0.}};
    m_maxCoords = {{1., 1., 1.}};
 
    m_nCells1D     = {{0, 0, 0}};
    m_nVertices1D  = {{0, 0, 0}};
    m_cellSpacings = {{0., 0., 0.}};
 
    for (int n = 0; n < 3; ++n) {
        std::vector<double>().swap(m_vertexCoords[n]);
        std::vector<double>().swap(m_cellCenters[n]);
    }
 
    m_nVertices   = 0;
    m_nCells      = 0;
    m_nInterfaces = 0;
}
 
void VolCartesian::resetInterfaces()
{
    PatchKernel::resetInterfaces();
}
 
 
void VolCartesian::_updateAdjacencies()
{
    // Partial updates are not supported
    CellConstIterator beginItr = cellConstBegin();
    CellConstIterator endItr   = cellConstEnd();
 
    bool partialUpdate = false;
    for (CellConstIterator cellIterator = beginItr; cellIterator != endItr; ++cellIterator) {
        long cellId = cellIterator.getId();
        if (!testCellAlterationFlags(cellId, FLAG_ADJACENCIES_DIRTY)) {
            partialUpdate = true;
            break;
        }
    }
 
    if (partialUpdate) {
        log::cout() << " It is not possible to partially update the adjacencies.";
        log::cout() << " All adjacencies will be updated.";
    }
 
    // Update adjacencies
    int nCellFaces = 2 * getDimension();
    for (Cell &cell : getCells()) {
        // Reset cell adjacencies
        if (cell.getAdjacencyCount() != 0) {
            cell.resetAdjacencies();
        }
 
        // Evaluate adjacencies
        long cellId = cell.getId();
        for (int face = 0; face < nCellFaces; ++face) {
            // Identify face neighbour
            long neighId = getCellFaceNeighsLinearId(cellId, face);
            if (neighId < 0) {
                continue;
            }
 
            // Add adjacency
            cell.pushAdjacency(face, neighId);
        }
    }
}
 
void VolCartesian::_updateInterfaces()
{
    // Partial updates are not supported
    CellConstIterator beginItr = cellConstBegin();
    CellConstIterator endItr   = cellConstEnd();
 
    bool partialUpdate = false;
    for (CellConstIterator cellIterator = beginItr; cellIterator != endItr; ++cellIterator) {
        long cellId = cellIterator.getId();
        if (!testCellAlterationFlags(cellId, FLAG_INTERFACES_DIRTY)) {
            partialUpdate = true;
            break;
        }
    }
 
    if (partialUpdate) {
        log::cout() << " It is not possible to partially update the interfaces.";
        log::cout() << " All interfaces will be updated.";
    }
 
    // Build interfaces
    if (getMemoryMode() == MemoryMode::MEMORY_NORMAL) {
        int nCellFaces = 2 * getDimension();
 
        // Info on the interfaces
        ElementType interfaceType = getInterfaceType();
 
        const ReferenceElementInfo &interfaceTypeInfo = ReferenceElementInfo::getInfo(interfaceType);
        const int nInterfaceVertices = interfaceTypeInfo.nVertices;
 
        // Enable manual adaption
        AdaptionMode previousAdaptionMode = getAdaptionMode();
        setAdaptionMode(ADAPTION_MANUAL);
 
        // Initialize interfaces
        for (Cell &cell : getCells()) {
            cell.setInterfaces(FlatVector2D<long>(nCellFaces, 1, Interface::NULL_ID));
        }
 
        // Build interfaces
        for (Cell &cell : getCells()) {
            long cellId = cell.getId();
            for (int face = 0; face < nCellFaces; ++face) {
                // Get neighbour information
                //
                // The interface between two cells needs to be processed only
                // once. In order to ensure this, we build an interface if the
                // faces is a border, or if the id of this cell is lower than
                // the id of its neighbour.
                long neighId;
                if (cell.getAdjacencyCount(face) > 0) {
                    neighId = cell.getAdjacencies(face)[0];
                    if (neighId < cellId) {
                        continue;
                    }
                } else {
                    neighId = Cell::NULL_ID;
                }
 
                // Create the interface
                InterfaceIterator interfaceIterator = VolumeKernel::addInterface(interfaceType);
                Interface &interface = *interfaceIterator;
                long interfaceId = interface.getId();
 
                // Set connectivity
                ConstProxyVector<long> faceConnect = cell.getFaceConnect(face);
                int connectSize = faceConnect.size();
 
                std::unique_ptr<long[]> connect = std::unique_ptr<long[]>(new long[nInterfaceVertices]);
                for (int k = 0; k < connectSize; ++k) {
                    connect[k] = faceConnect[k];
                }
                interface.setConnect(std::move(connect));
 
                // Set owner data
                interface.setOwner(cellId, face);
                cell.setInterface(face, 0, interfaceId);
 
                // Neighbour data
                if (neighId >= 0) {
                    Cell &neigh = getCell(neighId);
                    int neighFace = 2 * static_cast<int>(std::floor(face / 2.)) + (1 - face % 2);
 
                    interface.setNeigh(neighId, neighFace);
                    neigh.setInterface(neighFace, 0, interfaceId);
                } else {
                    interface.unsetNeigh();
                }
            }
        }
 
        // Restore previous adaption mode
        setAdaptionMode(previousAdaptionMode);
    }
}
 
void VolCartesian::initialize()
{
    // Normals
    int i = 0;
    for (int n = 0; n < 3; n++) {
        for (int k = -1; k <= 1; k += 2) {
            std::array<double, 3> normal = {{0.0, 0.0, 0.0}};
            normal[n] = k;
 
            m_normals[i++] = normal;
        }
    }
 
    // Deltas for the evaluation of the vertex neighbours
    m_vertexNeighDeltas = std::vector<std::array<int, 3>>(8);
    m_vertexNeighDeltas[0] = {{ 0,  0, 0}};
    m_vertexNeighDeltas[1] = {{-1,  0, 0}};
    m_vertexNeighDeltas[2] = {{ 0, -1, 0}};
    m_vertexNeighDeltas[3] = {{-1, -1, 0}};
    m_vertexNeighDeltas[4] = {{ 0,  0, -1}};
    m_vertexNeighDeltas[5] = {{-1,  0, -1}};
    m_vertexNeighDeltas[6] = {{ 0, -1, -1}};
    m_vertexNeighDeltas[7] = {{-1, -1, -1}};
 
    // Deltas for the evaluation of the edge neighbours
    m_edgeNeighDeltas = std::vector<std::array<int, 3>>(12);
    m_edgeNeighDeltas[ 0] = {{-1,  0, -1}};
    m_edgeNeighDeltas[ 1] = {{ 1,  0, -1}};
    m_edgeNeighDeltas[ 2] = {{ 0, -1, -1}};
    m_edgeNeighDeltas[ 3] = {{ 0,  1, -1}};
    m_edgeNeighDeltas[ 4] = {{-1, -1,  0}};
    m_edgeNeighDeltas[ 5] = {{ 1, -1,  0}};
    m_edgeNeighDeltas[ 6] = {{-1,  1,  0}};
    m_edgeNeighDeltas[ 7] = {{ 1,  1,  0}};
    m_edgeNeighDeltas[ 8] = {{-1,  0,  1}};
    m_edgeNeighDeltas[ 9] = {{ 1,  0,  1}};
    m_edgeNeighDeltas[10] = {{ 0, -1,  1}};
    m_edgeNeighDeltas[11] = {{ 0,  1,  1}};
 
    // Faces associated to the edges
    m_edgeFaces = std::vector<std::array<int, 2>>(12);
    m_edgeFaces[ 0] = {{ 0, 4}};
    m_edgeFaces[ 1] = {{ 1, 4}};
    m_edgeFaces[ 2] = {{ 2, 4}};
    m_edgeFaces[ 3] = {{ 3, 4}};
    m_edgeFaces[ 4] = {{ 0, 2}};
    m_edgeFaces[ 5] = {{ 1, 2}};
    m_edgeFaces[ 6] = {{ 0, 3}};
    m_edgeFaces[ 7] = {{ 1, 3}};
    m_edgeFaces[ 8] = {{ 0, 5}};
    m_edgeFaces[ 9] = {{ 1, 5}};
    m_edgeFaces[10] = {{ 2, 5}};
    m_edgeFaces[11] = {{ 3, 5}};
 
    // Set the bounding box as frozen
    setBoundingBoxFrozen(true);
 
    // Set the light memory mode
    setMemoryMode(MemoryMode::MEMORY_LIGHT);
 
    // Reset the patch
    reset();
}
 
void VolCartesian::initializeCellVolume()
{
    m_cellVolume = 1;
    for (int d = 0; d < getDimension(); ++d) {
        if (m_directionOrdering[d] >= 0) {
            m_cellVolume *= m_cellSpacings[d];
        }
    }
}
 
void VolCartesian::initializeCellSize()
{
    m_cellSize = std::pow(m_cellVolume, 1. / getDimension());
}
 
void VolCartesian::initializeInterfaceArea()
{
    for (int d = 0; d < getDimension(); ++d) {
        if (m_directionOrdering[d] >= 0) {
            m_interfaceArea[d] = m_cellVolume / m_cellSpacings[d];
        } else {
            m_interfaceArea[d] = 0.;
        }
    }
}
 
void VolCartesian::setDiscretization(const std::array<int, 3> &nCells)
{
    // Spacing
    for (int d = 0; d < getDimension(); ++d) {
        // Initialize cells
        if (nCells[d] >  0) {
            m_nCells1D[d]     = nCells[d];
            m_cellSpacings[d] = (m_maxCoords[d] - m_minCoords[d]) / m_nCells1D[d];
 
            m_cellCenters[d].resize(m_nCells1D[d]);
            for (int i = 0; i < m_nCells1D[d]; i++) {
                m_cellCenters[d][i] = m_minCoords[d] + (0.5 + i) * m_cellSpacings[d];
            }
        } else {
            m_nCells1D[d]     = 0;
            m_cellSpacings[d] = 0.;
        }
 
        log::cout() << "  - Cell count along direction " << d << " : " << m_nCells1D[d] << "\n";
 
        // Initialize vertices
        m_nVertices1D[d] = m_nCells1D[d] + 1;
 
        m_vertexCoords[d].resize(m_nVertices1D[d]);
        for (int i = 0; i < m_nVertices1D[d]; i++) {
            m_vertexCoords[d][i] = m_minCoords[d] + i * m_cellSpacings[d];
        }
    }
 
    for (int d = getDimension(); d < 3; ++d) {
        m_nCells1D[d]     = 0;
        m_cellSpacings[d] = 0.;
 
        m_nVertices1D[d] = 1;
    }
 
    log::cout() << std::endl;
 
    // Define direction ordering
    //
    // It is the ordering in which the different directions will be considered
    // when numbering cell, vertices, ... For each direction, it will contain
    // the order in which the direction should be considered or -1 if there
    // are no cells along the direction.
    for (int d = 0; d < 3; ++d) {
        if (d == 0) {
            m_directionOrdering[d] = 0;
        } else {
            int previousDirection = m_directionOrdering[d - 1];
            if (previousDirection < 0 || previousDirection >= 2) {
                m_directionOrdering[d] = -1;
                continue;
            }
 
            m_directionOrdering[d] = previousDirection  + 1;
        }
 
        while (m_nCells1D[m_directionOrdering[d]] == 0) {
            if (m_directionOrdering[d] == 2) {
                m_directionOrdering[d] = -1;
                break;
            }
 
            ++m_directionOrdering[d];
        }
    }
 
    // Count the total number of vertices
    m_nVertices = 1;
    for (int n = 0; n < getDimension(); ++n) {
        m_nVertices *= m_nVertices1D[n];
    }
    log::cout() << "  - Total vertex count: " << m_nVertices << "\n";
 
    // Count the total number of cells
    m_nCells = 1;
    for (int d = 0; d < getDimension(); ++d) {
        if (m_directionOrdering[d] >= 0) {
            m_nCells *= m_nCells1D[d];
        }
    }
 
    log::cout() << "  - Total cell count: " << m_nCells << "\n";
 
    // Count the total number of interfaces
    m_nInterfaces = 0;
    for (int d = 0; d < getDimension(); ++d) {
        int nDirectionInterfaces = 1;
        for (int n = 0; n < getDimension(); n++) {
            int nDirectionInterfaces1D = m_nCells1D[n];
            if (n == d) {
                ++nDirectionInterfaces1D;
            }
 
            nDirectionInterfaces *= nDirectionInterfaces1D;
        }
        m_nInterfaces += nDirectionInterfaces;
    }
 
    log::cout() << "  - Total interface count: " << m_nInterfaces << "\n";
 
    // Cell volume
    initializeCellVolume();
 
    // Cell size
    initializeCellSize();
 
    // Interface area
    initializeInterfaceArea();
 
    // Update patch data structures
    if (getMemoryMode() == MemoryMode::MEMORY_NORMAL) {
        // Switch to light mode to reset patchdata structures
        switchMemoryMode(MemoryMode::MEMORY_LIGHT);
 
        // Swich back to normal mode to rebuild patchdata structures
        switchMemoryMode(MemoryMode::MEMORY_NORMAL);
    }
}
 
long VolCartesian::getVertexCount() const
{
    return m_nVertices;
}
 
int VolCartesian::getVertexCount(int direction) const
{
    return m_nVertices1D[direction];
}
 
long VolCartesian::getCellCount() const
{
    return m_nCells;
}
 
int VolCartesian::getCellCount(int direction) const
{
    return m_nCells1D[direction];
}
 
ElementType VolCartesian::getCellType(long id) const
{
    BITPIT_UNUSED(id);
 
    return getCellType();
}
 
ElementType VolCartesian::getCellType() const
{
    if (isThreeDimensional()) {
        return ElementType::VOXEL;
    } else {
        return ElementType::PIXEL;
    }
}
 
long VolCartesian::getInterfaceCount() const
{
    return m_nInterfaces;
}
 
ElementType VolCartesian::getInterfaceType(long id) const
{
    BITPIT_UNUSED(id);
 
    return getInterfaceType();
}
 
ElementType VolCartesian::getInterfaceType() const
{
    if (isThreeDimensional()) {
        return ElementType::PIXEL;
    } else {
        return ElementType::LINE;
    }
}
 
std::array<double, 3> VolCartesian::evalVertexCoords(long id) const
{
    std::array<int, 3> ijk = getVertexCartesianId(id);
 
    return evalVertexCoords(ijk);
}
 
std::array<double, 3> VolCartesian::evalVertexCoords(const std::array<int, 3> &ijk) const
{
    std::array<double, 3> coords;
    for (int d = 0; d < 3; ++d) {
        if (ijk[d] >= 0) {
            coords[d] = m_vertexCoords[d][ijk[d]];
        } else {
            coords[d] = m_minCoords[d];
        }
    }
 
    return coords;
}
 
const std::vector<double> & VolCartesian::getVertexCoords(int direction) const
{
    return m_vertexCoords[direction];
}
 
double VolCartesian::evalCellVolume(long id) const
{
    BITPIT_UNUSED(id);
 
    return evalCellVolume();
}
 
double VolCartesian::evalCellVolume(const std::array<int, 3> &ijk) const
{
    BITPIT_UNUSED(ijk);
 
    return evalCellVolume();
}
 
double VolCartesian::evalCellVolume() const
{
    return m_cellVolume;
}
 
double VolCartesian::evalCellSize(long id) const
{
    BITPIT_UNUSED(id);
 
    return evalCellSize();
}
 
double VolCartesian::evalCellSize(const std::array<int, 3> &ijk) const
{
    BITPIT_UNUSED(ijk);
 
    return evalCellSize();
}
 
double VolCartesian::evalCellSize() const
{
    return m_cellSize;
}
 
void VolCartesian::evalCellBoundingBox(long id, std::array<double,3> *minPoint, std::array<double,3> *maxPoint) const
{
    std::array<double,3> cellCentroid = evalCellCentroid(id);
    for (int d = 0; d < 3; ++d) {
        (*minPoint)[d] = cellCentroid[d] - 0.5 * m_cellSpacings[d];
        (*maxPoint)[d] = cellCentroid[d] + 0.5 * m_cellSpacings[d];
    }
}
 
double VolCartesian::evalInterfaceArea(long id) const
{
    const Interface &interface = getInterface(id);
    int ownerFace = interface.getOwnerFace();
    int direction = static_cast<int>(std::floor(ownerFace / 2.));
 
    return m_interfaceArea[direction];
}
 
std::array<double, 3> VolCartesian::evalInterfaceNormal(long id) const
{
    const Interface &interface = getInterface(id);
    int ownerFace = interface.getOwnerFace();
 
    return m_normals[ownerFace];
}
 
std::array<double, 3> VolCartesian::getSpacing() const
{
    return m_cellSpacings;
}
 
void VolCartesian::switchMemoryMode(MemoryMode mode)
{
    // Early return if the current memory mode matches the requested one
    if (mode == getMemoryMode()) {
        return;
    }
 
    // Update patch data structures
    switch (mode) {
 
    case MemoryMode::MEMORY_NORMAL:
    {
        // Enable manual adaption
        AdaptionMode previousAdaptionMode = getAdaptionMode();
        setAdaptionMode(ADAPTION_MANUAL);
 
        // Create the mesh
        addVertices();
        addCells();
 
        // Restore previous adaption mode
        setAdaptionMode(previousAdaptionMode);
 
        break;
    }
 
    case MemoryMode::MEMORY_LIGHT:
    {
        // To put the patch in memory mode we need to reset the generic data
        // of the patch, therefore we can call the 'reset' implementation of
        // the kernel.
        VolumeKernel::reset();
 
        break;
    }
 
    }
 
    // Set the requested memory mode
    setMemoryMode(mode);
}
 
void VolCartesian::setMemoryMode(MemoryMode mode)
{
    m_memoryMode = mode;
}
 
VolCartesian::MemoryMode VolCartesian::getMemoryMode() const
{
    return m_memoryMode;
}
 
double VolCartesian::getSpacing(int direction) const
{
    return m_cellSpacings[direction];
}
 
void VolCartesian::addVertices()
{
    log::cout() << "  >> Creating vertices\n";
 
    log::cout() << "    - Vertex count: " << m_nVertices << "\n";
 
    m_vertices.reserve(m_nVertices);
    for (int k = 0; (isThreeDimensional()) ? (k < m_nVertices1D[Vertex::COORD_Z]) : (k <= 0); k++) {
        for (int j = 0; j < m_nVertices1D[Vertex::COORD_Y]; j++) {
            for (int i = 0; i < m_nVertices1D[Vertex::COORD_X]; i++) {
                // Linear index of the vertex
                long id_vertex = getVertexLinearId(i, j, k);
 
                // Vertex coordinates
                std::array<double, 3> coords = evalVertexCoords(id_vertex);
 
                // Add vertex
                VolumeKernel::addVertex(std::move(coords), id_vertex);
            }
        }
    }
}
 
void VolCartesian::addCells()
{
    log::cout() << "  >> Creating cells\n";
 
    // Info on the cells
    ElementType cellType = getCellType();
 
    // Create the cells
    log::cout() << "    - Cell count: " << m_nCells << "\n";
 
    m_cells.reserve(m_nCells);
    for (int k = 0; (isThreeDimensional()) ? (k < m_nCells1D[Vertex::COORD_Z]) : (k <= 0); k++) {
        for (int j = 0; j < m_nCells1D[Vertex::COORD_Y]; j++) {
            for (int i = 0; i < m_nCells1D[Vertex::COORD_X]; i++) {
                long id_cell = getCellLinearId(i, j, k);
                CellIterator cellIterator = VolumeKernel::addCell(cellType, id_cell);
                Cell &cell = *cellIterator;
 
                // Connettività
                long *cellConnect = cell.getConnect();
 
                cellConnect[0] = getVertexLinearId(i,     j,     k);
                cellConnect[1] = getVertexLinearId(i + 1, j,     k);
                cellConnect[2] = getVertexLinearId(i,     j + 1, k);
                cellConnect[3] = getVertexLinearId(i + 1, j + 1, k);
                if (cellType == ElementType::VOXEL) {
                    cellConnect[4] = getVertexLinearId(i,     j,     k + 1);
                    cellConnect[5] = getVertexLinearId(i + 1, j,     k + 1);
                    cellConnect[6] = getVertexLinearId(i,     j + 1, k + 1);
                    cellConnect[7] = getVertexLinearId(i + 1, j + 1, k + 1);
                }
            }
        }
    }
}
 
int VolCartesian::_getDumpVersion() const
{
    const int DUMP_VERSION = 2;
 
    return DUMP_VERSION;
}
 
void VolCartesian::_dump(std::ostream &stream) const
{
    utils::binary::write(stream, m_minCoords[0]);
    utils::binary::write(stream, m_minCoords[1]);
    utils::binary::write(stream, m_minCoords[2]);
 
    utils::binary::write(stream, m_maxCoords[0] - m_minCoords[0]);
    utils::binary::write(stream, m_maxCoords[1] - m_minCoords[1]);
    utils::binary::write(stream, m_maxCoords[2] - m_minCoords[2]);
 
    utils::binary::write(stream, m_nCells1D[0]);
    utils::binary::write(stream, m_nCells1D[1]);
    utils::binary::write(stream, m_nCells1D[2]);
 
    utils::binary::write(stream, m_memoryMode);
 
    if (m_memoryMode == MEMORY_NORMAL) {
        for (const Cell &cell : getCells()) {
            utils::binary::write(stream, cell.getPID());
        }
    }
 
    if (m_memoryMode == MEMORY_NORMAL) {
        for (const Interface &interface : getInterfaces()) {
            utils::binary::write(stream, interface.getPID());
        }
    }
}
 
void VolCartesian::_restore(std::istream &stream)
{
    // Origin
    std::array<double, 3> origin;
    utils::binary::read(stream, origin[0]);
    utils::binary::read(stream, origin[1]);
    utils::binary::read(stream, origin[2]);
 
    setOrigin(origin);
 
    // Lengths
    std::array<double, 3> lengths;
    utils::binary::read(stream, lengths[0]);
    utils::binary::read(stream, lengths[1]);
    utils::binary::read(stream, lengths[2]);
 
    setLengths(lengths);
 
    // Discretization
    std::array<int, 3> nCells;
    utils::binary::read(stream, nCells[0]);
    utils::binary::read(stream, nCells[1]);
    utils::binary::read(stream, nCells[2]);
 
    setDiscretization(nCells);
 
    // Memory mode
    MemoryMode memoryMode;
    utils::binary::read(stream, memoryMode);
    switchMemoryMode(memoryMode);
 
    // Cell data
    if (m_memoryMode == MEMORY_NORMAL) {
        for (Cell &cell : getCells()) {
            int PID;
            utils::binary::read(stream, PID);
            cell.setPID(PID);
        }
    }
 
    // Interface data
    if (m_memoryMode == MEMORY_NORMAL) {
        for (Interface &interface : getInterfaces()) {
            int PID;
            utils::binary::read(stream, PID);
            interface.setPID(PID);
        }
    }
}
 
bool VolCartesian::isPointInside(long id, const std::array<double, 3> &point) const
{
    std::array<int, 3> cellIjk = getCellCartesianId(id);
 
    const double EPS = getTol();
    for (int d = 0; d < 3; ++d){
        long index = cellIjk[d];
 
        double spacing = m_cellSpacings[d];
        double cellMin = m_cellCenters[d][index] - 0.5 * spacing;
        double cellMax = m_cellCenters[d][index] + 0.5 * spacing;
 
        if (point[d]< cellMin - EPS || point[d] > cellMax + EPS) {
            return false;
        }
    }
 
    return true;
}
 
bool VolCartesian::isPointInside(const std::array<double, 3> &point) const
{
    const double EPS = getTol();
 
    for (int n = 0; n < getDimension(); n++) {
        if (point[n] < m_minCoords[n] - EPS || point[n] > m_maxCoords[n] + EPS) {
            return false;
        }
    }
 
    return true;
}
 
long VolCartesian::locatePoint(const std::array<double, 3> &point) const
{
    std::array<int, 3> pointIjk = locatePointCartesian(point);
    if (isCellCartesianIdValid(pointIjk)) {
        return getCellLinearId(pointIjk);
    } else {
        return Element::NULL_ID;
    }
}
 
std::array<int, 3> VolCartesian::locatePointCartesian(const std::array<double, 3> &point) const
{
    std::array<int, 3> ijk;
    if (!isPointInside(point)) {
        ijk[0] = -1;
        ijk[1] = -1;
        ijk[2] = -1;
 
        return ijk;
    }
 
    ijk[0] = static_cast<int>(std::floor((point[Vertex::COORD_X] - m_minCoords[Vertex::COORD_X]) / m_cellSpacings[Vertex::COORD_X]));
    ijk[1] = static_cast<int>(std::floor((point[Vertex::COORD_Y] - m_minCoords[Vertex::COORD_Y]) / m_cellSpacings[Vertex::COORD_Y]));
    if (isThreeDimensional()) {
        ijk[2] = static_cast<int>(std::floor((point[Vertex::COORD_Z] - m_minCoords[Vertex::COORD_Z]) / m_cellSpacings[Vertex::COORD_Z]));
    } else {
        ijk[2] = -1;
    }
 
    return ijk;
}
 
long VolCartesian::locateClosestVertex(std::array<double,3> const &point) const
{
    return getVertexLinearId(locateClosestVertexCartesian(point));
}
 
std::array<int, 3> VolCartesian::locateClosestVertexCartesian(std::array<double,3> const &point) const
{
    std::array<int, 3> ijk;
    for (int n = 0; n < getDimension(); ++n) {
        ijk[n] = std::min(m_nVertices1D[n] - 1, std::max(0, (int) round((point[n] - m_minCoords[n]) / m_cellSpacings[n])));
    }
 
    if (!isThreeDimensional()) {
        ijk[Vertex::COORD_Z] = -1;
    }
 
    return ijk;
}
 
long VolCartesian::locateClosestCell(const std::array<double, 3> &point) const
{
    std::array<int, 3> pointIjk = locateClosestCellCartesian(point);
    return getCellLinearId(pointIjk);
}
 
std::array<int, 3> VolCartesian::locateClosestCellCartesian(const std::array<double, 3> &point) const
{
    std::array<int,3> ijk({{0,0,0}});
 
    for( int i=0; i<getDimension(); ++i){
        ijk[i] = static_cast<int>(std::floor((point[i] - m_minCoords[i]) / m_cellSpacings[i]));
        ijk[i] = std::max( ijk[i], 0 );
        ijk[i] = std::min( ijk[i], m_nCells1D[i]-1 );
    }
 
    return ijk;
}
 
long VolCartesian::getCellLinearId(int i, int j, int k) const
{
    return getCellLinearId(std::array<int, 3>{{i, j, k}});
}
 
long VolCartesian::getCellLinearId(const std::array<int, 3> &ijk) const
{
    int l = m_directionOrdering[0];
    int m = m_directionOrdering[1];
    int n = m_directionOrdering[2];
 
    if (l == -1) {
        return Cell::NULL_ID;
    }
 
    long id = ijk[l];
    if (m >= 0) {
        id += m_nCells1D[l] * ijk[m];
        if (n >= 0) {
            id += m_nCells1D[l] * m_nCells1D[m] * ijk[n];
        }
    }
 
    return id;
}
 
std::array<int, 3> VolCartesian::getCellCartesianId(long id) const
{
    int l = m_directionOrdering[0];
    int m = m_directionOrdering[1];
    int n = m_directionOrdering[2];
 
    std::array<int, 3> ijk = {{-1, -1, -1}};
    if (m == -1) {
        ijk[l] = id;
    } else if (l >= 0) {
        int offset_lm = m_nCells1D[l] * m_nCells1D[m];
        int index_n   = id / offset_lm;
 
        ijk[m] = (id - index_n * offset_lm) / m_nCells1D[l];
        ijk[l] = id - (id / m_nCells1D[l]) * m_nCells1D[l];
        if (n != -1) {
            ijk[n] = index_n;
        }
    }
 
    return ijk;
}
 
bool VolCartesian::isCellCartesianIdValid(const std::array<int, 3> &ijk) const
{
    for (int k = 0; k < getDimension(); ++k) {
        if (ijk[k] < 0 || ijk[k] >= m_nCells1D[k]) {
            return false;
        }
    }
 
    return true;
}
 
long VolCartesian::getVertexLinearId(int i, int j, int k) const
{
    return getVertexLinearId(std::array<int, 3>{{i, j, k}});
}
 
long VolCartesian::getVertexLinearId(const std::array<int, 3> &ijk) const
{
    int l = m_directionOrdering[0];
    int m = m_directionOrdering[1];
    int n = m_directionOrdering[2];
 
    if (l == -1) {
        return Vertex::NULL_ID;
    }
 
    long id = ijk[l];
    if (m >= 0) {
        id += m_nVertices1D[l] * ijk[m];
        if (n >= 0) {
            id += m_nVertices1D[l] * m_nVertices1D[m] * ijk[n];
        }
    }
 
    return id;
}
 
std::array<int, 3> VolCartesian::getVertexCartesianId(long id) const
{
    int l = m_directionOrdering[0];
    int m = m_directionOrdering[1];
    int n = m_directionOrdering[2];
 
    std::array<int, 3> ijk = {{-1, -1, -1}};
    if (m == -1) {
        ijk[l] = id;
    } else if (l >= 0) {
        int offset_lm = m_nVertices1D[l] * m_nVertices1D[m];
        int index_n   = id / offset_lm;
 
        ijk[m] = (id - index_n * offset_lm) / m_nVertices1D[l];
        ijk[l] = id - (id / m_nVertices1D[l]) * m_nVertices1D[l];
        if (n != -1) {
            ijk[n] = index_n;
        }
    }
 
    return ijk;
}
 
std::array<int, 3> VolCartesian::getVertexCartesianId(long cellId, int vertex) const
{
    return getVertexCartesianId(getCellCartesianId(cellId), vertex);
}
 
std::array<int, 3> VolCartesian::getVertexCartesianId(const std::array<int, 3> &cellIjk, int vertex) const
{
    std::bitset<3> vertexBitset(vertex);
    std::array<int, 3> vertexIjk(cellIjk);
    for (int k = 0; k < 3; ++k) {
        vertexIjk[k] += vertexBitset[k];
    }
 
    return vertexIjk;
}
 
bool VolCartesian::isVertexCartesianIdValid(const std::array<int, 3> &ijk) const
{
    for (int k = 0; k < getDimension(); ++k) {
        if (ijk[k] < 0 || ijk[k] >= m_nVertices1D[k]) {
            return false;
        }
    }
 
    return true;
}
 
void VolCartesian::_findCellFaceNeighs(long id, int face, const std::vector<long> *blackList, std::vector<long> *neighs) const
{
    long neighId = getCellFaceNeighsLinearId(id, face);
    if (neighId < 0) {
        return;
    } else if (blackList && utils::findInOrderedVector<long>(neighId, *blackList) != blackList->end()) {
        return;
    }
 
    neighs->push_back(neighId);
}
 
void VolCartesian::_findCellEdgeNeighs(long id, int edge, const std::vector<long> *blackList, std::vector<long> *neighs) const
{
    assert(isThreeDimensional());
    if (!isThreeDimensional()) {
        return;
    }
 
    // Diagonal neighbour
    std::array<int, 3> diagNeighIjk(getCellCartesianId(id) + m_edgeNeighDeltas[edge]);
    if (isCellCartesianIdValid(diagNeighIjk)) {
        long diagNeighId = getCellLinearId(diagNeighIjk);
        if (!blackList || utils::findInOrderedVector<long>(diagNeighId, *blackList) == blackList->end()) {
            utils::addToOrderedVector<long>(diagNeighId, *neighs);
        }
    }
 
    // Faces incident to the edge
    for (int face : m_edgeFaces[edge]) {
        _findCellFaceNeighs(id, face, blackList, neighs);
    }
}
 
void VolCartesian::_findCellVertexNeighs(long id, int vertex, const std::vector<long> *blackList, std::vector<long> *neighs) const
{
    std::array<int, 3> cellIjk   = getCellCartesianId(id);
    std::array<int, 3> vertexIjk = getVertexCartesianId(cellIjk, vertex);
 
    for (const auto &delta : m_vertexNeighDeltas) {
        // Get the Cartesian index
        std::array<int, 3> neighIjk(vertexIjk + delta);
        if (neighIjk == cellIjk) {
            continue;
        } else if (!isCellCartesianIdValid(neighIjk)) {
            continue;
        }
 
        // Get the linear neighbour index and, if it's not on the blacklist,
        // add it to the list of neighbours
        long neighId = getCellLinearId(neighIjk);
        if (!blackList || utils::findInOrderedVector<long>(neighId, *blackList) == blackList->end()) {
            utils::addToOrderedVector<long>(neighId, *neighs);
        }
    }
}
 
std::vector<long> VolCartesian::extractCellSubSet(std::array<int, 3> const &ijkMin, std::array<int, 3> const &ijkMax) const
{
    int nSubsetCells_x = ijkMax[0] - ijkMin[0] + 1;
    int nSubsetCells_y = ijkMax[1] - ijkMin[1] + 1;
    int nSubsetCells_z = ijkMax[2] - ijkMin[2] + 1;
 
    std::vector<long> ids;
    ids.resize(nSubsetCells_x * nSubsetCells_y * nSubsetCells_z);
 
    std::vector<long>::iterator it = ids.begin();
    for (int k = ijkMin[2]; k <= ijkMax[2]; ++k) {
        for (int j = ijkMin[1]; j <= ijkMax[1]; ++j) {
            for (int i = ijkMin[0]; i <= ijkMax[0]; ++i) {
                *it = getCellLinearId(i, j, k);
                ++it;
            }
        }
    }
 
    return ids;
}
 
std::vector<long> VolCartesian::extractCellSubSet(int idMin, int idMax) const
{
    return extractCellSubSet(getCellCartesianId(idMin), getCellCartesianId(idMax));
}
 
std::vector<long> VolCartesian::extractCellSubSet(std::array<double, 3> const &pointMin, std::array<double, 3> const &pointMax) const
{
    return extractCellSubSet(locatePointCartesian(pointMin), locatePointCartesian(pointMax));
}
 
std::vector<long> VolCartesian::extractVertexSubSet(std::array<int, 3> const &ijkMin, std::array<int, 3> const &ijkMax) const
{
    int nSubsetVertices_x = ijkMax[0] - ijkMin[0] + 1;
    int nSubsetVertices_y = ijkMax[1] - ijkMin[1] + 1;
    int nSubsetVertices_z = ijkMax[2] - ijkMin[2] + 1;
 
    std::vector<long> ids;
    ids.resize(nSubsetVertices_x * nSubsetVertices_y * nSubsetVertices_z);
 
    std::vector<long>::iterator it = ids.begin();
    for (int k = ijkMin[2]; k <= ijkMax[2]; ++k) {
        for (int j = ijkMin[1]; j <= ijkMax[1]; ++j) {
            for (int i = ijkMin[0]; i <= ijkMax[0]; ++i) {
                *it = getVertexLinearId(i, j, k);
                ++it;
            }
        }
    }
 
    return ids;
}
 
std::vector<long> VolCartesian::extractVertexSubSet(int idMin, int idMax) const
{
    return extractVertexSubSet(getVertexCartesianId(idMin), getVertexCartesianId(idMax));
}
 
std::vector<long> VolCartesian::extractVertexSubSet(std::array<double, 3> const &pointMin, std::array<double, 3> const &pointMax) const
{
    return extractVertexSubSet(locatePointCartesian(pointMin), locatePointCartesian(pointMax));
}
 
std::array<double, 3> VolCartesian::getOrigin() const
{
    return m_minCoords;
}
 
void VolCartesian::setOrigin(const std::array<double, 3> &origin)
{
    std::array<double, 3> translation = origin - getOrigin();
    translate(translation);
}
 
void VolCartesian::translate(const std::array<double, 3> &translation)
{
    for (int n = 0; n < 3; ++n) {
        m_minCoords[n] += translation[n];
        m_maxCoords[n] += translation[n];
 
        for (int i = 1; i < m_nVertices1D[n]; ++i) {
            m_vertexCoords[n][i] += translation[n];
        }
 
        for (int i = 1; i < m_nCells1D[n]; ++i) {
            m_cellCenters[n][i] += translation[n];
        }
    }
 
    setBoundingBox(m_minCoords, m_maxCoords);
 
    VolumeKernel::translate(translation);
}
 
void VolCartesian::rotate(const std::array<double, 3> &n0, const std::array<double, 3> &n1, double angle)
{
    BITPIT_UNUSED(n0);
    BITPIT_UNUSED(n1);
    BITPIT_UNUSED(angle);
 
    throw std::runtime_error("It is not possible to rotate a Cartesian patch.");
}
 
std::array<double, 3> VolCartesian::getLengths() const
{
    return (m_maxCoords - m_minCoords);
}
 
void VolCartesian::setLengths(const std::array<double, 3> &lengths)
{
    // Set the maximum coordinate
    m_maxCoords = m_minCoords + lengths;
 
    // Reset the bounding box
    setBoundingBox(m_minCoords, m_maxCoords);
 
    // If needed update the discretization
    if (m_nVertices > 0) {
        // Rebuild the discretization
        setDiscretization(m_nCells1D);
 
        // Rebuild vertex coordinates
        for (Vertex vertex : m_vertices) {
            long id = vertex.getId();
            vertex.setCoords(evalVertexCoords(id));
        }
    }
}
 
void VolCartesian::scale(const std::array<double, 3> &scaling, const std::array<double, 3> &center)
{
    for (int n = 0; n < 3; ++n) {
        m_minCoords[n] = center[n] + scaling[n] * (m_minCoords[n] - center[n]);
        m_maxCoords[n] = center[n] + scaling[n] * (m_maxCoords[n] - center[n]);
 
        for (int i = 1; i < m_nVertices1D[n]; ++i) {
            m_vertexCoords[n][i] = center[n] + scaling[n] * (m_vertexCoords[n][i] - center[n]);
        }
 
        for (int i = 1; i < m_nCells1D[n]; ++i) {
            m_cellCenters[n][i] = center[n] + scaling[n] * (m_cellCenters[n][i] - center[n]);
        }
 
        m_cellSpacings[n] *= scaling[n];
    }
 
    initializeCellVolume();
 
    initializeCellSize();
 
    initializeInterfaceArea();
 
    setBoundingBox(m_minCoords, m_maxCoords);
 
    VolumeKernel::scale(scaling, center);
}
 
std::vector<double> VolCartesian::convertToVertexData(const std::vector<double> &cellData) const
{
    int dimension = getDimension();
 
    int nContribs_x = 2;
    int nContribs_y = 2;
    int nContribs_z = dimension - 1;
 
    std::vector<int> nodeCounter(getVertexCount());
    std::vector<double> vertexData(getVertexCount());
    std::fill (vertexData.begin(), vertexData.end(), 0.);
 
    for (int k = 0; (isThreeDimensional()) ? (k < m_nCells1D[Vertex::COORD_Z]) : (k < 1); k++) {
        for (int j = 0; j < m_nCells1D[Vertex::COORD_Y]; ++j) {
            for (int i = 0; i < m_nCells1D[Vertex::COORD_X]; ++i) {
                // Get cell contribution
                long cellId = getCellLinearId(i, j, k);
                double cellContrib = cellData[cellId];
 
                // Eval vertex data
                for (int n = 0; n < nContribs_z; n++) {
                    int k_v = k + n;
                    for (int m = 0; m < nContribs_y; ++m) {
                        int j_v = j + m;
                        for (int l = 0; l < nContribs_x; ++l) {
                            int i_v = i + l;
 
                            // Get vertex index
                            long vertexId = getVertexLinearId(i_v, j_v, k_v);
 
                            // Sum cell contribution
                            vertexData[vertexId] += cellContrib;
                            nodeCounter[vertexId]++;
                        }
                    }
                }
            }
        }
    }
 
    // Average values
    for (int i = 0; i < getVertexCount(); ++i) {
        vertexData[i] /= nodeCounter[i];
    }
 
    return vertexData;
}
 
std::vector<double> VolCartesian::convertToCellData(const std::vector<double> &vertexData) const
{
    int dimension = getDimension();
 
    int nContribs_x = 2;
    int nContribs_y = 2;
    int nContribs_z = dimension - 1;
 
    std::vector<double> cellData(getCellCount());
    std::fill (cellData.begin(), cellData.end(), 0.);
 
    for (int k = 0; (isThreeDimensional()) ? (k < m_nCells1D[Vertex::COORD_Z]) : (k < 1); k++) {
        for (int j = 0; j < m_nCells1D[Vertex::COORD_Y]; ++j) {
            for (int i = 0; i < m_nCells1D[Vertex::COORD_X]; ++i) {
                // Eval cell data
                long cellId = getCellLinearId(i, j, k);
                for (int n = 0; n < nContribs_z; ++n) {
                    int k_v = k + n;
                    for (int m = 0; m < nContribs_y; ++m) {
                        int j_v = j + m;
                        for (int l = 0; l < nContribs_x; ++l) {
                            int i_v = i + l;
 
                            // Get vertex contribution
                            long vertexId = getVertexLinearId(i_v, j_v, k_v);
                            double vertexContrib = vertexData[vertexId];
 
                            // Sum vertex contribution
                            cellData[cellId] += vertexContrib;
                        }
                    }
                }
            }
        }
    }
 
    double weight = pow(0.5, dimension);
    for (int i = 0; i < getCellCount(); ++i) {
        cellData[i] *= weight;
    }
 
    return cellData;
}
 
 
int VolCartesian::linearCellInterpolation(const std::array<double,3> &point,
                                          std::vector<int> *stencil, std::vector<double> *weights) const
{
    std::array<int, 3> ijk_point = locatePointCartesian(point);
    if (ijk_point[0] < 0) {
        stencil->clear();
        weights->clear();
        return 0;
    }
 
    int dimension = getDimension();
 
    int nContribs_x = 1;
    int nContribs_y = 1;
    int nContribs_z = 1;
 
    std::array< std::array<int,2>, 3> cStencil;
    std::array< std::array<double,2>, 3> cWeights;
    for (int d = 0; d < dimension; ++d) {
        // Find cell index
        int index_point = ijk_point[d];
        if (point[d] < m_cellCenters[d][index_point]) {
            index_point = index_point - 1;
        }
 
        int index_next = index_point + 1;
 
        if (index_point < 0) {
            cStencil[d][0] = 0;
            cWeights[d][0] = 1.;
        } else if (index_next > m_nCells1D[d] - 1) {
            cStencil[d][0] = m_nCells1D[d] - 1;
            cWeights[d][0] = 1.;
        } else {
            if (d == 0) {
                nContribs_x = 2;
            } else if (d == 1) {
                nContribs_y = 2;
            } else {
                nContribs_z = 2;
            }
 
            cStencil[d][0] = index_point;
            cStencil[d][1] = index_next;
 
            cWeights[d][1] = (point[d] - m_cellCenters[d][index_point]) / m_cellSpacings[d]  ;
            cWeights[d][0] = 1.0 - cWeights[d][1];
        }
    }
 
    for (int d = dimension; d < 3; ++d) {
        cStencil[d][0] = 0;
        cWeights[d][0] = 1.;
    }
 
    int stencilSize = nContribs_x * nContribs_y * nContribs_z;
    stencil->resize(stencilSize);
    weights->resize(stencilSize);
 
    std::vector<int>::iterator itrStencil = stencil->begin();
    std::vector<double>::iterator itrWeights = weights->begin();
    for (int k = 0; k < nContribs_z; ++k) {
        for (int j = 0; j < nContribs_y; ++j) {
            for (int i = 0; i < nContribs_x; ++i) {
                int &is = cStencil[0][i];
                int &js = cStencil[1][j];
                int &ks = cStencil[2][k];
 
                double &iw = cWeights[0][i];
                double &jw = cWeights[1][j];
                double &kw = cWeights[2][k];
 
                *itrStencil = getCellLinearId(is, js, ks);
                *itrWeights = iw *jw * kw;
 
                ++itrStencil;
                ++itrWeights;
            }
        }
    }
 
    return stencilSize;
}
 
int VolCartesian::linearVertexInterpolation(const std::array<double,3> &point,
                                            std::vector<int> *stencil, std::vector<double> *weights) const
{
    std::array<int, 3> ijk_point = locatePointCartesian(point);
    if (ijk_point[0] < 0) {
        stencil->clear();
        weights->clear();
        return 0;
    }
 
    int dimension = getDimension();
 
    int nContribs_x = 2;
    int nContribs_y = 2;
    int nContribs_z = dimension - 1;
 
    std::array< std::array<int,2>, 3> cStencil;
    std::array< std::array<double,2>, 3> cWeights;
    for (int d = 0; d < dimension; ++d) {
        int index_point = ijk_point[d];
        int index_next  = index_point +1;
 
        cStencil[d][0] = index_point;
        cStencil[d][1] = index_next;
 
        cWeights[d][1] = (point[d] - m_vertexCoords[d][index_point]) / m_cellSpacings[d];
        cWeights[d][0] = 1.0 - cWeights[d][1];
    }
 
    for (int d = dimension; d < 3; ++d) {
        cStencil[d][0] = 0;
        cWeights[d][0] = 1.;
    }
 
    int stencilSize = uipow(2, dimension);
    stencil->resize(stencilSize);
    weights->resize(stencilSize);
 
    std::vector<int>::iterator itrStencil    = stencil->begin();
    std::vector<double>::iterator itrWeights = weights->begin();
    for (int k = 0; k < nContribs_z; ++k) {
        for (int j = 0; j < nContribs_y; ++j) {
            for (int i = 0; i < nContribs_x; ++i) {
                int &is = cStencil[0][i];
                int &js = cStencil[1][j];
                int &ks = cStencil[2][k];
 
                double &iw = cWeights[0][i];
                double &jw = cWeights[1][j];
                double &kw = cWeights[2][k];
 
                *itrStencil = getVertexLinearId(is, js, ks);
                *itrWeights = iw * jw * kw;
 
                ++itrStencil;
                ++itrWeights;
            }
        }
    }
 
    return stencilSize;
}
 
std::array<double, 3> VolCartesian::evalCellCentroid(long id) const
{
    std::array<int, 3> ijk = getCellCartesianId(id);
 
    return evalCellCentroid(ijk);
}
 
std::array<double, 3> VolCartesian::evalCellCentroid(const std::array<int, 3> &ijk) const
{
    std::array<double, 3> centroid;
    for (int d = 0; d < 3; ++d) {
        if (ijk[d] >= 0) {
            centroid[d] = m_cellCenters[d][ijk[d]];
        } else {
            centroid[d] = m_minCoords[d];
        }
    }
 
    return centroid;
}
 
const std::vector<double> & VolCartesian::getCellCentroids(int direction) const
{
    return m_cellCenters[direction];
}
 
std::array<int, 3> VolCartesian::getCellFaceNeighsCartesianId(long id, int face) const
{
    int neighSide      = face % 2;
    int neighDirection = static_cast<int>(std::floor(face / 2));
 
    std::array<int, 3> neighIjk(getCellCartesianId(id));
    if (neighSide == 0) {
        neighIjk[neighDirection]--;
    } else {
        neighIjk[neighDirection]++;
    }
 
    return neighIjk;
}
 
long VolCartesian::getCellFaceNeighsLinearId(long id, int face) const
{
    std::array<int, 3> neighIjk = getCellFaceNeighsCartesianId(id, face);
 
    long neighId;
    if (isCellCartesianIdValid(neighIjk)) {
        neighId = getCellLinearId(neighIjk);
    } else {
        neighId = Cell::NULL_ID;
    }
 
    return neighId;
}
 
}