patch_skd_tree.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 <queue>
 
#include "bitpit_CG.hpp"
 
#include "patch_skd_tree.hpp"
 
namespace bitpit {
 
SkdPatchInfo::SkdPatchInfo(const PatchKernel *patch, const std::vector<std::size_t> *cellRawIds)
    : m_patch(patch), m_cellRawIds(cellRawIds)
{
    if (!m_patch) {
        throw std::runtime_error("Unable to initialize the patch info. Provided patch is not valid.");
    }
}
 
void SkdPatchInfo::buildCache()
{
    PatchKernel::CellConstRange cellRange(m_patch->internalCellConstBegin(), m_patch->internalCellConstEnd());
 
    buildCache(cellRange);
}
 
void SkdPatchInfo::buildCache(const PatchKernel::CellConstRange &cellRange)
{
    m_cellBoxes = std::unique_ptr<BoxCache>(new BoxCache(1, &(m_patch->getCells())));
    for (auto itr = cellRange.cbegin(); itr != cellRange.cend(); ++itr) {
        std::size_t rawCellId = itr.getRawIndex();
 
        // Bounding box
        std::array<std::array<double, 3>, 2> *cellBox = m_cellBoxes->rawData(rawCellId);
        m_patch->evalElementBoundingBox(*itr, cellBox->data(), cellBox->data() + 1);
    }
}
 
void SkdPatchInfo::destroyCache()
{
    m_cellBoxes.reset();
}
 
const PatchKernel & SkdPatchInfo::getPatch() const
{
    return *m_patch;
}
 
const std::vector<std::size_t> & SkdPatchInfo::getCellRawIds() const
{
    return *m_cellRawIds;
}
 
std::size_t SkdPatchInfo::getCellRawId(std::size_t n) const
{
    return (*m_cellRawIds)[n];
}
 
const std::array<std::array<double, 3>, 2> & SkdPatchInfo::getCachedBox(std::size_t rawId) const
{
    return m_cellBoxes->rawAt(rawId);
}
 
std::array<double, 3> SkdPatchInfo::evalCachedBoxMean(std::size_t rawId) const
{
    const std::array<std::array<double, 3>, 2> &cellBox = m_cellBoxes->rawAt(rawId);
 
    return 0.5 * (cellBox[0] + cellBox[1]);
}
 
double SkdPatchInfo::evalCachedBoxMean(std::size_t rawId, int direction) const
{
    const std::array<std::array<double, 3>, 2> &cellBox = m_cellBoxes->rawAt(rawId);
 
    return 0.5 * (cellBox[0][direction] + cellBox[1][direction]);
}
 
SkdBox::SkdBox()
    : m_boxMin{  std::numeric_limits<double>::max(),   std::numeric_limits<double>::max(),   std::numeric_limits<double>::max()},
      m_boxMax{- std::numeric_limits<double>::max(), - std::numeric_limits<double>::max(), - std::numeric_limits<double>::max()}
{
}
 
SkdBox::SkdBox(const std::array<double,3> &boxMin, const std::array<double,3> &boxMax)
    : m_boxMin(boxMin),
      m_boxMax(boxMax)
{
}
 
bool SkdBox::isEmpty() const
{
    for (int d = 0; d < 3; ++d) {
        if (m_boxMin[d] > m_boxMax[d]) {
            return true;
        }
    }
 
    return false;
}
 
const std::array<double, 3> & SkdBox::getBoxMin() const
{
    return m_boxMin;
}
 
const std::array<double, 3> & SkdBox::getBoxMax() const
{
    return m_boxMax;
}
 
double SkdBox::evalPointMinDistance(const std::array<double, 3> &point, double emptyDistance) const
{
    if (isEmpty()) {
        return emptyDistance;
    }
 
    return evalPointMinDistance(point);
}
 
double SkdBox::evalPointMinDistance(const std::array<double, 3> &point) const
{
    return std::sqrt(evalPointMinSquareDistance(point));
}
 
double SkdBox::evalPointMinSquareDistance(const std::array<double, 3> &point, double emptySquareDistance) const
{
    if (isEmpty()) {
        return emptySquareDistance;
    }
 
    return evalPointMinSquareDistance(point);
}
 
double SkdBox::evalPointMinSquareDistance(const std::array<double, 3> &point) const
{
    double squareDistance = 0.;
    for (int d = 0; d < 3; ++d) {
        squareDistance += std::pow(std::max({0., m_boxMin[d] - point[d], point[d] - m_boxMax[d]}), 2);
    }
 
    return squareDistance;
}
 
double SkdBox::evalPointMaxDistance(const std::array<double, 3> &point, double emptyDistance) const
{
    if (isEmpty()) {
        return emptyDistance;
    }
 
    return evalPointMaxDistance(point);
}
 
double SkdBox::evalPointMaxDistance(const std::array<double, 3> &point) const
{
    return std::sqrt(evalPointMaxSquareDistance(point));
}
 
double SkdBox::evalPointMaxSquareDistance(const std::array<double, 3> &point, double emptySquareDistance) const
{
    if (isEmpty()) {
        return emptySquareDistance;
    }
 
    return evalPointMaxSquareDistance(point);
}
 
double SkdBox::evalPointMaxSquareDistance(const std::array<double, 3> &point) const
{
    double squareDistance = 0.;
    for (int d = 0; d < 3; ++d) {
        squareDistance += std::pow(std::max(point[d] - m_boxMin[d], m_boxMax[d] - point[d]), 2);
    }
 
    return squareDistance;
}
 
bool SkdBox::boxContainsPoint(const std::array<double, 3> &point, double offset) const
{
    if (isEmpty()) {
        return false;
    }
 
    for (int d = 0; d < 3; d++) {
        if (point[d] < (m_boxMin[d] - offset)) {
            return false;
        }
 
        if (point[d] > (m_boxMax[d] + offset)) {
            return false;
        }
    }
 
    return true;
}
 
bool SkdBox::boxIntersectsSphere(const std::array<double, 3> &center, double radius) const
{
    // If the box contains the center it will also intersect the sphere
    if (boxContainsPoint(center, 0.)) {
        return true;
    }
 
    // Check if the distance between the center and its projection on
    // the bounding box is smaller than the radius of the sphere.
    std::array<double, 3> delta;
    for (int d = 0; d < 3; d++) {
        delta[d] = std::min(std::max(center[d], m_boxMin[d]), m_boxMax[d]) - center[d];
    }
    double distance = dotProduct(delta, delta);
 
    return (distance < radius * radius);
}
 
const std::size_t SkdNode::NULL_ID = std::numeric_limits<std::size_t>::max();
 
SkdNode::SkdNode()
    : m_patchInfo(nullptr),
      m_cellRangeBegin(0), m_cellRangeEnd(0),
      m_children({{NULL_ID, NULL_ID}})
{
}
 
SkdNode::SkdNode(const SkdPatchInfo *patchInfo, std::size_t cellRangeBegin, std::size_t cellRangeEnd)
    : m_patchInfo(patchInfo),
      m_cellRangeBegin(cellRangeBegin), m_cellRangeEnd(cellRangeEnd),
      m_children({{NULL_ID, NULL_ID}})
{
    initializeBoundingBox();
}
 
void SkdNode::initializeBoundingBox()
{
    // Early return if there are no cells
    if (m_cellRangeBegin == m_cellRangeEnd){
        m_boxMin.fill(  std::numeric_limits<double>::max());
        m_boxMax.fill(- std::numeric_limits<double>::max());
 
        return;
    }
 
    // Evaluate the bounding box
    const std::vector<std::size_t> &cellRawIds = m_patchInfo->getCellRawIds();
 
    const std::array<std::array<double, 3>, 2> &firstCellBox = m_patchInfo->getCachedBox(cellRawIds[m_cellRangeBegin]);
    m_boxMin = firstCellBox[0];
    m_boxMax = firstCellBox[1];
    for (std::size_t n = m_cellRangeBegin + 1; n < m_cellRangeEnd; n++) {
        const std::array<std::array<double, 3>, 2> &cellBox = m_patchInfo->getCachedBox(cellRawIds[n]);
        for (int d = 0; d < 3; ++d) {
            m_boxMin[d] = std::min(cellBox[0][d], m_boxMin[d]);
            m_boxMax[d] = std::max(cellBox[1][d], m_boxMax[d]);
        }
    }
 
    // Inlfate the bounding box by a small epsilon
    //
    // The small inflation allows us to test intersection of a query
    // against a safety box and be assured that non-intersection of
    // with the box implies non-intersection with all the contents
    // in the box, even when computations are susceptible to round-off
    // error.
    double tolerance = m_patchInfo->getPatch().getTol();
 
    m_boxMin -= tolerance;
    m_boxMax += tolerance;
}
 
std::size_t SkdNode::getCellCount() const
{
    return (m_cellRangeEnd - m_cellRangeBegin);
}
 
long SkdNode::getCell(std::size_t n) const
{
    const PatchKernel &patch = m_patchInfo->getPatch();
    const PiercedKernel<long> &cellKernel = patch.getCells().getKernel();
 
    const std::size_t rawCellId = m_patchInfo->getCellRawId(m_cellRangeBegin + n);
    long cellId = cellKernel.rawFind(rawCellId).getId();
 
    return cellId;
}
 
std::vector<long> SkdNode::getCells() const
{
    const PatchKernel &patch = m_patchInfo->getPatch();
    const std::vector<std::size_t> &cellRawIds = m_patchInfo->getCellRawIds();
    const PiercedKernel<long> &cellKernel = patch.getCells().getKernel();
 
    std::size_t nCells = getCellCount();
    std::vector<long> cellList(nCells);
    for (std::size_t n = m_cellRangeBegin; n < m_cellRangeEnd; ++n) {
        const std::size_t rawCellId = cellRawIds[n];
        long cellId = cellKernel.rawFind(rawCellId).getId();
 
        cellList[n - m_cellRangeBegin] = cellId;
    }
 
    return cellList;
}
 
const SkdBox & SkdNode::getBoundingBox() const
{
    return *this;
}
 
std::array<double, 3> SkdNode::evalBoxWeightedMean() const
{
    const std::vector<std::size_t> &cellRawIds = m_patchInfo->getCellRawIds();
 
    std::array<double, 3> boxWeightedMean = m_patchInfo->evalCachedBoxMean(cellRawIds[m_cellRangeBegin]);
    for (std::size_t n = m_cellRangeBegin + 1; n < m_cellRangeEnd; ++n) {
        const std::size_t rawCellId = cellRawIds[n];
        boxWeightedMean += m_patchInfo->evalCachedBoxMean(rawCellId);
    }
    boxWeightedMean /= (double) getCellCount();
 
    return boxWeightedMean;
}
 
double SkdNode::evalBoxWeightedMean(int direction) const
{
    const std::vector<std::size_t> &cellRawIds = m_patchInfo->getCellRawIds();
 
    double boxWeightedMean = m_patchInfo->evalCachedBoxMean(cellRawIds[m_cellRangeBegin], direction);
    for (std::size_t n = m_cellRangeBegin + 1; n < m_cellRangeEnd; ++n) {
        const std::size_t rawCellId = cellRawIds[n];
        boxWeightedMean += m_patchInfo->evalCachedBoxMean(rawCellId, direction);
    }
    boxWeightedMean /= getCellCount();
 
    return boxWeightedMean;
}
 
bool SkdNode::isLeaf() const
{
    for (int i = CHILD_BEGIN; i != CHILD_END; ++i) {
        SkdNode::ChildLocation childLocation = static_cast<ChildLocation>(i);
        if (hasChild(childLocation)) {
            return false;
        }
    }
 
    return true;
}
 
bool SkdNode::hasChild(ChildLocation child) const
{
    return (m_children[static_cast<std::size_t>(child)] != NULL_ID);
}
 
std::size_t SkdNode::getChildId(ChildLocation child) const
{
    return m_children[static_cast<std::size_t>(child)];
}
 
double SkdNode::evalPointDistance(const std::array<double, 3> &point, bool interiorCellsOnly) const
{
    long id;
    double distance;
 
    findPointClosestCell(point, interiorCellsOnly, &id, &distance);
 
    return distance;
}
 
void SkdNode::findPointClosestCell(const std::array<double, 3> &point, bool interiorCellsOnly,
                                   long *id, double *closestDistance) const
{
    *id       = Cell::NULL_ID;
    *closestDistance = std::numeric_limits<double>::max();
 
    updatePointClosestCell(point, interiorCellsOnly, id, closestDistance);
}
 
void SkdNode::updatePointClosestCell(const std::array<double, 3> &point, bool interiorCellsOnly,
                                     long *closestId, double *closestDistance) const
{
    if (getCellCount() == 0) {
        return;
    }
 
    const PatchKernel &patch = m_patchInfo->getPatch();
    const PiercedVector<Cell> &cells = patch.getCells();
    const std::vector<std::size_t> &cellRawIds = m_patchInfo->getCellRawIds();
 
    for (std::size_t n = m_cellRangeBegin; n < m_cellRangeEnd; n++) {
        std::size_t cellRawId = cellRawIds[n];
        const Cell &cell = cells.rawAt(cellRawId);
        if (interiorCellsOnly && !cell.isInterior()) {
            continue;
        }
 
        // Evaluate point projection
        int nCellVertices = cell.getVertexCount();
        BITPIT_CREATE_WORKSPACE(cellVertexCoordinates, std::array<double BITPIT_COMMA 3>, nCellVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
        patch.getElementVertexCoordinates(cell, cellVertexCoordinates);
 
        double cellDistance;
        std::array<double, 3> cellProjection;
        cell.evalPointProjection(point, cellVertexCoordinates, &cellProjection, &cellDistance);
 
        // Check if the cell is closest that the current closest cell
        bool isCellClosest = false;
        if (utils::DoubleFloatingEqual()(cellDistance, *closestDistance, patch.getTol(), patch.getTol())) {
            if (*closestId == Cell::NULL_ID) {
                isCellClosest = true;
            } else {
                // Project point ont the closest cell
                const Cell &closest = patch.getCell(*closestId);
 
                int nClosestCellVertices = closest.getVertexCount();
                BITPIT_CREATE_WORKSPACE(closestCellVertexCoordinates, std::array<double BITPIT_COMMA 3>, nClosestCellVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
                patch.getElementVertexCoordinates(closest, closestCellVertexCoordinates);
 
                double closestCellDistance;
                std::array<double, 3> closestCellProjection;
                closest.evalPointProjection(point, closestCellVertexCoordinates, &closestCellProjection, &closestCellDistance);
 
                // Find normal of the cells at the centroids
                //
                // To be more precise, the normals should beevaluated on projection
                // points.
                std::array<double, 3> cellNormal    = cell.evalNormal(cellVertexCoordinates);
                std::array<double, 3> closestNormal = closest.evalNormal(closestCellVertexCoordinates);
 
                // Find cell more aligned with respect to the point
                //
                // We want to chhose the cell whose normal is most aligned towards
                // the point.
                double cellAligement    = dotProduct(cellNormal, point - cellProjection);
                double closestAligement = dotProduct(closestNormal, point - closestCellProjection);
 
                // The cell will be closest if it is more aligned than the current closest cell
                if (cellAligement > closestAligement) {
                    isCellClosest = true;
                }
            }
        } else if (cellDistance < *closestDistance) {
            isCellClosest = true;
        }
 
        // Update the closest cell
        if (isCellClosest) {
            *closestId       = cell.getId();
            *closestDistance = cellDistance;
        }
    }
}
 
void SkdNode::updatePointClosestCells(const std::array<double, 3> &point, bool interiorCellsOnly,
                                      std::vector<long> *closestIds, double *closestDistance) const
{
    if (getCellCount() == 0) {
        return;
    }
 
    const PatchKernel &patch = m_patchInfo->getPatch();
    const PiercedVector<Cell> &cells = patch.getCells();
    const std::vector<std::size_t> &cellRawIds = m_patchInfo->getCellRawIds();
 
    for (std::size_t n = m_cellRangeBegin; n < m_cellRangeEnd; n++) {
        std::size_t cellRawId = cellRawIds[n];
        const Cell &cell = cells.rawAt(cellRawId);
        if (interiorCellsOnly && !cell.isInterior()) {
            continue;
        }
 
        // Evaluate point projection
        int nCellVertices = cell.getVertexCount();
        BITPIT_CREATE_WORKSPACE(cellVertexCoordinates, std::array<double BITPIT_COMMA 3>, nCellVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
        patch.getElementVertexCoordinates(cell, cellVertexCoordinates);
 
        double cellDistance;
        std::array<double, 3> cellProjection;
        cell.evalPointProjection(point, cellVertexCoordinates, &cellProjection, &cellDistance);
 
        // Update the closest cell
        //
        // If the cell is closest that the current closest cells, the list should be cleard and
        // the current cell should be added to it. If the cell has the same distance of the
        // closest cells, it should be simply added to the list. If the cell i farther than the
        // closest one, there is nothing to do.
        if (utils::DoubleFloatingEqual()(cellDistance, *closestDistance, patch.getTol(), patch.getTol())) {
            closestIds->push_back(cell.getId());
        } else if (cellDistance < *closestDistance) {
            closestIds->assign(1, cell.getId());
            *closestDistance = cellDistance;
        }
    }
}
 
#if BITPIT_ENABLE_MPI
bool SkdGlobalCellDistance::m_MPIDatatypeInitialized = false;
MPI_Datatype SkdGlobalCellDistance::m_MPIDatatype;
 
bool SkdGlobalCellDistance::m_MPIMinOperationInitialized = false;
MPI_Op SkdGlobalCellDistance::m_MPIMinOperation;
 
MPI_Datatype SkdGlobalCellDistance::getMPIDatatype()
{
    if (!m_MPIDatatypeInitialized) {
        int blocklengths[3] = {1, 1, 1};
        MPI_Aint displacements[3] = {offsetof(SkdGlobalCellDistance, m_distance), offsetof(SkdGlobalCellDistance, m_id), offsetof(SkdGlobalCellDistance, m_rank)};
        MPI_Datatype types[3] = {MPI_DOUBLE, MPI_LONG, MPI_INT};
        MPI_Type_create_struct(3, blocklengths, displacements, types, &m_MPIDatatype);
        MPI_Type_commit(&m_MPIDatatype);
        m_MPIDatatypeInitialized = true;
    }
 
    return m_MPIDatatype;
}
 
MPI_Op SkdGlobalCellDistance::getMPIMinOperation()
{
    if (!m_MPIMinOperationInitialized) {
        MPI_Op_create((MPI_User_function *) executeMPIMinOperation, false, &m_MPIMinOperation);
        m_MPIMinOperationInitialized = true;
    }
 
    return m_MPIMinOperation;
}
 
void SkdGlobalCellDistance::executeMPIMinOperation(SkdGlobalCellDistance *in, SkdGlobalCellDistance *inout, int *len, MPI_Datatype *datatype)
{
    BITPIT_UNUSED(datatype);
 
    int nDistances = *len;
    for (int i = 0; i < nDistances; ++i) {
        bool updateDistance = false;
        if (in[i].m_id != Cell::NULL_ID) {
            if (inout[i].m_id == Cell::NULL_ID) {
                updateDistance = true;
            } else if (std::abs(in[i].m_distance) < std::abs(inout[i].m_distance)) {
                updateDistance = true;
            }
        }
 
        if (updateDistance) {
            inout[i] = in[i];
        }
    }
}
 
int & SkdGlobalCellDistance::getRank()
{
    return m_rank;
}
 
long & SkdGlobalCellDistance::getId()
{
    return m_id;
}
 
double & SkdGlobalCellDistance::getDistance()
{
    return m_distance;
}
 
void SkdGlobalCellDistance::exportData(int *rank, long *id, double *distance) const
{
    *rank     = m_rank;
    *id       = m_id;
    *distance = m_distance;
}
#endif
 
PatchSkdTree::PatchSkdTree(const PatchKernel *patch, bool interiorCellsOnly)
    : m_patchInfo(patch, &m_cellRawIds),
      m_cellRawIds(interiorCellsOnly ? patch->getInternalCellCount() : patch->getCellCount()),
      m_nLeafs(0), m_nMinLeafCells(0), m_nMaxLeafCells(0),
      m_interiorCellsOnly(interiorCellsOnly),
      m_threadSafeLookups(false)
#if BITPIT_ENABLE_MPI
    , m_rank(0), m_nProcessors(1), m_communicator(MPI_COMM_NULL)
#endif
{
 
}
 
std::size_t PatchSkdTree::getLeafMinCellCount() const
{
    return m_nMinLeafCells;
}
 
std::size_t PatchSkdTree::getLeafMaxCellCount() const
{
    return m_nMaxLeafCells;
}
 
void PatchSkdTree::build(std::size_t leafThreshold, bool squeezeStorage)
{
    const PatchKernel &patch = m_patchInfo.getPatch();
 
    // Clear existing tree
    clear();
 
    // Initialize list of cell raw ids
    std::size_t nCells;
    PatchKernel::CellConstRange cellRange;
    if (m_interiorCellsOnly) {
        nCells = patch.getInternalCellCount();
        cellRange.initialize(patch.internalCellConstBegin(), patch.internalCellConstEnd());
    } else {
        nCells = patch.getCellCount();
        cellRange.initialize(patch.cellConstBegin(), patch.cellConstEnd());
    }
 
    if (m_cellRawIds.size() != nCells) {
        m_cellRawIds.resize(nCells);
    }
 
    std::size_t k = 0;
    for (auto itr = cellRange.begin(); itr != cellRange.end(); ++itr) {
        m_cellRawIds[k++] = itr.getRawIndex();
    }
 
    // Build patch cache
    m_patchInfo.buildCache(cellRange);
 
    // Initialize node list
    std::size_t nodesCount = std::max(1, int(std::ceil(2. * nCells / leafThreshold - 1.)));
    m_nodes.reserve(nodesCount);
 
    // Create the root
    m_nodes.emplace_back(&m_patchInfo, 0, nCells);
 
    // Create the tree
    std::queue<std::size_t> nodeStack;
    nodeStack.push(0);
    while (!nodeStack.empty()) {
        std::size_t nodeId = nodeStack.front();
        nodeStack.pop();
 
        // Create the children of the node
        //
        // This function may add new nodes to the tree, this may invalidate
        // any reference or pointer to the nodes. To avoid potential problems
        // we pass to the function a node id.
        createChildren(nodeId, leafThreshold);
 
        // Add the newly created childrend to the stack
        const SkdNode &node = getNode(nodeId);
        for (int i = SkdNode::CHILD_BEGIN; i != SkdNode::CHILD_END; ++i) {
            std::size_t childId = node.getChildId(static_cast<SkdNode::ChildLocation>(i));
            if (childId != SkdNode::NULL_ID) {
                nodeStack.push(childId);
            }
        }
    }
 
    // Squeeze storage
    if (squeezeStorage) {
        m_nodes.shrink_to_fit();
    }
 
    // Patch cache is no longer needed
    m_patchInfo.destroyCache();
 
#if BITPIT_ENABLE_MPI
    // Set partition information
    if (patch.isPartitioned()){
        // Set communicator
        setCommunicator(patch.getCommunicator());
 
        // Build partition info with partition boxes if the patch is partitioned
        buildPartitionBoxes();
    } else {
        m_rank        = 0;
        m_nProcessors = 1;
    }
#endif
}
 
void PatchSkdTree::clear(bool release)
{
    m_nLeafs        = 0;
    m_nMinLeafCells = 0;
    m_nMaxLeafCells = 0;
 
    if (release) {
        std::vector<SkdNode>().swap(m_nodes);
        std::vector<std::size_t>().swap(m_cellRawIds);
    } else {
        m_nodes.clear();
        m_cellRawIds.clear();
    }
 
#if BITPIT_ENABLE_MPI
    freeCommunicator();
    if (release) {
        std::vector<SkdBox>().swap(m_partitionBoxes);
    } else {
        m_partitionBoxes.clear();
    }
#endif
}
 
const PatchKernel & PatchSkdTree::getPatch() const
{
      return m_patchInfo.getPatch();
}
 
std::size_t PatchSkdTree::getNodeCount() const
{
      return m_nodes.size();
}
 
std::size_t PatchSkdTree::getLeafCount() const
{
      return m_nLeafs;
}
 
const SkdNode & PatchSkdTree::getNode(std::size_t nodeId) const
{
      return m_nodes[nodeId];
}
 
SkdNode & PatchSkdTree::_getNode(std::size_t nodeId)
{
      return m_nodes[nodeId];
}
 
std::size_t PatchSkdTree::evalMaxDepth(std::size_t rootId) const
{
    if (rootId == SkdNode::NULL_ID) {
        return 0;
    }
 
    const SkdNode &node = getNode(rootId);
 
    std::size_t depth = 0;
    for (int i = SkdNode::CHILD_BEGIN; i != SkdNode::CHILD_END; ++i) {
        SkdNode::ChildLocation childLocation = static_cast<SkdNode::ChildLocation>(i);
        depth = std::max(evalMaxDepth(node.getChildId(childLocation)), depth);
    }
    ++depth;
 
    return depth;
}
 
void PatchSkdTree::createChildren(std::size_t parentId, std::size_t leafThreshold)
{
    const SkdNode &parent = getNode(parentId);
 
    // Check if the parent is a leaf
    std::size_t parentCellCount = parent.getCellCount();
    if (parentCellCount <= leafThreshold) {
        createLeaf(parentId);
        return;
    }
 
    // Evaluate the preferred direction along which elements will be split.
    //
    // The elements will be split along a plane normal to the direction
    // for which the bounding box has the maximum length.
    const std::array<double, 3> &parentBoxMin = parent.getBoxMin();
    const std::array<double, 3> &parentBoxMax = parent.getBoxMax();
 
    int largerDirection = 0;
    double boxMaximumLength = parentBoxMax[largerDirection] - parentBoxMin[largerDirection];
    for (int d = 1; d < 3; ++d) {
        double length = parentBoxMax[d] - parentBoxMin[d];
        if (length > boxMaximumLength) {
            largerDirection  = d;
            boxMaximumLength = length;
        }
    }
 
    // Split the elements.
    //
    // Preferred direction is tried first, if the split along this direction
    // will produce and empty child (e.g., all the cell's box centroids have
    // the same coordinate along the split direction and therefore all the
    // cells would be assigned to the left chiled) the split will be performed
    // along one of the other directions.
    for (int d = 0; d < 3; ++d) {
        // Update the split direction
        int splitDirection = (largerDirection + d) % 3;
 
        // Get the threshold for the split
        double splitThreshold = parent.evalBoxWeightedMean(splitDirection);
 
        // Order the elements
        //
        // All the elements with a centroid coordinate less or equal than the
        // threshold will be assigned to the left child, the others will be
        // assigned to the right child.
        std::size_t leftBegin  = parent.m_cellRangeBegin;
        std::size_t leftEnd    = parent.m_cellRangeEnd;
        std::size_t rightBegin = parent.m_cellRangeBegin;
        std::size_t rightEnd   = parent.m_cellRangeEnd;
        while (true) {
            // Update the right begin
            while (rightBegin != leftEnd && m_patchInfo.evalCachedBoxMean(m_cellRawIds[rightBegin], splitDirection) <= splitThreshold) {
                rightBegin++;
            }
 
            // Update the left end
            while (rightBegin != leftEnd && m_patchInfo.evalCachedBoxMean(m_cellRawIds[leftEnd - 1], splitDirection) > splitThreshold) {
                leftEnd--;
            }
 
            // If all the elements are in the right position we can exit
            if (rightBegin == leftEnd) {
                break;
            }
 
            // If left end and and right begin are not equal, that the two ids
            // point to misplaced elements. Swap the elements, advance the ids
            // and continue iterating.
            std::iter_swap(m_cellRawIds.begin() + (leftEnd - 1), m_cellRawIds.begin() + rightBegin);
        }
 
        if (leftEnd <= leftBegin || rightEnd <= rightBegin) {
            continue;
        }
 
        // Create the left child
        //
        // Adding new nodes may invalidate any pointer and reference to the
        // nodes, we cannot store a reference to the parent node.
        long leftId = m_nodes.size();
        m_nodes.emplace_back(&m_patchInfo, leftBegin, leftEnd);
        _getNode(parentId).m_children[static_cast<std::size_t>(SkdNode::CHILD_LEFT)] = leftId;
 
        // Create the right child
        //
        // Adding new nodes may invalidate any pointer and reference to the
        // nodes, we cannot store a reference to the parent node.
        long rightId = m_nodes.size();
        m_nodes.emplace_back(&m_patchInfo, rightBegin, rightEnd);
        _getNode(parentId).m_children[static_cast<std::size_t>(SkdNode::CHILD_RIGHT)] = rightId;
 
        // Done
        return;
    }
 
    // It was not possible to split the elements. They have all the same
    // characteristic position and therefore they will be clustered together
    // in a leaf node.
    createLeaf(parentId);
}
 
void PatchSkdTree::createLeaf(std::size_t nodeId)
{
    const SkdNode &node = getNode(nodeId);
    std::size_t nodeCellCount = node.getCellCount();
 
    ++m_nLeafs;
    m_nMinLeafCells = std::min(nodeCellCount, m_nMinLeafCells);
    m_nMaxLeafCells = std::max(nodeCellCount, m_nMaxLeafCells);
}
 
void PatchSkdTree::enableThreadSafeLookups(bool enable)
{
    m_threadSafeLookups = enable;
}
 
bool PatchSkdTree::areLookupsThreadSafe() const
{
    return m_threadSafeLookups;
}
 
#if BITPIT_ENABLE_MPI
void PatchSkdTree::setCommunicator(MPI_Comm communicator)
{
    // Communication can be set just once
    if (isCommunicatorSet()) {
        throw std::runtime_error("Skd-tree communicator can be set just once.");
    }
 
    // The communicator has to be valid
    if (communicator == MPI_COMM_NULL) {
        throw std::runtime_error("Skd-tree communicator is not valid.");
    }
 
    // Creat a copy of the user-specified communicator
    //
    // No library routine should use MPI_COMM_WORLD as the communicator;
    // instead, a duplicate of a user-specified communicator should always
    // be used.
    MPI_Comm_dup(communicator, &m_communicator);
 
    // Get MPI information
    MPI_Comm_size(m_communicator, &m_nProcessors);
    MPI_Comm_rank(m_communicator, &m_rank);
}
 
void PatchSkdTree::freeCommunicator()
{
    if (!isCommunicatorSet()) {
        return;
    }
 
    int finalizedCalled;
    MPI_Finalized(&finalizedCalled);
    if (finalizedCalled) {
        return;
    }
 
    MPI_Comm_free(&m_communicator);
}
 
bool PatchSkdTree::isCommunicatorSet() const
{
    return (getCommunicator() != MPI_COMM_NULL);
}
 
const MPI_Comm & PatchSkdTree::getCommunicator() const
{
    return m_communicator;
}
 
void PatchSkdTree::buildPartitionBoxes()
{
    // Recover local bounding box of the root node
    const SkdBox &box = getNode(0).getBoundingBox();
 
    // Collect minimum and maximum coordinates of bounding boxes
    std::vector<int> count(m_nProcessors, 3);
 
    std::vector<int> displacements(m_nProcessors);
    for (int i = 0; i < m_nProcessors; ++i) {
        displacements[i] = i*3;
    }
 
    std::vector<std::array<double, 3>> minBoxes(m_nProcessors);
    minBoxes[m_rank] = box.getBoxMin();
    MPI_Allgatherv(MPI_IN_PLACE, 3, MPI_DOUBLE, minBoxes.data(),
                   count.data(), displacements.data(), MPI_DOUBLE, m_communicator);
 
    std::vector<std::array<double, 3>> maxBoxes(m_nProcessors);
    maxBoxes[m_rank] = box.getBoxMax();
    MPI_Allgatherv(MPI_IN_PLACE, 3, MPI_DOUBLE, maxBoxes.data(),
                    count.data(), displacements.data(), MPI_DOUBLE, m_communicator);
 
    // Fill partition boxes
    m_partitionBoxes.resize(getPatch().getProcessorCount());
    for (int i = 0; i < m_nProcessors; ++i) {
        m_partitionBoxes[i] = SkdBox(minBoxes[i], maxBoxes[i]);
    }
}
 
 
const SkdBox & PatchSkdTree::getPartitionBox(int rank) const
{
    return m_partitionBoxes[rank];
}
#endif
 
}