Octant.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/>.
 *
\*---------------------------------------------------------------------------*/
 
// =================================================================================== //
// INCLUDES                                                                            //
// =================================================================================== //
#include "Octant.hpp"
#include "morton.hpp"
 
#include <algorithm>
#include <cassert>
#include <cmath>
#include <iostream>
 
namespace bitpit {
 
IBinaryStream& operator>>(IBinaryStream &buffer, Octant &octant)
{
    uint8_t dimensions;
    buffer >> dimensions;
 
    uint8_t level;
    buffer >> level;
 
    octant.initialize(dimensions, level, true);
 
    buffer >> octant.m_morton;
 
    buffer >> octant.m_marker;
 
    buffer >> octant.m_ghost;
 
    for(int i = 0; i < Octant::INFO_ITEM_COUNT; ++i){
        bool value;
        buffer >> value;
        octant.m_info[i] = value;
    }
 
    return buffer;
}
 
OBinaryStream& operator<<(OBinaryStream  &buffer, const Octant &octant)
{
    buffer << octant.m_dim;
    buffer << octant.m_level;
 
    buffer << octant.m_morton;
 
    buffer << octant.m_marker;
 
    buffer << octant.m_ghost;
 
    for(int i = 0; i < Octant::INFO_ITEM_COUNT; ++i){
        buffer << (bool) octant.m_info[i];
    }
 
    return buffer;
}
 
// =================================================================================== //
// NAME SPACES                                                                         //
// =================================================================================== //
using namespace std;
 
// =================================================================================== //
// CLASS IMPLEMENTATION                                                                    //
// =================================================================================== //
 
// =================================================================================== //
// STATIC AND CONSTANT
// =================================================================================== //
 
const TreeConstants::Instances Octant::sm_treeConstants = TreeConstants::instances();
 
// =================================================================================== //
// CONSTRUCTORS AND OPERATORS
// =================================================================================== //
 
Octant::Octant(){
    initialize();
};
 
Octant::Octant(uint8_t dim){
    initialize(dim, 0, true);
};
 
Octant::Octant(uint8_t dim, uint8_t level, uint64_t morton){
    initialize(dim, level, true);
 
    // Set the morton
    m_morton = morton;
 
};
 
Octant::Octant(uint8_t dim, uint8_t level, int32_t x, int32_t y, int32_t z){
    initialize(dim, level, true);
 
    // Set the morton
    m_morton = PABLO::computeMorton(m_dim, x, y, z);
 
};
 
Octant::Octant(bool bound, uint8_t dim, uint8_t level, int32_t x, int32_t y, int32_t z){
    initialize(dim, level, bound);
 
    // Set the morton
    m_morton = PABLO::computeMorton(m_dim, x, y, z);
 
};
 
bool Octant::operator ==(const Octant & oct2){
    bool check = true;
    check = check && (m_dim == oct2.m_dim);
    check = check && (m_morton == oct2.m_morton);
    check = check && (m_level == oct2.m_level);
    return check;
}
 
// =================================================================================== //
// METHODS
// =================================================================================== //
 
void
Octant::initialize() {
    initialize(0, 0, false);
}
 
void
Octant::initialize(uint8_t dim, uint8_t level, bool bound) {
    m_dim   = dim;
    m_level = level;
 
    // Reset the marker
    m_marker = 0;
 
    // Reset the morton
    m_morton = 0;
 
    // Initialize octant info
    m_info.reset();
    m_info[OctantInfo::INFO_BALANCED] = true;
    m_ghost = -1;
 
    // If this is the root octant we need to set the boundary condition bound
    // for faces
    if (m_dim >= 2 && m_level == 0) {
        uint8_t nf = m_dim*2;
        for (uint8_t i=0; i<nf; i++){
            m_info[i] = bound;
        }
    }
};
 
// =================================================================================== //
// BASIC GET/SET METHODS
// =================================================================================== //
 
uint32_t
Octant::getDim() const{return m_dim;};
 
u32array3
Octant::getLogicalCoordinates() const{
    return {{getLogicalCoordinates(0), getLogicalCoordinates(1), getLogicalCoordinates(2)}};
};
 
uint32_t
Octant::getLogicalCoordinates(int coord) const{
    return PABLO::computeCoordinate(m_dim, m_morton, coord);
};
 
uint8_t
Octant::getLevel() const{return m_level;};
 
int8_t
Octant::getMarker() const{return m_marker;};
 
bool
Octant::getBound(uint8_t face) const{
    return (m_info[OctantInfo::INFO_BOUNDFACE0 + face]);
};
 
bool
Octant::getEdgeBound(uint8_t edge) const{
    for (int face : sm_treeConstants[m_dim].edgeFace[edge]) {
        if (getBound(face)) {
            return true;
        }
    }
 
    return false;
};
 
bool
Octant::getNodeBound(uint8_t node) const{
    for (int face : sm_treeConstants[m_dim].nodeFace[node]) {
        if (getBound(face)) {
            return true;
        }
    }
 
    return false;
};
 
bool
Octant::getBound() const{
    for (int i = 0; i < sm_treeConstants[m_dim].nFaces; ++i) {
        if (getBound(i)) {
            return true;
        }
    }
 
    return false;
};
 
void
Octant::setBound(uint8_t face) {
    m_info[INFO_BOUNDFACE0 + face] = true;
};
 
bool
Octant::getPbound(uint8_t face) const{
    return m_info[INFO_PBOUNDFACE0 + face];
};
 
bool
Octant::getPbound() const{
    for (int i = 0; i < sm_treeConstants[m_dim].nFaces; ++i) {
        if (getPbound(i)) {
            return true;
        }
    }
 
    return false;
};
 
bool
Octant::getIsNewR() const{return m_info[OctantInfo::INFO_NEW4REFINEMENT];};
 
bool
Octant::getIsNewC() const{return m_info[OctantInfo::INFO_NEW4COARSENING];};
 
bool
Octant::getIsGhost() const{return (m_ghost >= 0);};
 
int
Octant::getGhostLayer() const{return m_ghost;};
 
bool
Octant::getBalance() const{return (m_info[OctantInfo::INFO_BALANCED]);};
 
void
Octant::setMarker(int8_t marker){
    this->m_marker = marker;
};
 
void
Octant::setBalance(bool balance){
    m_info[OctantInfo::INFO_BALANCED] = balance;
};
 
void
Octant::setLevel(uint8_t level){
    this->m_level = level;
};
 
void
Octant::setPbound(uint8_t face, bool flag){
    m_info[INFO_PBOUNDFACE0 + face] = flag;
};
 
void
Octant::setGhostLayer(int ghostLayer){
    m_ghost = ghostLayer;
};
 
 
// =================================================================================== //
// OTHER GET/SET METHODS
// =================================================================================== //
 
uint32_t
Octant::getLogicalSize() const{
    return sm_treeConstants[m_dim].lengths[m_level];
};
 
uint64_t
Octant::getLogicalArea() const{
    return sm_treeConstants[m_dim].areas[m_level];
};
 
uint64_t
Octant::getLogicalVolume() const{
    return sm_treeConstants[m_dim].volumes[m_level];
};
 
// =================================================================================== //
 
darray3
Octant::getLogicalCenter() const{
    double dh = double(getLogicalSize())*0.5;
    darray3 center = {{0., 0., 0.}};
    center[0] = getLogicalCoordinates(0) + dh;
    center[1] = getLogicalCoordinates(1) + dh;
    center[2] = getLogicalCoordinates(2) + double(m_dim-2)*dh;
    return center;
};
 
// =================================================================================== //
 
darray3
Octant::getLogicalFaceCenter(uint8_t iface) const{
    assert(iface < m_dim*2);
 
    double dh_2 = double(getLogicalSize())*0.5;
    darray3 center = {{0., 0., 0.}};
    center[0] = getLogicalCoordinates(0) + (double)sm_treeConstants[m_dim].faceDisplacements[iface][0] * dh_2;
    center[1] = getLogicalCoordinates(1) + (double)sm_treeConstants[m_dim].faceDisplacements[iface][1] * dh_2;
    center[2] = getLogicalCoordinates(2) + double(m_dim-2) * (double)sm_treeConstants[m_dim].faceDisplacements[iface][2] * dh_2;
 
    return center;
};
 
// =================================================================================== //
 
darray3
Octant::getLogicalEdgeCenter(uint8_t iedge) const{
    double  dh_2;
    darray3 center = {{0., 0., 0.}};
 
    dh_2 = double(getLogicalSize())*0.5;
    center[0] = getLogicalCoordinates(0) + (double)sm_treeConstants[m_dim].edgeDisplacements[iedge][0] * dh_2;
    center[1] = getLogicalCoordinates(1) + (double)sm_treeConstants[m_dim].edgeDisplacements[iedge][1] * dh_2;
    center[2] = getLogicalCoordinates(2) + double(m_dim-2) * (double)sm_treeConstants[m_dim].edgeDisplacements[iedge][2] * dh_2;
    return center;
};
 
// =================================================================================== //
 
void
Octant::getLogicalNodes(u32arr3vector & nodes) const{
    uint8_t     i;
    uint32_t    dh;
    uint8_t nn = uint8_t(1)<<m_dim;
 
    dh = getLogicalSize();
    nodes.resize(nn);
 
    u32array3 coords = getLogicalCoordinates();
 
    for (i = 0; i < nn; i++){
        nodes[i][0] = coords[0] + uint32_t(sm_treeConstants[m_dim].nodeCoordinates[i][0])*dh;
        nodes[i][1] = coords[1] + uint32_t(sm_treeConstants[m_dim].nodeCoordinates[i][1])*dh;
        nodes[i][2] = coords[2] + uint32_t(sm_treeConstants[m_dim].nodeCoordinates[i][2])*dh;
    }
};
 
u32arr3vector
Octant::getLogicalNodes() const{
    u32arr3vector nodes;
    getLogicalNodes(nodes);
    return nodes;
};
 
void        Octant::getLogicalNode(u32array3 & node, uint8_t inode) const{
    uint32_t    dh;
 
    dh = getLogicalSize();
 
    node = getLogicalCoordinates();
 
    node[0] += sm_treeConstants[m_dim].nodeCoordinates[inode][0]*dh;
    node[1] += sm_treeConstants[m_dim].nodeCoordinates[inode][1]*dh;
    node[2] += sm_treeConstants[m_dim].nodeCoordinates[inode][2]*dh;
 
};
 
u32array3       Octant::getLogicalNode(uint8_t inode) const{
    u32array3 node;
    getLogicalNode(node, inode);
    return node;
};
 
void        Octant::getNormal(uint8_t iface, i8array3 & normal, const int8_t (&normals)[6][3]) const{
    uint8_t     i;
    for (i = 0; i < 3; i++){
        normal[i] = normals[iface][i];
    }
};
 
 
uint64_t    Octant::getMorton() const{
    return m_morton;
};
 
uint64_t    Octant::computeLastDescMorton() const {
    u32array3 lastDescCoordinates = computeLastDescCoordinates();
    return PABLO::computeMorton(m_dim, lastDescCoordinates[0], lastDescCoordinates[1], lastDescCoordinates[2]);
};
 
u32array3   Octant::computeLastDescCoordinates() const {
    u32array3 lastDescCoords = getLogicalCoordinates();
    for (int i=0; i<m_dim; i++){
        lastDescCoords[i] += (uint32_t(1) << (sm_treeConstants[m_dim].maxLevel - m_level)) - 1;
    }
 
    return lastDescCoords;
};
 
uint64_t    Octant::computeFatherMorton() const {
    u32array3 fatherCoordinates = computeFatherCoordinates();
    return PABLO::computeMorton(m_dim, fatherCoordinates[0], fatherCoordinates[1], fatherCoordinates[2]);
};
 
u32array3   Octant::computeFatherCoordinates() const {
    u32array3 fatherCoordinates = getLogicalCoordinates();
    for (int i=0; i<m_dim; i++){
        fatherCoordinates[i] -= fatherCoordinates[i]%(uint32_t(1) << (sm_treeConstants[m_dim].maxLevel - max(0,(m_level-1))));
    }
    return fatherCoordinates;
};
 
uint64_t    Octant::computeNodePersistentKey(uint8_t inode) const{
 
    u32array3 node = getLogicalNode(inode);
 
    return computeNodePersistentKey(node);
};
 
uint64_t    Octant::computeNodePersistentKey(const u32array3 &node) const{
 
    return PABLO::computeXYZKey(m_dim, node[0], node[1], node[2]);
};
 
unsigned int Octant::getBinarySize()
{
    unsigned int binarySize = 0;
    binarySize += sizeof(uint8_t); // dimensions
    binarySize += sizeof(uint8_t); // level
    binarySize += sizeof(uint64_t); // morton
    binarySize += sizeof(int8_t); // marker
    binarySize += sizeof(int); // ghost layer
    binarySize += INFO_ITEM_COUNT * sizeof(bool); // info
 
    return binarySize;
}
 
// =================================================================================== //
// OTHER METHODS
// =================================================================================== //
 
Octant  Octant::buildLastDesc() const {
    uint64_t lastDescMorton = computeLastDescMorton();
    Octant lastDesc(m_dim, sm_treeConstants[m_dim].maxLevel, lastDescMorton);
    return lastDesc;
};
 
// =================================================================================== //
 
Octant  Octant::buildFather() const {
    uint64_t fatherMorton = computeFatherMorton();
    Octant father(m_dim, max(0,m_level-1), fatherMorton);
    return father;
};
 
// =================================================================================== //
 
uint8_t Octant::countChildren() const {
    if (this->m_level < sm_treeConstants[m_dim].maxLevel){
        return sm_treeConstants[m_dim].nChildren;
    } else {
        return 0;
    }
}
 
// =================================================================================== //
 
vector< Octant >    Octant::buildChildren() const {
 
    uint8_t nchildren = countChildren();
    std::vector< Octant > children(nchildren);
    buildChildren(children.data());
 
    return children;
}
 
// =================================================================================== //
 
void    Octant::buildChildren(Octant *children) const {
 
    u32array3 coords = getLogicalCoordinates();
 
    int nChildren = countChildren();
    for (int i=0; i<nChildren; ++i){
        // Octant information
        uint8_t xf;
        uint8_t yf;
        uint8_t zf;
        uint8_t dx;
        uint8_t dy;
        uint8_t dz;
        switch (i) {
 
        case 0 :
            dx = 0;
            dy = 0;
            dz = 0;
 
            xf = 1;
            yf = 3;
            zf = 5;
 
            break;
 
        case 1 :
            dx = 1;
            dy = 0;
            dz = 0;
 
            xf = 0;
            yf = 3;
            zf = 5;
 
            break;
 
        case 2 :
            dx = 0;
            dy = 1;
            dz = 0;
 
            xf = 1;
            yf = 2;
            zf = 5;
 
            break;
 
        case 3 :
            dx = 1;
            dy = 1;
            dz = 0;
 
            xf = 0;
            yf = 2;
            zf = 5;
 
            break;
 
        case 4 :
            dx = 0;
            dy = 0;
            dz = 1;
 
            xf = 1;
            yf = 3;
            zf = 4;
 
            break;
 
        case 5 :
            dx = 1;
            dy = 0;
            dz = 1;
 
            xf = 0;
            yf = 3;
            zf = 4;
 
            break;
 
        case 6 :
            dx = 0;
            dy = 1;
            dz = 1;
 
            xf = 1;
            yf = 2;
            zf = 4;
 
            break;
 
        case 7 :
            dx = 1;
            dy = 1;
            dz = 1;
 
            xf = 0;
            yf = 2;
            zf = 4;
 
            break;
 
        default:
            BITPIT_UNREACHABLE("The maximum number of children is 8.");
 
        }
 
        // Create octant
        children[i] = Octant(*this);
        Octant &oct = children[i];
 
        oct.setMarker(std::max(0, oct.m_marker - 1));
        oct.setLevel(oct.m_level + 1);
 
        uint32_t dh = oct.getLogicalSize();
        oct.m_morton = PABLO::computeMorton(m_dim, coords[0] + dh * dx, coords[1] + dh * dy, coords[2] + dh * dz);
 
        oct.m_info[OctantInfo::INFO_NEW4REFINEMENT] = true;
 
        oct.m_info[INFO_BOUNDFACE0 + xf] = false;
        oct.m_info[INFO_BOUNDFACE0 + yf] = false;
        oct.m_info[INFO_BOUNDFACE0 + zf] = false;
 
        oct.m_info[INFO_PBOUNDFACE0 + xf] = false;
        oct.m_info[INFO_PBOUNDFACE0 + yf] = false;
        oct.m_info[INFO_PBOUNDFACE0 + zf] = false;
    }
};
 
uint8_t Octant::getFamilySplittingNode() const {
 
    u32array3 xx = getLogicalCoordinates();
 
    //Delta to build father (use the boolean negation to identify the splitting node)
    bool delta[3];
    delta[2] = 0;
    for (int i=0; i<m_dim; i++){
        delta[i] = ((xx[i]%(uint32_t(1) << (sm_treeConstants[m_dim].maxLevel - max(0,(m_level-1))))) == 0);
    }
    return sm_treeConstants[m_dim].nodeFromCoordinates[delta[0]][delta[1]][delta[2]];
};
 
 
 
}