patch_kernel.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 <sstream>
#include <typeinfo>
#include <unordered_map>
#include <unordered_set>
#if BITPIT_ENABLE_MPI==1
#   include <mpi.h>
#endif
 
#include "bitpit_CG.hpp"
#include "bitpit_common.hpp"
 
#include "patch_info.hpp"
#include "patch_kernel.hpp"
#include "patch_manager.hpp"
 
namespace bitpit {

 
#if BITPIT_ENABLE_MPI==1
PatchKernel::PatchKernel(MPI_Comm communicator, std::size_t haloSize,
                         AdaptionMode adaptionMode, PartitioningMode partitioningMode)
#else
PatchKernel::PatchKernel(AdaptionMode adaptionMode)
#endif
    : m_adaptionMode(adaptionMode)
#if BITPIT_ENABLE_MPI==1
      , m_partitioningMode(partitioningMode)
#endif
{
    // Initialize the patch
#if BITPIT_ENABLE_MPI==1
    initialize(communicator, haloSize);
#else
    initialize();
#endif
 
    // Register the patch
    patch::manager().registerPatch(this);
}
 
#if BITPIT_ENABLE_MPI==1
PatchKernel::PatchKernel(int dimension, MPI_Comm communicator, std::size_t haloSize,
                         AdaptionMode adaptionMode, PartitioningMode partitioningMode)
#else
PatchKernel::PatchKernel(int dimension, AdaptionMode adaptionMode)
#endif
    : m_adaptionMode(adaptionMode)
#if BITPIT_ENABLE_MPI==1
      , m_partitioningMode(partitioningMode)
#endif
{
    // Initialize the patch
#if BITPIT_ENABLE_MPI==1
    initialize(communicator, haloSize);
#else
    initialize();
#endif
 
    // Register the patch
    patch::manager().registerPatch(this);
 
    // Set the dimension
    setDimension(dimension);
}
 
#if BITPIT_ENABLE_MPI==1
PatchKernel::PatchKernel(int id, int dimension, MPI_Comm communicator, std::size_t haloSize,
                         AdaptionMode adaptionMode, PartitioningMode partitioningMode)
#else
PatchKernel::PatchKernel(int id, int dimension, AdaptionMode adaptionMode)
#endif
    : m_adaptionMode(adaptionMode)
#if BITPIT_ENABLE_MPI==1
      , m_partitioningMode(partitioningMode)
#endif
{
    // Initialize the patch
#if BITPIT_ENABLE_MPI==1
    initialize(communicator, haloSize);
#else
    initialize();
#endif
 
    // Register the patch
    patch::manager().registerPatch(this, id);
 
    // Set the dimension
    //
    // Here we can only call the base function, if needed, every derived patch
    // has to call its derived function.
    PatchKernel::setDimension(dimension);
}

PatchKernel::PatchKernel(const PatchKernel &other)
    : VTKBaseStreamer(other),
      m_vertices(other.m_vertices),
      m_cells(other.m_cells),
      m_interfaces(other.m_interfaces),
      m_alteredCells(other.m_alteredCells),
      m_alteredInterfaces(other.m_alteredInterfaces),
      m_nInternalVertices(other.m_nInternalVertices),
#if BITPIT_ENABLE_MPI==1
      m_nGhostVertices(other.m_nGhostVertices),
#endif
      m_lastInternalVertexId(other.m_lastInternalVertexId),
#if BITPIT_ENABLE_MPI==1
      m_firstGhostVertexId(other.m_firstGhostVertexId),
#endif
      m_nInternalCells(other.m_nInternalCells),
#if BITPIT_ENABLE_MPI==1
      m_nGhostCells(other.m_nGhostCells),
#endif
      m_lastInternalCellId(other.m_lastInternalCellId),
#if BITPIT_ENABLE_MPI==1
      m_firstGhostCellId(other.m_firstGhostCellId),
#endif
      m_vtk(other.m_vtk),
      m_vtkWriteTarget(other.m_vtkWriteTarget),
      m_vtkVertexMap(other.m_vtkVertexMap),
      m_boxFrozen(other.m_boxFrozen),
      m_boxDirty(other.m_boxDirty),
      m_boxMinPoint(other.m_boxMinPoint),
      m_boxMaxPoint(other.m_boxMaxPoint),
      m_boxMinCounter(other.m_boxMinCounter),
      m_boxMaxCounter(other.m_boxMaxCounter),
      m_adjacenciesBuildStrategy(other.m_adjacenciesBuildStrategy),
      m_interfacesBuildStrategy(other.m_interfacesBuildStrategy),
      m_adaptionMode(other.m_adaptionMode),
      m_adaptionStatus(other.m_adaptionStatus),
      m_dimension(other.m_dimension),
      m_toleranceCustom(other.m_toleranceCustom),
      m_tolerance(other.m_tolerance)
#if BITPIT_ENABLE_MPI==1
      , m_partitioningMode(other.m_partitioningMode),
      m_partitioningStatus(other.m_partitioningStatus),
      m_owner(other.m_owner),
      m_haloSize(other.m_haloSize),
      m_partitioningCellsTag(other.m_partitioningCellsTag),
      m_partitioningVerticesTag(other.m_partitioningVerticesTag),
      m_partitioningSerialization(other.m_partitioningSerialization),
      m_partitioningOutgoings(other.m_partitioningOutgoings),
      m_partitioningGlobalExchanges(other.m_partitioningGlobalExchanges),
      m_partitioningInfoDirty(other.m_partitioningInfoDirty),
      m_ghostVertexInfo(other.m_ghostVertexInfo),
      m_ghostVertexExchangeTargets(other.m_ghostVertexExchangeTargets),
      m_ghostVertexExchangeSources(other.m_ghostVertexExchangeSources),
      m_ghostCellInfo(other.m_ghostCellInfo),
      m_ghostCellExchangeTargets(other.m_ghostCellExchangeTargets),
      m_ghostCellExchangeSources(other.m_ghostCellExchangeSources)
#endif
{
#if BITPIT_ENABLE_MPI==1
    // Initialize the communicator
    initializeCommunicator(other.getCommunicator());
#else
    // Initialize serial communicator
    initializeSerialCommunicator();
#endif
 
    // Create index generators
    importVertexIndexGenerator(other);
    importInterfaceIndexGenerator(other);
    importCellIndexGenerator(other);
 
    // Register the patch
    patch::manager().registerPatch(this);
 
    // Update the VTK streamer
    replaceVTKStreamer(&other, this);
}

PatchKernel::PatchKernel(PatchKernel &&other)
    : VTKBaseStreamer(std::move(other)),
      m_vertices(std::move(other.m_vertices)),
      m_cells(std::move(other.m_cells)),
      m_interfaces(std::move(other.m_interfaces)),
      m_alteredCells(std::move(other.m_alteredCells)),
      m_alteredInterfaces(std::move(other.m_alteredInterfaces)),
      m_vertexIdGenerator(std::move(other.m_vertexIdGenerator)),
      m_interfaceIdGenerator(std::move(other.m_interfaceIdGenerator)),
      m_cellIdGenerator(std::move(other.m_cellIdGenerator)),
      m_nInternalVertices(std::move(other.m_nInternalVertices)),
#if BITPIT_ENABLE_MPI==1
      m_nGhostVertices(std::move(other.m_nGhostVertices)),
#endif
      m_lastInternalVertexId(std::move(other.m_lastInternalVertexId)),
#if BITPIT_ENABLE_MPI==1
      m_firstGhostVertexId(std::move(other.m_firstGhostVertexId)),
#endif
      m_nInternalCells(std::move(other.m_nInternalCells)),
#if BITPIT_ENABLE_MPI==1
      m_nGhostCells(std::move(other.m_nGhostCells)),
#endif
      m_lastInternalCellId(std::move(other.m_lastInternalCellId)),
#if BITPIT_ENABLE_MPI==1
      m_firstGhostCellId(std::move(other.m_firstGhostCellId)),
#endif
      m_vtk(std::move(other.m_vtk)),
      m_vtkWriteTarget(std::move(other.m_vtkWriteTarget)),
      m_vtkVertexMap(std::move(other.m_vtkVertexMap)),
      m_boxFrozen(std::move(other.m_boxFrozen)),
      m_boxDirty(std::move(other.m_boxDirty)),
      m_boxMinPoint(std::move(other.m_boxMinPoint)),
      m_boxMaxPoint(std::move(other.m_boxMaxPoint)),
      m_boxMinCounter(std::move(other.m_boxMinCounter)),
      m_boxMaxCounter(std::move(other.m_boxMaxCounter)),
      m_adjacenciesBuildStrategy(std::move(other.m_adjacenciesBuildStrategy)),
      m_interfacesBuildStrategy(std::move(other.m_interfacesBuildStrategy)),
      m_adaptionMode(std::move(other.m_adaptionMode)),
      m_adaptionStatus(std::move(other.m_adaptionStatus)),
      m_id(std::move(other.m_id)),
      m_dimension(std::move(other.m_dimension)),
      m_toleranceCustom(std::move(other.m_toleranceCustom)),
      m_tolerance(std::move(other.m_tolerance)),
      m_rank(std::move(other.m_rank)),
      m_nProcessors(std::move(other.m_nProcessors))
#if BITPIT_ENABLE_MPI==1
      , m_communicator(std::move(MPI_COMM_NULL)),
      m_partitioningMode(other.m_partitioningMode),
      m_partitioningStatus(std::move(other.m_partitioningStatus)),
      m_owner(std::move(other.m_owner)),
      m_haloSize(std::move(other.m_haloSize)),
      m_partitioningCellsTag(std::move(other.m_partitioningCellsTag)),
      m_partitioningVerticesTag(std::move(other.m_partitioningVerticesTag)),
      m_partitioningSerialization(std::move(other.m_partitioningSerialization)),
      m_partitioningOutgoings(std::move(other.m_partitioningOutgoings)),
      m_partitioningGlobalExchanges(std::move(other.m_partitioningGlobalExchanges)),
      m_partitioningInfoDirty(std::move(other.m_partitioningInfoDirty)),
      m_ghostVertexInfo(std::move(other.m_ghostVertexInfo)),
      m_ghostVertexExchangeTargets(std::move(other.m_ghostVertexExchangeTargets)),
      m_ghostVertexExchangeSources(std::move(other.m_ghostVertexExchangeSources)),
      m_ghostCellInfo(std::move(other.m_ghostCellInfo)),
      m_ghostCellExchangeTargets(std::move(other.m_ghostCellExchangeTargets)),
      m_ghostCellExchangeSources(std::move(other.m_ghostCellExchangeSources))
#endif
{
    // Handle patch regstration
    patch::manager().unregisterPatch(&other);
 
    patch::manager().registerPatch(this, m_id);
    patch::manager().registerPatch(&other);
 
    // Update the VTK streamer
    //
    // Pointers to VTK streamers has been copied, we need to replace all the
    // pointer to the other object with pointer to this object.
    replaceVTKStreamer(&other, this);
 
#if BITPIT_ENABLE_MPI==1
    // Handle the communication
    std::swap(m_communicator, other.m_communicator);
#endif
}

PatchKernel & PatchKernel::operator=(PatchKernel &&other)
{
    VTKBaseStreamer::operator=(std::move(other));
    m_vertices = std::move(other.m_vertices);
    m_cells = std::move(other.m_cells);
    m_interfaces = std::move(other.m_interfaces);
    m_alteredCells = std::move(other.m_alteredCells);
    m_alteredInterfaces = std::move(other.m_alteredInterfaces);
    m_vertexIdGenerator = std::move(other.m_vertexIdGenerator);
    m_interfaceIdGenerator = std::move(other.m_interfaceIdGenerator);
    m_cellIdGenerator = std::move(other.m_cellIdGenerator);
    m_nInternalVertices = std::move(other.m_nInternalVertices);
#if BITPIT_ENABLE_MPI==1
    m_nGhostVertices = std::move(other.m_nGhostVertices);
#endif
    m_lastInternalVertexId = std::move(other.m_lastInternalVertexId);
#if BITPIT_ENABLE_MPI==1
    m_firstGhostVertexId = std::move(other.m_firstGhostVertexId);
#endif
    m_nInternalCells = std::move(other.m_nInternalCells);
#if BITPIT_ENABLE_MPI==1
    m_nGhostCells = std::move(other.m_nGhostCells);
#endif
    m_lastInternalCellId = std::move(other.m_lastInternalCellId);
#if BITPIT_ENABLE_MPI==1
    m_firstGhostCellId = std::move(other.m_firstGhostCellId);
#endif
    m_vtk = std::move(other.m_vtk);
    m_vtkWriteTarget = std::move(other.m_vtkWriteTarget);
    m_vtkVertexMap = std::move(other.m_vtkVertexMap);
    m_boxFrozen = std::move(other.m_boxFrozen);
    m_boxDirty = std::move(other.m_boxDirty);
    m_boxMinPoint = std::move(other.m_boxMinPoint);
    m_boxMaxPoint = std::move(other.m_boxMaxPoint);
    m_boxMinCounter = std::move(other.m_boxMinCounter);
    m_boxMaxCounter = std::move(other.m_boxMaxCounter);
    m_adjacenciesBuildStrategy = std::move(other.m_adjacenciesBuildStrategy);
    m_interfacesBuildStrategy = std::move(other.m_interfacesBuildStrategy);
    m_adaptionMode = std::move(other.m_adaptionMode);
    m_adaptionStatus = std::move(other.m_adaptionStatus);
    m_id = std::move(other.m_id);
    m_dimension = std::move(other.m_dimension);
    m_toleranceCustom = std::move(other.m_toleranceCustom);
    m_tolerance = std::move(other.m_tolerance);
    m_rank = std::move(other.m_rank);
    m_nProcessors = std::move(other.m_nProcessors);
#if BITPIT_ENABLE_MPI==1
    m_communicator = std::move(MPI_COMM_NULL);
    m_partitioningMode = std::move(other.m_partitioningMode);
    m_partitioningStatus = std::move(other.m_partitioningStatus);
    m_owner = std::move(other.m_owner);
    m_haloSize = std::move(other.m_haloSize);
    m_partitioningCellsTag = std::move(other.m_partitioningCellsTag);
    m_partitioningVerticesTag = std::move(other.m_partitioningVerticesTag);
    m_partitioningSerialization = std::move(other.m_partitioningSerialization);
    m_partitioningOutgoings = std::move(other.m_partitioningOutgoings);
    m_partitioningGlobalExchanges = std::move(other.m_partitioningGlobalExchanges);
    m_partitioningInfoDirty = std::move(other.m_partitioningInfoDirty);
    m_ghostVertexInfo = std::move(other.m_ghostVertexInfo);
    m_ghostVertexExchangeTargets = std::move(other.m_ghostVertexExchangeTargets);
    m_ghostVertexExchangeSources = std::move(other.m_ghostVertexExchangeSources);
    m_ghostCellInfo = std::move(other.m_ghostCellInfo);
    m_ghostCellExchangeTargets = std::move(other.m_ghostCellExchangeTargets);
    m_ghostCellExchangeSources = std::move(other.m_ghostCellExchangeSources);
#endif
 
    // Handle patch regstration
    patch::manager().unregisterPatch(this);
    patch::manager().unregisterPatch(&other);
 
    patch::manager().registerPatch(this, m_id);
    patch::manager().registerPatch(&other);
 
    // Update the VTK streamer
    //
    // Pointers to VTK streamers has been moved, we need to replace all the
    // pointer to the other object with pointer to this object.
    replaceVTKStreamer(&other, this);
 
#if BITPIT_ENABLE_MPI==1
    // Handle the communication
    std::swap(m_communicator, other.m_communicator);
#endif
 
    return *this;
}

#if BITPIT_ENABLE_MPI==1
void PatchKernel::initialize(MPI_Comm communicator, std::size_t haloSize)
#else
void PatchKernel::initialize()
#endif
{
    // Id
    m_id = PatchManager::AUTOMATIC_ID;
 
    // Vertex count
    m_nInternalVertices = 0;
#if BITPIT_ENABLE_MPI==1
    m_nGhostVertices = 0;
#endif
 
    m_lastInternalVertexId = Vertex::NULL_ID;
#if BITPIT_ENABLE_MPI==1
    m_firstGhostVertexId = Vertex::NULL_ID;
#endif
 
    // Cell count
    m_nInternalCells = 0;
#if BITPIT_ENABLE_MPI==1
    m_nGhostCells    = 0;
#endif
 
    m_lastInternalCellId = Cell::NULL_ID;
#if BITPIT_ENABLE_MPI==1
    m_firstGhostCellId   = Cell::NULL_ID;
#endif
 
    // Dimension
    m_dimension = -1;
 
    // Index generators
    setVertexAutoIndexing(true);
    setInterfaceAutoIndexing(true);
    setCellAutoIndexing(true);
 
    // Set adjacencies build strategy
    setAdjacenciesBuildStrategy(ADJACENCIES_NONE);
 
    // Set interfaces build strategy
    setInterfacesBuildStrategy(INTERFACES_NONE);
 
    // Set the adaption as clean
    setAdaptionStatus(ADAPTION_CLEAN);
 
#if BITPIT_ENABLE_MPI==1
    // Initialize communicator
    initializeCommunicator(communicator);
 
    // Set halo size
    initializeHaloSize(haloSize);
 
    // Mark patch as partioned
    setPartitioningStatus(PARTITIONING_CLEAN);
 
    // Initialize partitioning tags
    m_partitioningCellsTag    = -1;
    m_partitioningVerticesTag = -1;
 
    // Update partitioning information
    if (isPartitioned()) {
        updatePartitioningInfo(true);
    }
#else
    // Initialize serial communicator
    initializeSerialCommunicator();
#endif
 
    // Initialize the geometrical tolerance
    resetTol();
 
    // Initializes the bounding box
    setBoundingBoxFrozen(false);
    clearBoundingBox();
 
    // Set VTK write target
    m_vtkWriteTarget = WRITE_TARGET_CELLS_ALL;
 
    // Set VTK information
    std::ostringstream convert;
    convert << getId();
 
    m_vtk.setName(convert.str());
    m_vtk.setCodex(VTKFormat::APPENDED);
 
    // Set VTK Geom Data
    m_vtk.setGeomData<double>(VTKUnstructuredField::POINTS, this);
    m_vtk.setGeomData<long>(VTKUnstructuredField::OFFSETS, this);
    m_vtk.setGeomData<int>(VTKUnstructuredField::TYPES, this);
    m_vtk.setGeomData<long>(VTKUnstructuredField::CONNECTIVITY, this);
    m_vtk.setGeomData<long>(VTKUnstructuredField::FACE_STREAMS, this);
    m_vtk.setGeomData<long>(VTKUnstructuredField::FACE_OFFSETS, this);
 
    // Add VTK basic patch data
    m_vtk.addData<long>("cellIndex", VTKFieldType::SCALAR, VTKLocation::CELL, this);
    m_vtk.addData<int>("PID", VTKFieldType::SCALAR, VTKLocation::CELL, this);
    m_vtk.addData<long>("vertexIndex", VTKFieldType::SCALAR, VTKLocation::POINT, this);
#if BITPIT_ENABLE_MPI==1
    m_vtk.addData<long>("cellGlobalIndex", VTKFieldType::SCALAR, VTKLocation::CELL, this);
    m_vtk.addData<int>("cellRank", VTKFieldType::SCALAR, VTKLocation::CELL, this);
    m_vtk.addData<int>("cellHaloLayer", VTKFieldType::SCALAR, VTKLocation::CELL, this);
    m_vtk.addData<int>("vertexRank", VTKFieldType::SCALAR, VTKLocation::POINT, this);
#endif
}
 
#if BITPIT_ENABLE_MPI!=1
void PatchKernel::initializeSerialCommunicator()
{
    m_rank        = 0;
    m_nProcessors = 1;
}
#endif

PatchKernel::~PatchKernel()
{
#if BITPIT_ENABLE_MPI==1
    freeCommunicator();
#endif
 
    try {
        patch::manager().unregisterPatch(this);
    } catch (const std::runtime_error &e) {
        log::cout() << "Unable to unregister the patch" << std::endl;
        log::cout() << " Error message: " << e.what() << std::endl;
    }
}

std::vector<adaption::Info> PatchKernel::update(bool trackAdaption, bool squeezeStorage)
{
    std::vector<adaption::Info> adaptionData;
 
    // Early return if the patch is not dirty
    //
    // If we need to squeeze the storage we need to perform the update also
    // if the patch is not dirty.
    if (!squeezeStorage && !isDirty(true)) {
        return adaptionData;
    }
 
    // Finalize alterations
    finalizeAlterations(squeezeStorage);
 
    // Adaption
    bool adaptionDirty = (getAdaptionStatus(true) == ADAPTION_DIRTY);
    if (adaptionDirty) {
        mergeAdaptionInfo(adaption(trackAdaption, squeezeStorage), adaptionData);
    }
 
    return adaptionData;
}

void PatchKernel::simulateCellUpdate(const long id, adaption::Marker marker, std::vector<Cell> *virtualCells, PiercedVector<Vertex, long> *virtualVertices) const
{
    BITPIT_UNUSED(id);
    BITPIT_UNUSED(marker);
    BITPIT_UNUSED(virtualCells);
    BITPIT_UNUSED(virtualVertices);
 
    throw std::runtime_error ("This function has not been implemented for the specified patch.");
}

std::vector<adaption::Info> PatchKernel::adaption(bool trackAdaption, bool squeezeStorage)
{
    std::vector<adaption::Info> adaptionData;
 
    // Early return if adaption cannot be performed
    AdaptionMode adaptionMode = getAdaptionMode();
    if (adaptionMode == ADAPTION_DISABLED) {
        return adaptionData;
    }
 
    AdaptionStatus adaptionStatus = getAdaptionStatus(true);
    if (adaptionStatus == ADAPTION_CLEAN) {
        return adaptionData;
    } else if (adaptionStatus != ADAPTION_DIRTY) {
        throw std::runtime_error ("An adaption is already in progress.");
    }
 
    // Run adaption
    adaptionPrepare(false);
 
    adaptionData = adaptionAlter(trackAdaption, squeezeStorage);
 
    adaptionCleanup();
 
    return adaptionData;
}

std::vector<adaption::Info> PatchKernel::adaptionPrepare(bool trackAdaption)
{
    std::vector<adaption::Info> adaptionData;
 
    // Early return if adaption cannot be performed
    AdaptionMode adaptionMode = getAdaptionMode();
    if (adaptionMode == ADAPTION_DISABLED) {
        return adaptionData;
    }
 
    AdaptionStatus adaptionStatus = getAdaptionStatus(true);
    if (adaptionStatus == ADAPTION_CLEAN) {
        return adaptionData;
    } else if (adaptionStatus != ADAPTION_DIRTY) {
        throw std::runtime_error ("An adaption is already in progress.");
    }
 
    // Execute the adaption preparation
    adaptionData = _adaptionPrepare(trackAdaption);
 
    // Update the status
    setAdaptionStatus(ADAPTION_PREPARED);
 
    return adaptionData;
}

std::vector<adaption::Info> PatchKernel::adaptionAlter(bool trackAdaption, bool squeezeStorage)
{
    std::vector<adaption::Info> adaptionData;
 
    // Early return if adaption cannot be performed
    AdaptionMode adaptionMode = getAdaptionMode();
    if (adaptionMode == ADAPTION_DISABLED) {
        return adaptionData;
    }
 
    AdaptionStatus adaptionStatus = getAdaptionStatus();
    if (adaptionStatus == ADAPTION_CLEAN) {
        return adaptionData;
    } else if (adaptionStatus != ADAPTION_PREPARED) {
        throw std::runtime_error ("The prepare function has not been called.");
    }
 
    // Adapt the patch
    adaptionData = _adaptionAlter(trackAdaption);
 
    // Finalize patch alterations
    finalizeAlterations(squeezeStorage);
 
    // Update the status
    setAdaptionStatus(ADAPTION_ALTERED);
 
    return adaptionData;
}

void PatchKernel::adaptionCleanup()
{
    // Early return if adaption cannot be performed
    AdaptionMode adaptionMode = getAdaptionMode();
    if (adaptionMode == ADAPTION_DISABLED) {
        return;
    }
 
    AdaptionStatus adaptionStatus = getAdaptionStatus();
    if (adaptionStatus == ADAPTION_CLEAN) {
        return;
    } else if (adaptionStatus == ADAPTION_PREPARED) {
        throw std::runtime_error ("It is not yet possible to abort an adaption.");
    } else if (adaptionStatus != ADAPTION_ALTERED) {
        throw std::runtime_error ("The alter function has not been called.");
    }
 
    // Complete the adaption
    _adaptionCleanup();
 
    // Update the status
    setAdaptionStatus(ADAPTION_CLEAN);
}

void PatchKernel::settleAdaptionMarkers()
{
    // Nothing to do
}

void PatchKernel::finalizeAlterations(bool squeezeStorage)
{
    // Flush vertex data structures
    m_vertices.flush();
 
    // Update bounding box
    bool boundingBoxDirty = isBoundingBoxDirty();
    if (boundingBoxDirty) {
        updateBoundingBox();
    }
 
    // Flush cell data structures
    m_cells.flush();
 
    // Update adjacencies
    bool adjacenciesDirty = areAdjacenciesDirty();
    if (adjacenciesDirty) {
        updateAdjacencies();
    }
 
    // Update interfaces
    bool interfacesDirty = areInterfacesDirty();
    if (interfacesDirty) {
        updateInterfaces();
    }
 
    // Flush interfaces data structures
    m_interfaces.flush();
 
#if BITPIT_ENABLE_MPI==1
    // Update partitioning information
    bool partitioningInfoDirty = arePartitioningInfoDirty();
    if (partitioningInfoDirty) {
        updatePartitioningInfo(true);
    }
#endif
 
    // Clear alteration flags
    m_alteredCells.clear();
    m_alteredInterfaces.clear();
 
    // Squeeze the patch
    if (squeezeStorage) {
        squeeze();
    }
 
    // Synchronize storage
    m_cells.sync();
    m_interfaces.sync();
    m_vertices.sync();
}

void PatchKernel::markCellForRefinement(long id)
{
    bool updated = _markCellForRefinement(id);
 
    if (updated) {
        setAdaptionStatus(ADAPTION_DIRTY);
    }
}

void PatchKernel::markCellForCoarsening(long id)
{
    bool updated = _markCellForCoarsening(id);
 
    if (updated) {
        setAdaptionStatus(ADAPTION_DIRTY);
    }
}

void PatchKernel::resetCellAdaptionMarker(long id)
{
    bool updated = _resetCellAdaptionMarker(id);
 
    if (updated) {
        setAdaptionStatus(ADAPTION_DIRTY);
    }
}

adaption::Marker PatchKernel::getCellAdaptionMarker(long id)
{
    return _getCellAdaptionMarker(id);
}

void PatchKernel::enableCellBalancing(long id, bool enabled)
{
    bool updated = _enableCellBalancing(id, enabled);
 
    if (updated) {
        setAdaptionStatus(ADAPTION_DIRTY);
    }
}

void PatchKernel::reset()
{
    resetVertices();
    resetCells();
    resetInterfaces();
}

void PatchKernel::resetVertices()
{
    m_vertices.clear();
    if (m_vertexIdGenerator) {
        m_vertexIdGenerator->reset();
    }
    m_nInternalVertices = 0;
#if BITPIT_ENABLE_MPI==1
    m_nGhostVertices = 0;
#endif
    m_lastInternalVertexId = Vertex::NULL_ID;
#if BITPIT_ENABLE_MPI==1
    m_firstGhostVertexId = Vertex::NULL_ID;
#endif
 
    for (auto &cell : m_cells) {
        cell.unsetConnect();
    }
}

void PatchKernel::resetCells()
{
    m_cells.clear();
    if (m_cellIdGenerator) {
        m_cellIdGenerator->reset();
    }
    m_nInternalCells = 0;
#if BITPIT_ENABLE_MPI==1
    m_nGhostCells = 0;
#endif
    m_lastInternalCellId = Cell::NULL_ID;
#if BITPIT_ENABLE_MPI==1
    m_firstGhostCellId = Cell::NULL_ID;
#endif
 
#if BITPIT_ENABLE_MPI==1
    clearGhostCellsInfo();
    clearGhostVerticesInfo();
#endif
 
    resetAdjacencies();
 
    for (auto &interface : m_interfaces) {
        interface.unsetNeigh();
        interface.unsetOwner();
    }
 
    m_alteredCells.clear();
}

void PatchKernel::resetInterfaces()
{
    // Early return if no interfaces have been built
    if (getInterfacesBuildStrategy() == INTERFACES_NONE) {
        return;
    }
 
    // Reset the interfaces
    _resetInterfaces(false);
 
    // Mark cell interfaces as dirty
    setCellAlterationFlags(FLAG_INTERFACES_DIRTY);
 
    // Clear list of altered interfaces
    m_alteredInterfaces.clear();
}

void PatchKernel::_resetInterfaces(bool release)
{
    for (auto &cell : m_cells) {
        cell.resetInterfaces(!release);
    }
 
    m_interfaces.clear(release);
    if (m_interfaceIdGenerator) {
        m_interfaceIdGenerator->reset();
    }
}

bool PatchKernel::reserveVertices(size_t nVertices)
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    m_vertices.reserve(nVertices);
 
    return true;
}

bool PatchKernel::reserveCells(size_t nCells)
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    m_cells.reserve(nCells);
 
    return true;
}

bool PatchKernel::reserveInterfaces(size_t nInterfaces)
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    m_interfaces.reserve(nInterfaces);
 
    return true;
}

void PatchKernel::write(const std::string &filename, VTKWriteMode mode)
{
    std::string oldFilename = m_vtk.getName();
 
    m_vtk.setName(filename);
    write(mode);
    m_vtk.setName(oldFilename);
}

void PatchKernel::write(const std::string &filename, VTKWriteMode mode, double time)
{
    std::string oldFilename = m_vtk.getName();
 
    m_vtk.setName(filename);
    write(mode, time);
    m_vtk.setName(oldFilename);
}

void PatchKernel::write(VTKWriteMode mode)
{
    _writePrepare();
 
    m_vtk.write(mode);
 
    _writeFinalize();
}

void PatchKernel::write(VTKWriteMode mode, double time)
{
    _writePrepare();
 
    m_vtk.write(mode, time);
 
    _writeFinalize();
}

void PatchKernel::_writePrepare()
{
    // Get VTK cell count
    long vtkCellCount = 0;
    if (m_vtkWriteTarget == WRITE_TARGET_CELLS_ALL) {
        vtkCellCount = getCellCount();
#if BITPIT_ENABLE_MPI==1
    } else if (m_vtkWriteTarget == WRITE_TARGET_CELLS_INTERNAL) {
        vtkCellCount = getInternalCellCount();
#endif
    }
 
    // Set the dimensions of the mesh
    //
    // Even if there is just a single process that needs to write the VTK face stream, than
    // all the processes need to write the face stream as well (even the processes whose
    // local cells don't require the face steam).
    PiercedStorage<bool, long> vertexWriteFlag(1, &m_vertices);
    vertexWriteFlag.fill(false);
 
    bool vtkFaceStreamNeeded = false;
    for (const Cell &cell : m_cells) {
        if (cell.getDimension() > 2 && !cell.hasInfo()) {
            vtkFaceStreamNeeded = true;
            break;
        }
    }
 
#if BITPIT_ENABLE_MPI==1
    if (isPartitioned()) {
        const auto &communicator = getCommunicator();
        MPI_Allreduce(MPI_IN_PLACE, &vtkFaceStreamNeeded, 1, MPI_C_BOOL, MPI_LOR, communicator);
    }
#endif
 
    long vtkConnectSize    = 0;
    long vtkFaceStreamSize = 0;
    for (const Cell &cell : getVTKCellWriteRange()) {
        ConstProxyVector<long> cellVertexIds = cell.getVertexIds();
        const int nCellVertices = cellVertexIds.size();
        for (int k = 0; k < nCellVertices; ++k) {
            long vertexId = cellVertexIds[k];
            vertexWriteFlag.at(vertexId) = true;
        }
 
        vtkConnectSize += nCellVertices;
        if (vtkFaceStreamNeeded) {
            if (cell.getDimension() <= 2 || cell.hasInfo()) {
                vtkFaceStreamSize += 1;
            } else {
                vtkFaceStreamSize += cell.getFaceStreamSize();
            }
        }
    }
 
    int vtkVertexCount = 0;
    m_vtkVertexMap.unsetKernel();
    m_vtkVertexMap.setStaticKernel(&m_vertices);
    VertexConstIterator endItr = vertexConstEnd();
    for (VertexConstIterator itr = vertexConstBegin(); itr != endItr; ++itr) {
        std::size_t vertexRawId = itr.getRawIndex();
        if (vertexWriteFlag.rawAt(vertexRawId)) {
            m_vtkVertexMap.rawAt(vertexRawId) = vtkVertexCount++;
        } else {
            m_vtkVertexMap.rawAt(vertexRawId) = Vertex::NULL_ID;
        }
    }
 
    m_vtk.setDimensions(vtkCellCount, vtkVertexCount, vtkConnectSize, vtkFaceStreamSize);
}

void PatchKernel::_writeFinalize()
{
    // Nothing to do
}

bool PatchKernel::isAdaptionSupported() const
{
    return (getAdaptionMode() != ADAPTION_DISABLED);
}

PatchKernel::AdaptionMode PatchKernel::getAdaptionMode() const
{
    return m_adaptionMode;
}

void PatchKernel::setAdaptionMode(AdaptionMode mode)
{
    m_adaptionMode = mode;
}

PatchKernel::AdaptionStatus PatchKernel::getAdaptionStatus(bool global) const
{
    int adaptionStatus = static_cast<int>(m_adaptionStatus);
 
#if BITPIT_ENABLE_MPI==1
    if (global && isPartitioned()) {
        const auto &communicator = getCommunicator();
        MPI_Allreduce(MPI_IN_PLACE, &adaptionStatus, 1, MPI_INT, MPI_MAX, communicator);
    }
#else
    BITPIT_UNUSED(global);
#endif
 
    return static_cast<AdaptionStatus>(adaptionStatus);
}

void PatchKernel::setAdaptionStatus(AdaptionStatus status)
{
    m_adaptionStatus = status;
}

bool PatchKernel::isDirty(bool global) const
{
#if BITPIT_ENABLE_MPI==0
    BITPIT_UNUSED(global);
#endif
 
    bool isDirty = false;
 
    if (!isDirty) {
        isDirty |= areAdjacenciesDirty(false);
    }
 
    if (!isDirty) {
        isDirty |= !m_alteredCells.empty();
    }
 
    if (!isDirty) {
        isDirty |= areInterfacesDirty(false);
        assert(isDirty || m_alteredInterfaces.empty());
    }
 
    if (!isDirty) {
        isDirty |= (getAdaptionStatus(false) == ADAPTION_DIRTY);
    }
 
    if (!isDirty) {
        isDirty |= isBoundingBoxDirty(false);
    }
 
#if BITPIT_ENABLE_MPI==1
    if (!isDirty) {
        isDirty |= arePartitioningInfoDirty(false);
    }
#endif
 
    if (!isDirty) {
        isDirty |= !m_vertices.isSynced();
    }
 
    if (!isDirty) {
        isDirty |= !m_cells.isSynced();
    }
 
    if (!isDirty) {
        isDirty |= !m_interfaces.isSynced();
    }
 
#if BITPIT_ENABLE_MPI==1
    // Get the global flag
    if (global && isPartitioned()) {
        const auto &communicator = getCommunicator();
        MPI_Allreduce(MPI_IN_PLACE, &isDirty, 1, MPI_C_BOOL, MPI_LOR, communicator);
    }
#endif
 
    return isDirty;
}

void PatchKernel::setExpert(bool expert)
{
    static AdaptionMode initialAdaptionMode;
    if (!expert) {
        initialAdaptionMode = m_adaptionMode;
        m_adaptionMode = ADAPTION_MANUAL;
    } else {
        m_adaptionMode = initialAdaptionMode;
    }
}

bool PatchKernel::isExpert() const
{
    return (m_adaptionMode == ADAPTION_MANUAL);
}

void PatchKernel::setId(int id)
{
    if (id == m_id) {
        return;
    }
 
    patch::manager().unregisterPatch(this);
    patch::manager().registerPatch(this, id);
}

void PatchKernel::_setId(int id)
{
    m_id = id;
}

int PatchKernel::getId() const
{
    return m_id;
}

void PatchKernel::setDimension(int dimension)
{
    if (dimension == m_dimension) {
        return;
    }
 
    // If the dimension was already assigned, reset the patch
    if (m_dimension > 0) {
        reset();
    }
 
    // Set the dimension
    m_dimension = dimension;
}

int PatchKernel::getDimension() const
{
    return m_dimension;
}

bool PatchKernel::isThreeDimensional() const
{
    return (m_dimension == 3);
}

bool PatchKernel::isVertexAutoIndexingEnabled() const
{
    return static_cast<bool>(m_vertexIdGenerator);
}

void PatchKernel::setVertexAutoIndexing(bool enabled)
{
    if (isVertexAutoIndexingEnabled() == enabled) {
        return;
    }
 
    if (enabled) {
        createVertexIndexGenerator(true);
    } else {
        m_vertexIdGenerator.reset();
    }
}

void PatchKernel::dumpVertexAutoIndexing(std::ostream &stream) const
{
    m_vertexIdGenerator->dump(stream);
}

void PatchKernel::restoreVertexAutoIndexing(std::istream &stream)
{
    createVertexIndexGenerator(false);
    m_vertexIdGenerator->restore(stream);
}

void PatchKernel::createVertexIndexGenerator(bool populate)
{
    // Create or reset index generator
    if (!m_vertexIdGenerator) {
        m_vertexIdGenerator = std::unique_ptr<IndexGenerator<long>>(new IndexGenerator<long>());
    } else {
        m_vertexIdGenerator->reset();
    }
 
    // Populate index generator with current ids
    if (populate) {
        VertexConstIterator beginItr = m_vertices.cbegin();
        VertexConstIterator endItr   = m_vertices.cend();
        for (VertexConstIterator itr = beginItr; itr != endItr; ++itr) {
            m_vertexIdGenerator->setAssigned(itr.getId());
        }
    }
}

void PatchKernel::importVertexIndexGenerator(const PatchKernel &source)
{
    if (source.m_vertexIdGenerator) {
        m_vertexIdGenerator = std::unique_ptr<IndexGenerator<long>>(new IndexGenerator<long>(*(source.m_vertexIdGenerator)));
    } else {
        m_vertexIdGenerator.reset();
    }
}

long PatchKernel::getVertexCount() const
{
    return m_vertices.size();
}

long PatchKernel::getInternalVertexCount() const
{
    return m_nInternalVertices;
}

PiercedVector<Vertex> & PatchKernel::getVertices()
{
    return m_vertices;
}

const PiercedVector<Vertex> & PatchKernel::getVertices() const
{
    return m_vertices;
}

Vertex & PatchKernel::getVertex(long id)
{
    return m_vertices[id];
}

const Vertex & PatchKernel::getVertex(long id) const
{
    return m_vertices[id];
}

Vertex & PatchKernel::getLastInternalVertex()
{
    return m_vertices[m_lastInternalVertexId];
}

const Vertex & PatchKernel::getLastInternalVertex() const
{
    return m_vertices[m_lastInternalVertexId];
}

PatchKernel::VertexIterator PatchKernel::getVertexIterator(long id)
{
    return m_vertices.find(id);
}

PatchKernel::VertexIterator PatchKernel::vertexBegin()
{
    return m_vertices.begin();
}

PatchKernel::VertexIterator PatchKernel::vertexEnd()
{
    return m_vertices.end();
}

PatchKernel::VertexIterator PatchKernel::internalVertexBegin()
{
    return m_vertices.begin();
}

PatchKernel::VertexIterator PatchKernel::internalVertexEnd()
{
    if (m_nInternalVertices > 0) {
        return ++m_vertices.find(m_lastInternalVertexId);
    } else {
        return m_vertices.end();
    }
}

PatchKernel::VertexConstIterator PatchKernel::getVertexConstIterator(long id) const
{
    return m_vertices.find(id);
}

PatchKernel::VertexConstIterator PatchKernel::vertexConstBegin() const
{
    return m_vertices.cbegin();
}

PatchKernel::VertexConstIterator PatchKernel::vertexConstEnd() const
{
    return m_vertices.cend();
}

PatchKernel::VertexConstIterator PatchKernel::internalVertexConstBegin() const
{
    return m_vertices.cbegin();
}

PatchKernel::VertexConstIterator PatchKernel::internalVertexConstEnd() const
{
    if (m_nInternalVertices > 0) {
        return ++m_vertices.find(m_lastInternalVertexId);
    } else {
        return m_vertices.end();
    }
}

PatchKernel::VertexIterator PatchKernel::addVertex(const Vertex &source, long id)
{
    Vertex vertex = source;
    vertex.setId(id);
 
    return addVertex(std::move(vertex), id);
}

PatchKernel::VertexIterator PatchKernel::addVertex(Vertex &&source, long id)
{
    // Get the id
    if (id < 0) {
        id = source.getId();
    }
 
    // Add a dummy vertex
    //
    // It is not possible to directly add the source into the storage. First a
    // dummy vertex is created and then that vertex is replaced with the source.
    //
    // The corrdinates of the dummy vertex shuold match the coordinates of the
    // source, because the bounding box of the patch is updated every time a
    // new vertex is added.
    VertexIterator iterator = _addInternalVertex(source.getCoords(), id);
 
    // Replace the newly created vertex with the source
    //
    // Before replacing the newly created vertex we need to set the id
    // of the source to the id that has been assigned to the newly created
    // vertex, we also need to mark the vertex as internal.
    source.setId(iterator->getId());
    source.setInterior(true);
 
    Vertex &vertex = (*iterator);
    vertex = std::move(source);
 
    return iterator;
}

PatchKernel::VertexIterator PatchKernel::addVertex(const std::array<double, 3> &coords, long id)
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return vertexEnd();
    }
 
    // Add the vertex
    VertexIterator iterator = _addInternalVertex(coords, id);
 
    return iterator;
}

PatchKernel::VertexIterator PatchKernel::_addInternalVertex(const std::array<double, 3> &coords, long id)
{
    // Get the id
    if (m_vertexIdGenerator) {
        if (id < 0) {
            id = m_vertexIdGenerator->generate();
        } else {
            m_vertexIdGenerator->setAssigned(id);
        }
    } else if (id < 0) {
        throw std::runtime_error("No valid id has been provided for the vertex.");
    }
 
    // Get the id of the vertex before which the new vertex should be inserted
#if BITPIT_ENABLE_MPI==1
    //
    // If there are ghosts vertices, the internal vertex should be inserted
    // before the first ghost vertex.
#endif
    long referenceId;
#if BITPIT_ENABLE_MPI==1
    referenceId = m_firstGhostVertexId;
#else
    referenceId = Vertex::NULL_ID;
#endif
 
    // Create the vertex
    VertexIterator iterator;
    if (referenceId == Vertex::NULL_ID) {
        iterator = m_vertices.emreclaim(id, id, coords, true);
    } else {
        iterator = m_vertices.emreclaimBefore(referenceId, id, id, coords, true);
    }
    m_nInternalVertices++;
 
    // Update the id of the last internal vertex
    if (m_lastInternalVertexId < 0) {
        m_lastInternalVertexId = id;
    } else if (m_vertices.rawIndex(m_lastInternalVertexId) < m_vertices.rawIndex(id)) {
        m_lastInternalVertexId = id;
    }
 
    // Update the bounding box
    addPointToBoundingBox(iterator->getCoords());
 
#if BITPIT_ENABLE_MPI==1
    // Set partitioning information as dirty
    setPartitioningInfoDirty(true);
#endif
 
    return iterator;
}
 
#if BITPIT_ENABLE_MPI==0
PatchKernel::VertexIterator PatchKernel::restoreVertex(const std::array<double, 3> &coords, long id)
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return vertexEnd();
    }
 
    VertexIterator iterator = m_vertices.find(id);
    if (iterator == m_vertices.end()) {
        throw std::runtime_error("Unable to restore the specified vertex: the kernel doesn't contain an entry for that vertex.");
    }
 
    _restoreInternalVertex(iterator, coords);
 
    return iterator;
}
#endif

void PatchKernel::_restoreInternalVertex(const VertexIterator &iterator, const std::array<double, 3> &coords)
{
    // Restore the vertex
    //
    // There is no need to set the id of the vertex as assigned, because
    // also the index generator will be restored.
    Vertex &vertex = *iterator;
    vertex.initialize(iterator.getId(), coords, true);
    m_nInternalVertices++;
 
    // Update the bounding box
    addPointToBoundingBox(vertex.getCoords());
}

bool PatchKernel::deleteVertex(long id)
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    // Delete the vertex
#if BITPIT_ENABLE_MPI==1
    const Vertex &vertex = m_vertices[id];
    bool isInternal = vertex.isInterior();
    if (isInternal) {
        _deleteInternalVertex(id);
    } else {
        _deleteGhostVertex(id);
    }
#else
    _deleteInternalVertex(id);
#endif
 
    return true;
}

void PatchKernel::_deleteInternalVertex(long id)
{
    // Update the bounding box
    const Vertex &vertex = m_vertices[id];
    removePointFromBoundingBox(vertex.getCoords());
 
    // Delete vertex
    m_vertices.erase(id, true);
    m_nInternalVertices--;
    if (id == m_lastInternalVertexId) {
        updateLastInternalVertexId();
    }
 
    // Vertex id is no longer used
    if (m_vertexIdGenerator) {
        m_vertexIdGenerator->trash(id);
    }
}

long PatchKernel::countFreeVertices() const
{
    return countBorderVertices();
}

long PatchKernel::countBorderVertices() const
{
    std::unordered_set<long> borderVertices;
    for (const Cell &cell : m_cells) {
        int nCellFaces = cell.getFaceCount();
        for (int i = 0; i < nCellFaces; ++i) {
            if (!cell.isFaceBorder(i)) {
                continue;
            }
 
            int nFaceVertices = cell.getFaceVertexCount(i);
            for (int k = 0; k < nFaceVertices; ++k) {
                long faceVertexId = cell.getFaceVertexId(i, k);
                borderVertices.insert(faceVertexId);
            }
        }
    }
 
    return borderVertices.size();
}

long PatchKernel::countOrphanVertices() const
{
    std::unordered_set<long> usedVertices;
    for (const Cell &cell : m_cells) {
        ConstProxyVector<long> cellVertexIds = cell.getVertexIds();
        int nCellVertices = cellVertexIds.size();
        for (int i = 0; i < nCellVertices; ++i) {
            usedVertices.insert(cellVertexIds[i]);
        }
    }
 
    return (getVertexCount() - usedVertices.size());
}

std::vector<long> PatchKernel::findOrphanVertices()
{
    VertexConstIterator beginItr = m_vertices.cbegin();
    VertexConstIterator endItr   = m_vertices.cend();
 
    // Detect used vertices
    PiercedStorage<bool, long> vertexUsedFlag(1, &m_vertices);
    vertexUsedFlag.fill(false);
    for (const Cell &cell : m_cells) {
        ConstProxyVector<long> cellVertexIds = cell.getVertexIds();
        int nCellVertices = cellVertexIds.size();
        for (int j = 0; j < nCellVertices; ++j) {
            long vertexId = cellVertexIds[j];
            vertexUsedFlag[vertexId] = true;
        }
    }
 
    // Count the orphan vertices
    std::size_t nOrhpanVertices = 0;
    for (VertexConstIterator itr = beginItr; itr != endItr; ++itr) {
        std::size_t vertexRawIndex = itr.getRawIndex();
        if (!vertexUsedFlag.rawAt(vertexRawIndex)) {
            ++nOrhpanVertices;
        }
    }
 
    // Build a list of unused vertices
    std::vector<long> orhpanVertices(nOrhpanVertices);
 
    std::size_t orphanVertexIndex = 0;
    for (VertexConstIterator itr = beginItr; itr != endItr; ++itr) {
        std::size_t vertexRawIndex = itr.getRawIndex();
        if (!vertexUsedFlag.rawAt(vertexRawIndex)) {
            orhpanVertices[orphanVertexIndex] = itr.getId();
            ++orphanVertexIndex;
        }
    }
 
    return orhpanVertices;
}

bool PatchKernel::deleteOrphanVertices()
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    std::vector<long> list = findOrphanVertices();
    deleteVertices(list);
    updateBoundingBox();
 
    return true;
}

std::vector<long> PatchKernel::collapseCoincidentVertices()
{
    std::vector<long> collapsedVertices;
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return collapsedVertices;
    }
 
    long nVertices = getVertexCount();
    if (nVertices == 0) {
        return collapsedVertices;
    }
 
    // Update bounding box
    updateBoundingBox();
 
    //
    // Collapse double vertices
    //
    Vertex::Less vertexLess(10 * std::numeric_limits<double>::epsilon());
    auto rawVertexLess = [this, &vertexLess](const std::size_t &i, const std::size_t &j)
    {
        return vertexLess(this->m_vertices.rawAt(i), this->m_vertices.rawAt(j));
    };
    std::set<std::size_t, decltype(rawVertexLess)> vertexTree(rawVertexLess);
 
    std::unordered_map<long, long> vertexMap;
    for (VertexConstIterator vertexItr = m_vertices.cbegin(); vertexItr != m_vertices.cend(); ++vertexItr) {
        std::size_t vertexRawId = vertexItr.getRawIndex();
        auto vertexTreeItr = vertexTree.find(vertexRawId);
        if (vertexTreeItr == vertexTree.end()) {
            vertexTree.insert(vertexRawId);
        } else {
            long vertexId = vertexItr.getId();
            long vertexCoincidentId = m_vertices.rawFind(*vertexTreeItr).getId();
            vertexMap.insert({vertexId, vertexCoincidentId});
        }
    }
 
    // Collapse vertices
    if (!vertexMap.empty()) {
        // Update cells
        //
        // If the adjacencies are currently up-to-date, they should kept up-to-date.
        AdjacenciesBuildStrategy adjacenciesBuildStrategy = getAdjacenciesBuildStrategy();
 
        bool keepAdjacenciesUpToDate;
        if (adjacenciesBuildStrategy != ADJACENCIES_NONE) {
            keepAdjacenciesUpToDate = (!areAdjacenciesDirty());
        } else {
            keepAdjacenciesUpToDate = false;
        }
 
        for (Cell &cell : m_cells) {
            // Renumber cell vertices
            int nRenumberedVertices = cell.renumberVertices(vertexMap);
 
            // Mark adjacencies are dirty
            //
            // If some vertices have been renumbered, the adjacencies of the cells
            // are now dirty.
            if ((adjacenciesBuildStrategy != ADJACENCIES_NONE) && (nRenumberedVertices > 0)) {
                setCellAlterationFlags(cell.getId(), FLAG_ADJACENCIES_DIRTY);
            }
        }
 
        if (keepAdjacenciesUpToDate) {
            updateAdjacencies();
        }
 
        // Update interface
        for (Interface &interface : m_interfaces) {
            // Renumber interface vertices
            interface.renumberVertices(vertexMap);
        }
 
        // Create the list of collapsed vertices
        collapsedVertices.resize(vertexMap.size());
 
        std::size_t k = 0;
        for (const auto &entry : vertexMap) {
            collapsedVertices[k++] = entry.first;
        }
    }
 
    return collapsedVertices;
}

bool PatchKernel::deleteCoincidentVertices()
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    std::vector<long> verticesToDelete = collapseCoincidentVertices();
    deleteVertices(verticesToDelete);
 
    return true;
}

const std::array<double, 3> & PatchKernel::getVertexCoords(long id) const
{
    return getVertex(id).getCoords();
}

void PatchKernel::getVertexCoords(std::size_t nVertices, const long *ids, std::unique_ptr<std::array<double, 3>[]> *coordinates) const
{
    *coordinates = std::unique_ptr<std::array<double, 3>[]>(new std::array<double, 3>[nVertices]);
    getVertexCoords(nVertices, ids, coordinates->get());
}

void PatchKernel::getVertexCoords(std::size_t nVertices, const long *ids, std::array<double, 3> *coordinates) const
{
    for (std::size_t i = 0; i < nVertices; ++i) {
        coordinates[i] = getVertex(ids[i]).getCoords();
    }
}

bool PatchKernel::empty(bool global) const
{
    bool isEmpty = (getCellCount() == 0);
#if BITPIT_ENABLE_MPI==1
    if (global && isPartitioned()) {
        MPI_Allreduce(MPI_IN_PLACE, &isEmpty, 1, MPI_C_BOOL, MPI_LAND, getCommunicator());
    }
#else
    BITPIT_UNUSED(global);
#endif
 
    return isEmpty;
}

void PatchKernel::updateLastInternalVertexId()
{
    if (m_nInternalVertices == 0) {
        m_lastInternalVertexId = Vertex::NULL_ID;
        return;
    }
 
    VertexIterator lastInternalVertexItr;
#if BITPIT_ENABLE_MPI==1
    if (m_nGhostVertices == 0) {
        lastInternalVertexItr = --m_vertices.end();
        m_lastInternalVertexId = lastInternalVertexItr->getId();
    } else {
        m_lastInternalVertexId = m_vertices.getSizeMarker(m_nInternalVertices - 1, Vertex::NULL_ID);
    }
#else
    lastInternalVertexItr = --m_vertices.end();
    m_lastInternalVertexId = lastInternalVertexItr->getId();
#endif
}

bool PatchKernel::isCellAutoIndexingEnabled() const
{
    return static_cast<bool>(m_cellIdGenerator);
}

void PatchKernel::setCellAutoIndexing(bool enabled)
{
    if (isCellAutoIndexingEnabled() == enabled) {
        return;
    }
 
    if (enabled) {
        createCellIndexGenerator(true);
    } else {
        m_cellIdGenerator.reset();
    }
}

void PatchKernel::dumpCellAutoIndexing(std::ostream &stream) const
{
    m_cellIdGenerator->dump(stream);
}

void PatchKernel::restoreCellAutoIndexing(std::istream &stream)
{
    createCellIndexGenerator(false);
    m_cellIdGenerator->restore(stream);
}

void PatchKernel::createCellIndexGenerator(bool populate)
{
    // Create or reset index generator
    if (!m_cellIdGenerator) {
        m_cellIdGenerator = std::unique_ptr<IndexGenerator<long>>(new IndexGenerator<long>());
    } else {
        m_cellIdGenerator->reset();
    }
 
    // Populate index generator with current ids
    if (populate) {
        CellConstIterator beginItr = m_cells.cbegin();
        CellConstIterator endItr   = m_cells.cend();
        for (CellConstIterator itr = beginItr; itr != endItr; ++itr) {
            m_cellIdGenerator->setAssigned(itr.getId());
        }
    }
}

void PatchKernel::importCellIndexGenerator(const PatchKernel &source)
{
    if (source.m_cellIdGenerator) {
        m_cellIdGenerator = std::unique_ptr<IndexGenerator<long>>(new IndexGenerator<long>(*(source.m_cellIdGenerator)));
    } else {
        m_cellIdGenerator.reset();
    }
}

long PatchKernel::getCellCount() const
{
    return m_cells.size();
}

long PatchKernel::getInternalCellCount() const
{
    return m_nInternalCells;
}

long PatchKernel::getInternalCount() const
{
    return getInternalCellCount();
}

PiercedVector<Cell> & PatchKernel::getCells()
{
    return m_cells;
}

const PiercedVector<Cell> & PatchKernel::getCells() const
{
    return m_cells;
}

Cell & PatchKernel::getCell(long id)
{
    return m_cells[id];
}

const Cell & PatchKernel::getCell(long id) const
{
    return m_cells[id];
}

ElementType PatchKernel::getCellType(long id) const
{
    return m_cells[id].getType();
}

Cell & PatchKernel::getLastInternalCell()
{
    return m_cells[m_lastInternalCellId];
}

const Cell & PatchKernel::getLastInternalCell() const
{
    return m_cells[m_lastInternalCellId];
}

PatchKernel::CellIterator PatchKernel::getCellIterator(long id)
{
    return m_cells.find(id);
}

PatchKernel::CellIterator PatchKernel::cellBegin()
{
    return m_cells.begin();
}

PatchKernel::CellIterator PatchKernel::cellEnd()
{
    return m_cells.end();
}

PatchKernel::CellIterator PatchKernel::internalCellBegin()
{
    return m_cells.begin();
}

PatchKernel::CellIterator PatchKernel::internalBegin()
{
    return internalCellBegin();
}

PatchKernel::CellIterator PatchKernel::internalCellEnd()
{
    if (m_nInternalCells > 0) {
        return ++m_cells.find(m_lastInternalCellId);
    } else {
        return m_cells.end();
    }
}

PatchKernel::CellIterator PatchKernel::internalEnd()
{
    return internalCellEnd();
}

PatchKernel::CellConstIterator PatchKernel::getCellConstIterator(long id) const
{
    return m_cells.find(id);
}

PatchKernel::CellConstIterator PatchKernel::cellConstBegin() const
{
    return m_cells.cbegin();
}

PatchKernel::CellConstIterator PatchKernel::cellConstEnd() const
{
    return m_cells.cend();
}

PatchKernel::CellConstIterator PatchKernel::internalCellConstBegin() const
{
    return m_cells.cbegin();
}

PatchKernel::CellConstIterator PatchKernel::internalConstBegin() const
{
    return internalCellConstBegin();
}

PatchKernel::CellConstIterator PatchKernel::internalCellConstEnd() const
{
    if (m_nInternalCells > 0) {
        return ++m_cells.find(m_lastInternalCellId);
    } else {
        return m_cells.cend();
    }
}

PatchKernel::CellConstIterator PatchKernel::internalConstEnd() const
{
    return internalCellConstEnd();
}

PatchKernel::CellIterator PatchKernel::addCell(const Cell &source, long id)
{
    Cell cell = source;
    cell.setId(id);
 
    return addCell(std::move(cell), id);
}

PatchKernel::CellIterator PatchKernel::addCell(Cell &&source, long id)
{
    // Get the id
    if (id < 0) {
        id = source.getId();
    }
 
    // Add a dummy cell
    //
    // It is not possible to directly add the source into the storage. First a
    // dummy cell is created and then that cell is replaced with the source.
    std::unique_ptr<long[]> dummyConnectStorage = std::unique_ptr<long[]>(nullptr);
 
    CellIterator iterator = _addInternalCell(ElementType::UNDEFINED, std::move(dummyConnectStorage), id);
 
    // Replace the newly created cell with the source
    //
    // Before replacing the newly created cell we need to set the id of the
    // source to the id that has been assigned to the newly created cell.
    source.setId(iterator->getId());
 
    Cell &cell = (*iterator);
    cell = std::move(source);
 
    return iterator;
}

PatchKernel::CellIterator PatchKernel::addCell(ElementType type, long id)
{
    std::unique_ptr<long[]> connectStorage;
    if (ReferenceElementInfo::hasInfo(type)) {
        int connectSize = ReferenceElementInfo::getInfo(type).nVertices;
        connectStorage = std::unique_ptr<long[]>(new long[connectSize]);
    } else {
        connectStorage = std::unique_ptr<long[]>(nullptr);
    }
 
    return addCell(type, std::move(connectStorage), id);
}

PatchKernel::CellIterator PatchKernel::addCell(ElementType type, const std::vector<long> &connectivity,
                                               long id)
{
    int connectSize = connectivity.size();
    std::unique_ptr<long[]> connectStorage = std::unique_ptr<long[]>(new long[connectSize]);
    std::copy(connectivity.data(), connectivity.data() + connectSize, connectStorage.get());
 
    return addCell(type, std::move(connectStorage), id);
}

PatchKernel::CellIterator PatchKernel::addCell(ElementType type, std::unique_ptr<long[]> &&connectStorage,
                                               long id)
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return cellEnd();
    }
 
    if (Cell::getDimension(type) > getDimension()) {
        return cellEnd();
    }
 
    CellIterator iterator = _addInternalCell(type, std::move(connectStorage), id);
 
    return iterator;
}

PatchKernel::CellIterator PatchKernel::_addInternalCell(ElementType type, std::unique_ptr<long[]> &&connectStorage,
                                                    long id)
{
    // Get the id of the cell
    if (m_cellIdGenerator) {
        if (id < 0) {
            id = m_cellIdGenerator->generate();
        } else {
            m_cellIdGenerator->setAssigned(id);
        }
    } else if (id < 0) {
        throw std::runtime_error("No valid id has been provided for the cell.");
    }
 
    // Get the id of the cell before which the new cell should be inserted
#if BITPIT_ENABLE_MPI==1
    //
    // If there are ghosts cells, the internal cell should be inserted
    // before the first ghost cell.
#endif
    long referenceId;
#if BITPIT_ENABLE_MPI==1
    referenceId = m_firstGhostCellId;
#else
    referenceId = Cell::NULL_ID;
#endif
 
    // Create the cell
    bool storeInterfaces  = (getInterfacesBuildStrategy() != INTERFACES_NONE);
    bool storeAdjacencies = storeInterfaces || (getAdjacenciesBuildStrategy() != ADJACENCIES_NONE);
 
    CellIterator iterator;
    if (referenceId == Cell::NULL_ID) {
        iterator = m_cells.emreclaim(id, id, type, std::move(connectStorage), true, storeInterfaces, storeAdjacencies);
    } else {
        iterator = m_cells.emreclaimBefore(referenceId, id, id, type, std::move(connectStorage), true, storeInterfaces, storeAdjacencies);
    }
    m_nInternalCells++;
 
    // Update the id of the last internal cell
    if (m_lastInternalCellId < 0) {
        m_lastInternalCellId = id;
    } else if (m_cells.rawIndex(m_lastInternalCellId) < m_cells.rawIndex(id)) {
        m_lastInternalCellId = id;
    }
 
    // Set the alteration flags of the cell
    setAddedCellAlterationFlags(id);
 
#if BITPIT_ENABLE_MPI==1
    // Set partitioning information as dirty
    setPartitioningInfoDirty(true);
#endif
 
    return iterator;
}

void PatchKernel::setAddedCellAlterationFlags(long id)
{
    AlterationFlags flags = FLAG_NONE;
    if (getAdjacenciesBuildStrategy() != ADJACENCIES_NONE) {
        flags |= FLAG_ADJACENCIES_DIRTY;
    }
    if (getInterfacesBuildStrategy() != INTERFACES_NONE) {
        flags |= FLAG_INTERFACES_DIRTY;
    }
 
    setCellAlterationFlags(id, flags);
}
 
#if BITPIT_ENABLE_MPI==0
PatchKernel::CellIterator PatchKernel::restoreCell(ElementType type, std::unique_ptr<long[]> &&connectStorage,
                                                   long id)
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return cellEnd();
    }
 
    if (Cell::getDimension(type) > getDimension()) {
        return cellEnd();
    }
 
    CellIterator iterator = m_cells.find(id);
    if (iterator == m_cells.end()) {
        throw std::runtime_error("Unable to restore the specified cell: the kernel doesn't contain an entry for that cell.");
    }
 
    _restoreInternalCell(iterator, type, std::move(connectStorage));
 
    return iterator;
}
#endif

void PatchKernel::_restoreInternalCell(const CellIterator &iterator, ElementType type,
                                   std::unique_ptr<long[]> &&connectStorage)
{
    // Restore the cell
    //
    // There is no need to set the id of the cell as assigned, because
    // also the index generator will be restored.
    bool storeInterfaces  = (getInterfacesBuildStrategy() != INTERFACES_NONE);
    bool storeAdjacencies = storeInterfaces || (getAdjacenciesBuildStrategy() != ADJACENCIES_NONE);
 
    long cellId = iterator.getId();
    Cell &cell = *iterator;
    cell.initialize(cellId, type, std::move(connectStorage), true, storeInterfaces, storeAdjacencies);
    m_nInternalCells++;
 
    // Set the alteration flags of the cell
    setRestoredCellAlterationFlags(cellId);
}

void PatchKernel::setRestoredCellAlterationFlags(long id)
{
    AlterationFlags flags = FLAG_NONE;
    if (getAdjacenciesBuildStrategy() != ADJACENCIES_NONE) {
        flags |= FLAG_ADJACENCIES_DIRTY;
    }
    if (getInterfacesBuildStrategy() != INTERFACES_NONE) {
        flags |= FLAG_INTERFACES_DIRTY;
    }
 
    setCellAlterationFlags(id, flags);
}

bool PatchKernel::deleteCell(long id)
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    // Delete cell
#if BITPIT_ENABLE_MPI==1
    const Cell &cell = m_cells[id];
    bool isInternalCell = cell.isInterior();
    if (isInternalCell) {
        _deleteInternalCell(id);
    } else {
        _deleteGhostCell(id);
    }
#else
    _deleteInternalCell(id);
#endif
 
    return true;
}

void PatchKernel::_deleteInternalCell(long id)
{
    // Set the alteration flags of the cell
    setDeletedCellAlterationFlags(id);
 
    // Delete cell
    m_cells.erase(id, true);
    m_nInternalCells--;
    if (id == m_lastInternalCellId) {
        updateLastInternalCellId();
    }
 
    // Cell id is no longer used
    if (m_cellIdGenerator) {
        m_cellIdGenerator->trash(id);
    }
}

void PatchKernel::setDeletedCellAlterationFlags(long id)
{
    const Cell &cell = getCell(id);
 
    // Set the alteration flags of the cell
    resetCellAlterationFlags(id, FLAG_DELETED);
 
    // Set the alteration flags of the adjacent cells
    const int nCellAdjacencies = cell.getAdjacencyCount();
    const long *cellAdjacencies = cell.getAdjacencies();
    for (int k = 0; k < nCellAdjacencies; ++k) {
        long adjacencyId = cellAdjacencies[k];
        if (!testCellAlterationFlags(adjacencyId, FLAG_DELETED)) {
            AlterationFlags flags = FLAG_DANGLING;
            if (getAdjacenciesBuildStrategy() != ADJACENCIES_NONE) {
                flags |= FLAG_ADJACENCIES_DIRTY;
            }
            if (getInterfacesBuildStrategy() != INTERFACES_NONE) {
                flags |= FLAG_INTERFACES_DIRTY;
            }
 
            setCellAlterationFlags(adjacencyId, flags);
        }
    }
 
    // Set the alteration flags of the interfaces
    const int nCellInterfaces = cell.getInterfaceCount();
    const long *cellInterfaces = cell.getInterfaces();
    for (int k = 0; k < nCellInterfaces; ++k) {
        long interfaceId = cellInterfaces[k];
        if (!testInterfaceAlterationFlags(interfaceId, FLAG_DELETED)) {
            setInterfaceAlterationFlags(interfaceId, FLAG_DANGLING);
        }
    }
}

long PatchKernel::countFreeCells() const
{
    return countBorderCells();
}

long PatchKernel::countBorderCells() const
{
    long nBorderCells = 0;
    for (const Cell &cell : m_cells) {
        int nCellFaces = cell.getFaceCount();
        for (int i = 0; i < nCellFaces; ++i) {
            if (cell.isFaceBorder(i)) {
                ++nBorderCells;
                break;
            }
        }
    }
 
    return nBorderCells;
}

long PatchKernel::countOrphanCells() const
{
    return findOrphanCells().size();
}

std::vector<long> PatchKernel::findOrphanCells() const
{
    // Compute vertex valence
    std::unordered_map<long, short> vertexValence;
    for (const Cell &cell : m_cells) {
        ConstProxyVector<long> cellVertexIds = cell.getVertexIds();
        int nCellVertices = cellVertexIds.size();
        for (int j = 0; j < nCellVertices; j++) {
            long vertexId = cellVertexIds[j];
            vertexValence[vertexId] += 1;
        }
    }
 
    // Loop over cells
    std::vector<long> orphanCells;
    for (const Cell &cell : m_cells) {
        long isIsolated = true;
        ConstProxyVector<long> cellVertexIds = cell.getVertexIds();
        int nCellVertices = cellVertexIds.size();
        for (int j = 0; j < nCellVertices; j++) {
            long vertexId = cellVertexIds[j];
            if (vertexValence[vertexId] > 1) {
                isIsolated = false;
                break;
            }
        }
 
        if (isIsolated) {
            orphanCells.push_back(cell.getId());
        }
    }
 
    return orphanCells;
}

long PatchKernel::countDuplicateCells() const
{
    return findDuplicateCells().size();
}

std::vector<long> PatchKernel::findDuplicateCells() const
{
    // Define the hasher and the predicate to be used for lists of vertex ids
    //
    // These operators should be able to identify that two lists of vertex ids
    // are the same regardless of the order in which the ids are listed.
    auto hasher = [](const ConstProxyVector<long> &ids)
    {
        std::size_t hash = std::hash<long>{}(0);
        for (long id : ids) {
            hash = hash + std::hash<long>{}(id);
        }
 
        return hash;
    };
 
    std::unordered_map<long, std::size_t > counters;
    auto predicate = [&counters](const ConstProxyVector<long> &ids_A, const ConstProxyVector<long> &ids_B)
    {
        counters.clear();
        for (long id : ids_A) {
            counters[id]++;
        }
        for (long id : ids_B) {
            counters[id]--;
        }
 
        for (const auto &entry : counters) {
            if (entry.second != 0) {
                return false;
            }
        }
 
        return true;
    };
 
    // Detect if there are cells that share the same list of vertices
    //
    // For each cell, the list of vertex ids is added into a set. If a collision
    // is detected, we have found a duplicate cell.
    std::vector<long> duplicateCells;
    std::unordered_set<ConstProxyVector<long>, decltype(hasher), decltype(predicate)> vertexStash(getCellCount(), hasher, predicate);
    for (const Cell &cell : m_cells) {
        ConstProxyVector<long> cellVertexIds = cell.getVertexIds();
        auto vertexStashItr = vertexStash.find(cellVertexIds);
        if (vertexStashItr == vertexStash.end()) {
            vertexStash.insert(std::move(cellVertexIds));
        } else {
            duplicateCells.push_back(cell.getId());
        }
    }
 
    return duplicateCells;
}

void PatchKernel::updateLastInternalCellId()
{
    if (m_nInternalCells == 0) {
        m_lastInternalCellId = Cell::NULL_ID;
        return;
    }
 
    CellIterator lastInternalCellItr;
#if BITPIT_ENABLE_MPI==1
    if (m_nGhostCells == 0) {
        lastInternalCellItr = --m_cells.end();
        m_lastInternalCellId = lastInternalCellItr->getId();
    } else {
        m_lastInternalCellId = m_cells.getSizeMarker(m_nInternalCells - 1, Cell::NULL_ID);
    }
#else
    lastInternalCellItr = --m_cells.end();
    m_lastInternalCellId = lastInternalCellItr->getId();
#endif
}

long PatchKernel::_getCellNativeIndex(long id) const
{
    return id;
}

std::vector<long> PatchKernel::findCellNeighs(long id) const
{
    std::vector<long> neighs;
    findCellNeighs(id, &neighs);
 
    return neighs;
}

void PatchKernel::findCellNeighs(long id, std::vector<long> *neighs) const
{
    _findCellNeighs(id, nullptr, neighs);
}

std::vector<long> PatchKernel::findCellNeighs(long id, int codimension, bool complete) const
{
    std::vector<long> neighs;
    findCellNeighs(id, codimension, complete, &neighs);
 
    return neighs;
}

void PatchKernel::findCellNeighs(long id, int codimension, bool complete, std::vector<long> *neighs) const
{
    assert(codimension >= 1 && codimension <= getDimension());
 
    if (codimension == 1) {
        findCellFaceNeighs(id, neighs);
    } else if (codimension == getDimension()) {
        findCellVertexNeighs(id, complete, neighs);
    } else if (codimension == 2) {
        findCellEdgeNeighs(id, complete, neighs);
    }
}

std::vector<long> PatchKernel::findCellFaceNeighs(long id) const
{
    std::vector<long> neighs;
    findCellFaceNeighs(id, &neighs);
 
    return neighs;
}

void PatchKernel::_findCellNeighs(long id, const std::vector<long> *blackList, std::vector<long> *neighs) const
{
    // Some patches can work (at least partially) without initializing the
    // cell list. To handle those patches, if there are no cells the vertex
    // count is evaluated using the ReferenceElementInfo associated to the cell.
    int nCellVertices;
    if (m_cells.size() == 0) {
        ElementType cellType = getCellType(id);
        const ReferenceElementInfo &cellTypeInfo = ReferenceElementInfo::getInfo(cellType);
        nCellVertices = cellTypeInfo.nVertices;
    } else {
        const Cell &cell = getCell(id);
        nCellVertices = cell.getVertexCount();
    }
 
    for (int i = 0; i < nCellVertices; ++i) {
        _findCellVertexNeighs(id, i, blackList, neighs);
    }
}

void PatchKernel::findCellFaceNeighs(long id, std::vector<long> *neighs) const
{
    // Some patches can work (at least partially) without initializing the
    // cell list. To handle those patches, if there are no cells the face
    // count is evaluated using the ReferenceElementInfo associated to the cell.
    int nCellFaces;
    if (m_cells.size() == 0) {
        ElementType cellType = getCellType(id);
        const ReferenceElementInfo &cellTypeInfo = ReferenceElementInfo::getInfo(cellType);
        nCellFaces = cellTypeInfo.nFaces;
    } else {
        const Cell &cell = getCell(id);
        nCellFaces = cell.getFaceCount();
    }
 
    // Get the neighbours
    for (int i = 0; i < nCellFaces; ++i) {
        _findCellFaceNeighs(id, i, nullptr, neighs);
    }
}

std::vector<long> PatchKernel::findCellFaceNeighs(long id, int face) const
{
    std::vector<long> neighs;
    _findCellFaceNeighs(id, face, nullptr, &neighs);
 
    return neighs;
}

void PatchKernel::findCellFaceNeighs(long id, int face, std::vector<long> *neighs) const
{
    _findCellFaceNeighs(id, face, nullptr, neighs);
}

void PatchKernel::_findCellFaceNeighs(long id, int face, const std::vector<long> *blackList, std::vector<long> *neighs) const
{
    const Cell &cell = getCell(id);
 
    int nFaceAdjacencies = cell.getAdjacencyCount(face);
    const long *faceAdjacencies = cell.getAdjacencies(face);
    for (int k = 0; k < nFaceAdjacencies; ++k) {
        long neighId = faceAdjacencies[k];
        if (!blackList || utils::findInOrderedVector<long>(neighId, *blackList) == blackList->end()) {
            utils::addToOrderedVector<long>(neighId, *neighs);
        }
    }
}

std::vector<long> PatchKernel::findCellEdgeNeighs(long id, bool complete) const
{
    std::vector<long> neighs;
    findCellEdgeNeighs(id, complete, &neighs);
 
    return neighs;
}

void PatchKernel::findCellEdgeNeighs(long id, bool complete, std::vector<long> *neighs) const
{
    assert(isThreeDimensional());
    if (!isThreeDimensional()) {
        return;
    }
 
    // Some patches can work (at least partially) without initializing the
    // cell list. To handle those patches, if there are no cells the edge
    // count is evaluated using the ReferenceElementInfo associated to the cell.
    int nCellEdges;
    if (m_cells.size() == 0) {
        ElementType cellType = getCellType(id);
        const ReferenceElementInfo &cellTypeInfo = ReferenceElementInfo::getInfo(cellType);
        nCellEdges = cellTypeInfo.nEdges;
    } else {
        const Cell &cell = getCell(id);
        nCellEdges = cell.getEdgeCount();
    }
 
    // Get the neighbours
    std::unique_ptr<std::vector<long>> blackList;
    if (!complete) {
        blackList = std::unique_ptr<std::vector<long>>(new std::vector<long>());
        findCellFaceNeighs(id, blackList.get());
    }
 
    for (int i = 0; i < nCellEdges; ++i) {
        _findCellEdgeNeighs(id, i, blackList.get(), neighs);
    }
}

std::vector<long> PatchKernel::findCellEdgeNeighs(long id, int edge) const
{
    std::vector<long> neighs;
    findCellEdgeNeighs(id, edge, &neighs);
 
    return neighs;
}

void PatchKernel::findCellEdgeNeighs(long id, int edge, std::vector<long> *neighs) const
{
    _findCellEdgeNeighs(id, edge, nullptr, neighs);
}

void PatchKernel::_findCellEdgeNeighs(long id, int edge, const std::vector<long> *blackList, std::vector<long> *neighs) const
{
    assert(isThreeDimensional());
    if (!isThreeDimensional()) {
        return;
    }
 
    const Cell &cell = getCell(id);
    ConstProxyVector<int> edgeVertices = cell.getEdgeLocalConnect(edge);
    std::size_t nEdgeVertices = edgeVertices.size();
    if (nEdgeVertices < 2) {
        return;
    }
 
    // Find candidates neighbours
    //
    // These are the negihbours of the first vertex of the edge.
    const int GUESS_NEIGHS_COUNT = 3;
 
    std::vector<long> candidateIds;
    candidateIds.reserve(GUESS_NEIGHS_COUNT);
    _findCellVertexNeighs(id, edgeVertices[0], blackList, &candidateIds);
 
    // Discard candidates that doesn't contain the last vertex of the edge
    long lastEdgeVertexId = cell.getEdgeVertexId(edge, nEdgeVertices - 1);
 
    ConstProxyVector<long> candidateVertexIds;
    for (long candidateId : candidateIds) {
        const Cell &candidateNeigh = m_cells.at(candidateId);
        candidateVertexIds = candidateNeigh.getVertexIds();
        std::size_t nCandidateVertices = candidateVertexIds.size();
 
        bool isEdgeNeighbour = false;
        for (std::size_t k = 0; k < nCandidateVertices; ++k) {
            if (candidateVertexIds[k] == lastEdgeVertexId) {
                isEdgeNeighbour = true;
                break;
            }
        }
 
        if (isEdgeNeighbour) {
            utils::addToOrderedVector<long>(candidateId, *neighs);
        }
    }
}

std::vector<long> PatchKernel::findCellVertexNeighs(long id, bool complete) const
{
    std::vector<long> neighs;
    findCellVertexNeighs(id, complete, &neighs);
 
    return neighs;
}

void PatchKernel::findCellVertexNeighs(long id, bool complete, std::vector<long> *neighs) const
{
    // Some patches can work (at least partially) without initializing the
    // cell list. To handle those patches, if there are no cells the vertex
    // count is evaluated using the ReferenceElementInfo associated to the cell.
    int nCellVertices;
    if (m_cells.size() == 0) {
        ElementType cellType = getCellType(id);
        const ReferenceElementInfo &cellTypeInfo = ReferenceElementInfo::getInfo(cellType);
        nCellVertices = cellTypeInfo.nVertices;
    } else {
        const Cell &cell = getCell(id);
        nCellVertices = cell.getVertexCount();
    }
 
    // Get the neighbours
    std::unique_ptr<std::vector<long>> blackList;
    if (!complete) {
        blackList = std::unique_ptr<std::vector<long>>(new std::vector<long>());
        if (isThreeDimensional()) {
            findCellEdgeNeighs(id, true, blackList.get());
        } else {
            findCellFaceNeighs(id, true, blackList.get());
        }
    }
 
    for (int i = 0; i < nCellVertices; ++i) {
        _findCellVertexNeighs(id, i, blackList.get(), neighs);
    }
}

std::vector<long> PatchKernel::findCellVertexNeighs(long id, int vertex) const
{
    std::vector<long> neighs;
    findCellVertexNeighs(id, vertex, &neighs);
 
    return neighs;
}

void PatchKernel::findCellVertexNeighs(long id, int vertex, std::vector<long> *neighs) const
{
    _findCellVertexNeighs(id, vertex, nullptr, neighs);
}

void PatchKernel::_findCellVertexNeighs(long id, int vertex, const std::vector<long> *blackList, std::vector<long> *neighs) const
{
    const Cell &cell = getCell(id);
    long vertexId = cell.getVertexId(vertex);
 
    // Since we are processing a small number of cells it's more efficient to
    // store the list of already processed cells in a vector instead of using
    // a set or an unordered_set. To speed-up the lookup, the vector is kept
    // sorted.
    const int GUESS_NEIGHS_COUNT = 8;
 
    std::vector<long> alreadyProcessed;
    alreadyProcessed.reserve(GUESS_NEIGHS_COUNT);
 
    std::vector<long> scanQueue;
    scanQueue.reserve(GUESS_NEIGHS_COUNT);
    scanQueue.push_back(id);
 
    while (!scanQueue.empty()) {
        // Pop a cell to process
        long scanCellId = scanQueue.back();
        const Cell &scanCell = getCell(scanCellId);
        scanQueue.pop_back();
        utils::addToOrderedVector<long>(scanCell.getId(), alreadyProcessed);
 
        // Get cell information
        const ReferenceElementInfo *scanCellInfo = nullptr;
        const long *scanCellConnectivity = nullptr;
        if (scanCell.hasInfo()) {
            scanCellInfo = &(scanCell.getInfo());
            scanCellConnectivity = scanCell.getConnect();
        }
 
        // Use face adjacencies to find vertex negihbours
        int nFaces;
        if (scanCellInfo) {
            nFaces = scanCellInfo->nFaces;
        } else {
            nFaces = scanCell.getFaceCount();
        }
 
        for (int i = 0; i < nFaces; ++i) {
            // Discard faces with no neighbours
            if (scanCell.isFaceBorder(i)) {
                continue;
            }
 
            // Discard faces that don't own the vertex
            //
            // We use differents algorithm depending on whether the cell has a
            // reference element or not. If the cell has a reference element,
            // accessing the local connectivity of the face is cheap. If the
            // cell has no reference element, it is better to avoid using the
            // local connectivity of the face.
            //
            // Moreover, to optimize the search for cells that has a reference
            // element, we work directly on the connectivity inseatd of on the
            // list of vertex ids. That's because for reference elements the
            // connectivity is just a list of vertices.
            bool faceOwnsVertex = false;
            if (scanCellInfo) {
                const ReferenceElementInfo &faceInfo = ReferenceElementInfo::getInfo(scanCellInfo->faceTypeStorage[i]);
                const int *faceLocalConnectivity = scanCellInfo->faceConnectStorage[i].data();
                const int nFaceVertices = faceInfo.nVertices;
                for (int k = 0; k < nFaceVertices; ++k) {
                    long faceVertexId = scanCellConnectivity[faceLocalConnectivity[k]];
                    if (faceVertexId == vertexId) {
                        faceOwnsVertex = true;
                        break;
                    }
                }
            } else {
                ConstProxyVector<long> faceVertexIds = scanCell.getFaceVertexIds(i);
                int nFaceVertices = faceVertexIds.size();
                for (int k = 0; k < nFaceVertices; ++k) {
                    long faceVertexId = faceVertexIds[k];
                    if (faceVertexId == vertexId) {
                        faceOwnsVertex = true;
                        break;
                    }
                }
            }
 
            if (!faceOwnsVertex) {
                continue;
            }
 
            // Loop through the adjacencies
            //
            // Non-manifold patches may have faces with multiple adjacencies.
            int nFaceAdjacencies = scanCell.getAdjacencyCount(i);
            const long *faceAdjacencies = scanCell.getAdjacencies(i);
            for (int k = 0; k < nFaceAdjacencies; ++k) {
                long faceNeighId = faceAdjacencies[k];
 
                // Discard neighbours that have already been processed
                if (utils::findInOrderedVector<long>(faceNeighId, alreadyProcessed) != alreadyProcessed.end()) {
                    continue;
                }
 
                // Update list of vertex neighbours
                if (!blackList || utils::findInOrderedVector<long>(faceNeighId, *blackList) == blackList->end()) {
                    utils::addToOrderedVector<long>(faceNeighId, *neighs);
                }
 
                // Update scan list
                scanQueue.push_back(faceNeighId);
            }
        }
    }
}

std::vector<long> PatchKernel::findCellVertexOneRing(long id, int vertex) const
{
    std::vector<long> ring;
    findCellVertexOneRing(id, vertex, &ring);
 
    return ring;
}

void PatchKernel::findCellVertexOneRing(long id, int vertex, std::vector<long> *ring) const
{
    findCellVertexNeighs(id, vertex, ring);
    utils::addToOrderedVector<long>(id, *ring);
}

bool PatchKernel::findFaceNeighCell(long cellId, long neighId, int *cellFace, int *cellAdjacencyId) const
{
    const Cell &cell = getCell(cellId);
 
    int nFaces = cell.getFaceCount();
    for (int i = 0; i < nFaces; ++i) {
        int nFaceAdjacencies = cell.getAdjacencyCount(i);
        const long *faceAdjacencies = cell.getAdjacencies(i);
        for (int k = 0; k < nFaceAdjacencies; ++k) {
            long faceNeighId = faceAdjacencies[k];
            if (faceNeighId < 0) {
                continue;
            } else if (faceNeighId == neighId) {
                *cellFace        = i;
                *cellAdjacencyId = k;
                return true;
            }
        }
    }
 
    *cellFace        = -1;
    *cellAdjacencyId = -1;
 
    return false;
}

std::set<int> PatchKernel::getInternalCellPIDs()
{
    std::set<int> list;
    CellConstIterator endItr = internalCellConstEnd();
    for (CellConstIterator itr = internalCellConstBegin(); itr != endItr; ++itr) {
        list.insert(itr->getPID());
    }
 
    return list;
}

std::vector<long> PatchKernel::getInternalCellsByPID(int pid)
{
    std::vector<long> cells;
    CellConstIterator endItr = internalCellConstEnd();
    for (CellConstIterator itr = internalCellConstBegin(); itr != endItr; ++itr) {
        if (itr->getPID() == pid){
            cells.push_back(itr.getId());
        }
    }
 
    return cells;
}

std::vector<long> PatchKernel::findVertexOneRing(long vertexId) const
{
    std::vector<long> ring;
    findVertexOneRing(vertexId, &ring);
 
    return ring;
}

void PatchKernel::findVertexOneRing(long vertexId, std::vector<long> *ring) const
{
    // Find local id of the vertex
    //
    // Coincident vertices with different ids are not supported, this case is
    // special because a cell may contain a point with the vertex coordinates,
    // but may not contain the requested vertex (because it contains one of
    // the other coincident vertices).
    const std::array<double, 3> &coords = getVertexCoords(vertexId);
    long cellId = locatePoint(coords);
    if (cellId == bitpit::Element::NULL_ID) {
        return;
    }
 
    int vertexLocalId = getCell(cellId).findVertex(vertexId);
    assert(vertexLocalId >= 0);
 
    // Find vertex one-ring
    findCellVertexOneRing(cellId, vertexLocalId, ring);
}

bool PatchKernel::isInterfaceAutoIndexingEnabled() const
{
    return static_cast<bool>(m_interfaceIdGenerator);
}

void PatchKernel::setInterfaceAutoIndexing(bool enabled)
{
    if (isInterfaceAutoIndexingEnabled() == enabled) {
        return;
    }
 
    if (enabled) {
        createInterfaceIndexGenerator(true);
    } else {
        if (getInterfacesBuildStrategy() == INTERFACES_AUTOMATIC) {
            throw std::runtime_error("Auto-indexing cannot be disabled if interfaces build strategy is set to automatic.");
        }
 
        m_interfaceIdGenerator.reset();
    }
}

void PatchKernel::dumpInterfaceAutoIndexing(std::ostream &stream) const
{
    m_interfaceIdGenerator->dump(stream);
}

void PatchKernel::restoreInterfaceAutoIndexing(std::istream &stream)
{
    createInterfaceIndexGenerator(false);
    m_interfaceIdGenerator->restore(stream);
}

void PatchKernel::createInterfaceIndexGenerator(bool populate)
{
    // Create or reset index generator
    if (!m_interfaceIdGenerator) {
        m_interfaceIdGenerator = std::unique_ptr<IndexGenerator<long>>(new IndexGenerator<long>());
    } else {
        m_interfaceIdGenerator->reset();
    }
 
    // Populate index generator with current ids
    if (populate) {
        InterfaceConstIterator beginItr = m_interfaces.cbegin();
        InterfaceConstIterator endItr   = m_interfaces.cend();
        for (InterfaceConstIterator itr = beginItr; itr != endItr; ++itr) {
            m_interfaceIdGenerator->setAssigned(itr.getId());
        }
    }
}

void PatchKernel::importInterfaceIndexGenerator(const PatchKernel &source)
{
    if (source.m_interfaceIdGenerator) {
        m_interfaceIdGenerator = std::unique_ptr<IndexGenerator<long>>(new IndexGenerator<long>(*(source.m_interfaceIdGenerator)));
    } else {
        m_interfaceIdGenerator.reset();
    }
}

long PatchKernel::getInterfaceCount() const
{
    return m_interfaces.size();
}

PiercedVector<Interface> & PatchKernel::getInterfaces()
{
    return m_interfaces;
}

const PiercedVector<Interface> & PatchKernel::getInterfaces() const
{
    return m_interfaces;
}

Interface & PatchKernel::getInterface(long id)
{
    return m_interfaces[id];
}

const Interface & PatchKernel::getInterface(long id) const
{
    return m_interfaces[id];
}

ElementType PatchKernel::getInterfaceType(long id) const
{
    return m_interfaces[id].getType();
}

PatchKernel::InterfaceIterator PatchKernel::getInterfaceIterator(long id)
{
    return m_interfaces.find(id);
}

PatchKernel::InterfaceIterator PatchKernel::interfaceBegin()
{
    return m_interfaces.begin();
}

PatchKernel::InterfaceIterator PatchKernel::interfaceEnd()
{
    return m_interfaces.end();
}

PatchKernel::InterfaceConstIterator PatchKernel::getInterfaceConstIterator(long id) const
{
    return m_interfaces.find(id);
}

PatchKernel::InterfaceConstIterator PatchKernel::interfaceConstBegin() const
{
    return m_interfaces.cbegin();
}

PatchKernel::InterfaceConstIterator PatchKernel::interfaceConstEnd() const
{
    return m_interfaces.cend();
}

PatchKernel::InterfaceIterator PatchKernel::addInterface(const Interface &source, long id)
{
    Interface interface = source;
    interface.setId(id);
 
    return addInterface(std::move(interface), id);
}

PatchKernel::InterfaceIterator PatchKernel::addInterface(Interface &&source, long id)
{
    // Get the id
    if (id < 0) {
        id = source.getId();
    }
 
    // Add a dummy interface
    //
    // It is not possible to directly add the source into the storage. First
    // a dummy interface is created and then that interface is replaced with
    // the source.
    std::unique_ptr<long[]> dummyConnectStorage = std::unique_ptr<long[]>(nullptr);
 
    InterfaceIterator iterator = _addInterface(ElementType::UNDEFINED, std::move(dummyConnectStorage), id);
 
    // Replace the newly created interface with the source
    //
    // Before replacing the newly created interface we need to set the id
    // of the source to the id that has been assigned to the newly created
    // interface.
    source.setId(iterator->getId());
 
    Interface &interface = (*iterator);
    interface = std::move(source);
 
    return iterator;
}

PatchKernel::InterfaceIterator PatchKernel::addInterface(ElementType type, long id)
{
    int connectSize = ReferenceElementInfo::getInfo(type).nVertices;
    std::unique_ptr<long[]> connectStorage = std::unique_ptr<long[]>(new long[connectSize]);
 
    return addInterface(type, std::move(connectStorage), id);
}

PatchKernel::InterfaceIterator PatchKernel::addInterface(ElementType type,
                                                         const std::vector<long> &connectivity,
                                                         long id)
{
    int connectSize = connectivity.size();
    std::unique_ptr<long[]> connectStorage = std::unique_ptr<long[]>(new long[connectSize]);
    std::copy(connectivity.data(), connectivity.data() + connectSize, connectStorage.get());
 
    return addInterface(type, std::move(connectStorage), id);
}

PatchKernel::InterfaceIterator PatchKernel::addInterface(ElementType type,
                                                         std::unique_ptr<long[]> &&connectStorage,
                                                         long id)
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return interfaceEnd();
    }
 
    if (Interface::getDimension(type) > (getDimension() - 1)) {
        return interfaceEnd();
    }
 
    InterfaceIterator iterator = _addInterface(type, std::move(connectStorage), id);
 
    return iterator;
}

PatchKernel::InterfaceIterator PatchKernel::_addInterface(ElementType type,
                                                          std::unique_ptr<long[]> &&connectStorage,
                                                          long id)
{
    // Get the id
    if (m_interfaceIdGenerator) {
        if (id < 0) {
            id = m_interfaceIdGenerator->generate();
        } else {
            m_interfaceIdGenerator->setAssigned(id);
        }
    } else if (id < 0) {
        throw std::runtime_error("No valid id has been provided for the interface.");
    }
 
    // Create the interface
    PiercedVector<Interface>::iterator iterator = m_interfaces.emreclaim(id, id, type, std::move(connectStorage));
 
    // Set the alteration flags
    setAddedInterfaceAlterationFlags(id);
 
    return iterator;
}

void PatchKernel::setAddedInterfaceAlterationFlags(long id)
{
    BITPIT_UNUSED(id);
}

PatchKernel::InterfaceIterator PatchKernel::restoreInterface(ElementType type,
                                                             std::unique_ptr<long[]> &&connectStorage,
                                                             long id)
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return interfaceEnd();
    }
 
    if (Interface::getDimension(type) > (getDimension() - 1)) {
        return interfaceEnd();
    }
 
    InterfaceIterator iterator = m_interfaces.find(id);
    if (iterator == m_interfaces.end()) {
        throw std::runtime_error("Unable to restore the specified interface: the kernel doesn't contain an entry for that interface.");
    }
 
    _restoreInterface(iterator, type, std::move(connectStorage));
 
    return iterator;
}

void PatchKernel::_restoreInterface(const InterfaceIterator &iterator, ElementType type,
                                    std::unique_ptr<long[]> &&connectStorage)
{
    // Restore the interface
    //
    // There is no need to set the id of the interfaces as assigned, because
    // also the index generator will be restored.
    long interfaceId = iterator.getId();
    Interface &interface = *iterator;
    interface.initialize(interfaceId, type, std::move(connectStorage));
 
    // Set the alteration flags
    setRestoredInterfaceAlterationFlags(interfaceId);
}

void PatchKernel::setRestoredInterfaceAlterationFlags(long id)
{
    BITPIT_UNUSED(id);
}

bool PatchKernel::deleteInterface(long id)
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    _deleteInterface(id);
 
    return true;
}

void PatchKernel::_deleteInterface(long id)
{
    // Set the alteration flags
    setDeletedInterfaceAlterationFlags(id);
 
    // Delete interface
    m_interfaces.erase(id, true);
 
    // Interface id is no longer used
    if (m_interfaceIdGenerator) {
        m_interfaceIdGenerator->trash(id);
    }
}

void PatchKernel::setDeletedInterfaceAlterationFlags(long id)
{
    resetInterfaceAlterationFlags(id, FLAG_DELETED);
}

long PatchKernel::countFreeInterfaces() const
{
    return countBorderInterfaces();
}

long PatchKernel::countBorderInterfaces() const
{
    long nBorderInterfaces = 0;
    for (const Interface &interface : m_interfaces) {
        if (interface.getNeigh() < 0) {
            ++nBorderInterfaces;
        }
        }
 
    return nBorderInterfaces;
}

long PatchKernel::countOrphanInterfaces() const
{
    long nOrphanInterfaces = 0;
    InterfaceConstIterator endItr = interfaceConstEnd();
    for (InterfaceConstIterator itr = interfaceConstBegin(); itr != endItr; ++itr) {
        const long interfaceId = itr.getId();
        if (isInterfaceOrphan(interfaceId)) {
            ++nOrphanInterfaces;
        }
    }
 
    return nOrphanInterfaces;
}

std::vector<long> PatchKernel::findOrphanInterfaces() const
{
    std::vector<long> orphanInterfaces;
    InterfaceConstIterator endItr = interfaceConstEnd();
    for (InterfaceConstIterator itr = interfaceConstBegin(); itr != endItr; ++itr) {
        const long interfaceId = itr.getId();
        if (isInterfaceOrphan(interfaceId)) {
            orphanInterfaces.push_back(interfaceId);
        }
    }
 
    return orphanInterfaces;
}

bool PatchKernel::deleteOrphanInterfaces()
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    std::vector<long> list = findOrphanInterfaces();
    deleteInterfaces(list);
 
    return true;
}

bool PatchKernel::isInterfaceOrphan(long id) const
{
    const Interface &interface = getInterface(id);
 
    // Interface is not orphan if it has an owner or a neighbour
    long ownerId = interface.getOwner();
    long neighId = interface.getNeigh();
    if (ownerId >= 0 && neighId >= 0) {
        return false;
    }
 
    // Interface is not orphan if it's on a border face of an internal cell
    if (ownerId >= 0 && neighId < 0) {
        const Cell &owner = getCell(ownerId);
        if (owner.isInterior()) {
            return false;
        }
    }
 
    return true;
}

long PatchKernel::countFaces() const
{
    double nFaces = 0;
    for (const Cell &cell : m_cells) {
        int nCellFaces = cell.getFaceCount();
        for (int i = 0; i < nCellFaces; ++i) {
            if (!cell.isFaceBorder(i)) {
                nFaces += 1. / (cell.getAdjacencyCount(i) + 1);
            } else {
                nFaces += 1.;
            }
        }
    }
 
    return ((long) round(nFaces));
}

long PatchKernel::countFreeFaces() const
{
    return countBorderFaces();
}

long PatchKernel::countBorderFaces() const
{
    long nBorderFaces = 0;
    for (const Cell &cell : m_cells) {
        int nCellFaces = cell.getFaceCount();
        for (int i = 0; i < nCellFaces; ++i) {
            if (cell.isFaceBorder(i)) {
                ++nBorderFaces;
            }
        }
    }
 
    return nBorderFaces;
}

void PatchKernel::dumpVertices(std::ostream &stream) const
{
    // Dump kernel
    m_vertices.dumpKernel(stream);
 
    // Dump vertices
    for (const Vertex &vertex : m_vertices) {
        utils::binary::write(stream, vertex.getId());
#if BITPIT_ENABLE_MPI==1
        utils::binary::write(stream, getVertexOwner(vertex.getId()));
#else
        int dummyOwner = 0;
        utils::binary::write(stream, dummyOwner);
#endif
 
        const std::array<double, 3> &coords = vertex.getCoords();
        utils::binary::write(stream, coords[0]);
        utils::binary::write(stream, coords[1]);
        utils::binary::write(stream, coords[2]);
    }
 
    // Dump ghost/internal subdivision
#if BITPIT_ENABLE_MPI==1
    utils::binary::write(stream, m_firstGhostVertexId);
    utils::binary::write(stream, m_lastInternalVertexId);
#else
    utils::binary::write(stream, Vertex::NULL_ID);
    utils::binary::write(stream, Vertex::NULL_ID);
#endif
}

void PatchKernel::restoreVertices(std::istream &stream)
{
    // Restore kernel
    m_vertices.restoreKernel(stream);
 
    // Enable manual adaption
    AdaptionMode previousAdaptionMode = getAdaptionMode();
    setAdaptionMode(ADAPTION_MANUAL);
 
    // Restore vertices
    long nVertices = m_vertices.size();
    for (long i = 0; i < nVertices; ++i) {
        long id;
        utils::binary::read(stream, id);
 
#if BITPIT_ENABLE_MPI==1
        int owner;
        utils::binary::read(stream, owner);
#else
        int dummyOwner;
        utils::binary::read(stream, dummyOwner);
#endif
 
        std::array<double, 3> coords;
        utils::binary::read(stream, coords[0]);
        utils::binary::read(stream, coords[1]);
        utils::binary::read(stream, coords[2]);
 
#if BITPIT_ENABLE_MPI==1
        restoreVertex(std::move(coords), owner, id);
#else
        restoreVertex(std::move(coords), id);
#endif
    }
 
    // Restore ghost/internal subdivision
#if BITPIT_ENABLE_MPI==1
    utils::binary::read(stream, m_firstGhostVertexId);
    utils::binary::read(stream, m_lastInternalVertexId);
#else
    long dummyFirstGhostVertexId;
    long dummyLastInternalVertexId;
    utils::binary::read(stream, dummyFirstGhostVertexId);
    utils::binary::read(stream, dummyLastInternalVertexId);
#endif
 
    // Restore previous adaption mode
    setAdaptionMode(previousAdaptionMode);
}

void PatchKernel::dumpCells(std::ostream &stream) const
{
    // Dump kernel
    m_cells.dumpKernel(stream);
 
    // Dump cells
    for (const Cell &cell: m_cells) {
        utils::binary::write(stream, cell.getId());
        utils::binary::write(stream, cell.getPID());
        utils::binary::write(stream, cell.getType());
 
#if BITPIT_ENABLE_MPI==1
        utils::binary::write(stream, getCellOwner(cell.getId()));
        utils::binary::write(stream, getCellHaloLayer(cell.getId()));
#else
        int dummyOwner = 0;
        utils::binary::write(stream, dummyOwner);
 
        int dummyHaloLayer = 0;
        utils::binary::write(stream, dummyHaloLayer);
#endif
 
        int cellConnectSize = cell.getConnectSize();
        utils::binary::write(stream, cellConnectSize);
 
        const long *cellConnect = cell.getConnect();
        for (int i = 0; i < cellConnectSize; ++i) {
            utils::binary::write(stream, cellConnect[i]);
        }
    }
 
    // Dump ghost/internal subdivision
#if BITPIT_ENABLE_MPI==1
    utils::binary::write(stream, m_firstGhostCellId);
    utils::binary::write(stream, m_lastInternalCellId);
#else
    utils::binary::write(stream, Cell::NULL_ID);
    utils::binary::write(stream, Cell::NULL_ID);
#endif
}

void PatchKernel::restoreCells(std::istream &stream)
{
    // Restore kernel
    m_cells.restoreKernel(stream);
 
    // Enable manual adaption
    AdaptionMode previousAdaptionMode = getAdaptionMode();
    setAdaptionMode(ADAPTION_MANUAL);
 
    // Restore cells
    long nCells = m_cells.size();
    for (long i = 0; i < nCells; ++i) {
        long id;
        utils::binary::read(stream, id);
 
        int PID;
        utils::binary::read(stream, PID);
 
        ElementType type;
        utils::binary::read(stream, type);
 
#if BITPIT_ENABLE_MPI==1
        int owner;
        utils::binary::read(stream, owner);
 
        int haloLayer;
        utils::binary::read(stream, haloLayer);
#else
        int dummyOwner;
        utils::binary::read(stream, dummyOwner);
 
        int dummyHaloLayer;
        utils::binary::read(stream, dummyHaloLayer);
#endif
 
        int cellConnectSize;
        utils::binary::read(stream, cellConnectSize);
 
        std::unique_ptr<long[]> cellConnect = std::unique_ptr<long[]>(new long[cellConnectSize]);
        for (int k = 0; k < cellConnectSize; ++k) {
            utils::binary::read(stream, cellConnect[k]);
        }
 
        CellIterator iterator;
#if BITPIT_ENABLE_MPI==1
        iterator = restoreCell(type, std::move(cellConnect), owner, haloLayer, id);
#else
        iterator = restoreCell(type, std::move(cellConnect), id);
#endif
        iterator->setPID(PID);
    }
 
    // Restore ghost/internal subdivision
#if BITPIT_ENABLE_MPI==1
    utils::binary::read(stream, m_firstGhostCellId);
    utils::binary::read(stream, m_lastInternalCellId);
#else
    long dummyFirstGhostCellId;
    long dummyLastInternalCellId;
    utils::binary::read(stream, dummyFirstGhostCellId);
    utils::binary::read(stream, dummyLastInternalCellId);
#endif
 
    // Update adjacencies
    updateAdjacencies();
 
    // Restore previous adaption mode
    setAdaptionMode(previousAdaptionMode);
}

void PatchKernel::dumpInterfaces(std::ostream &stream) const
{
    // Early return if the interfaces are not build
    if (getInterfacesBuildStrategy() == INTERFACES_NONE) {
        return;
    }
 
    // Dump kernel
    m_interfaces.dumpKernel(stream);
 
    // Dump interfaces
    for (const Interface &interface : getInterfaces()) {
        utils::binary::write(stream, interface.getId());
 
        utils::binary::write(stream, interface.getOwner());
        utils::binary::write(stream, interface.getOwnerFace());
 
        long neighId = interface.getNeigh();
        utils::binary::write(stream, interface.getNeigh());
        if (neighId >= 0) {
            utils::binary::write(stream, interface.getNeighFace());
        }
 
        utils::binary::write(stream, interface.getPID());
    }
}

void PatchKernel::restoreInterfaces(std::istream &stream)
{
    // Early return if the interfaces are not build
    if (getInterfacesBuildStrategy() == INTERFACES_NONE) {
        return;
    }
 
    // Interfaces need up-to-date adjacencies
    updateAdjacencies();
 
    // Restore kernel
    m_interfaces.restoreKernel(stream);
 
    // Enable manual adaption
    AdaptionMode previousAdaptionMode = getAdaptionMode();
    setAdaptionMode(ADAPTION_MANUAL);
 
    // Restore interfaces
    long nInterfaces = m_interfaces.size();
    for (long n = 0; n < nInterfaces; ++n) {
        long interfaceId;
        utils::binary::read(stream, interfaceId);
 
        long ownerId;
        int ownerFace;
        utils::binary::read(stream, ownerId);
        utils::binary::read(stream, ownerFace);
        Cell *owner = &(m_cells.at(ownerId));
 
        long neighId;
        int neighFace;
        utils::binary::read(stream, neighId);
        Cell *neigh;
        if (neighId >= 0) {
            utils::binary::read(stream, neighFace);
            neigh = &(m_cells.at(neighId));
        } else {
            neighFace = -1;
            neigh = nullptr;
        }
 
        int pid;
        utils::binary::read(stream, pid);
 
        InterfaceIterator interfaceIterator = buildCellInterface(owner, ownerFace, neigh, neighFace, interfaceId);
        interfaceIterator->setPID(pid);
    }
 
    // Interfaces are now updated
    unsetCellAlterationFlags(FLAG_INTERFACES_DIRTY);
    m_alteredInterfaces.clear();
 
    // Restore previous adaption mode
    setAdaptionMode(previousAdaptionMode);
}

std::vector<adaption::Info> PatchKernel::_adaptionPrepare(bool trackAdaption)
{
    BITPIT_UNUSED(trackAdaption);
 
    return std::vector<adaption::Info>();
}

std::vector<adaption::Info> PatchKernel::_adaptionAlter(bool trackAdaption)
{
    BITPIT_UNUSED(trackAdaption);
 
    assert(false && "The patch needs to implement _adaptionAlter");
 
    return std::vector<adaption::Info>();
}

void PatchKernel::_adaptionCleanup()
{
}

bool PatchKernel::_markCellForRefinement(long id)
{
    BITPIT_UNUSED(id);
 
    return false;
}

bool PatchKernel::_markCellForCoarsening(long id)
{
    BITPIT_UNUSED(id);
 
    return false;
}

bool PatchKernel::_resetCellAdaptionMarker(long id)
{
    BITPIT_UNUSED(id);
 
    return false;
}

adaption::Marker PatchKernel::_getCellAdaptionMarker(long id)
{
    BITPIT_UNUSED(id);
 
    return adaption::MARKER_UNDEFINED;
}

bool PatchKernel::_enableCellBalancing(long id, bool enabled)
{
    BITPIT_UNUSED(id);
    BITPIT_UNUSED(enabled);
 
    return false;
}

bool PatchKernel::sortVertices()
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    // Sort internal vertices
    if (m_nInternalVertices > 0) {
        m_vertices.sortBefore(m_lastInternalVertexId, true);
        updateLastInternalVertexId();
    }
 
#if BITPIT_ENABLE_MPI==1
    // Sort ghost cells
    if (m_nGhostVertices > 0) {
        m_vertices.sortAfter(m_firstGhostVertexId, true);
        updateFirstGhostVertexId();
    }
#endif
 
    // Synchronize storage
    m_vertices.sync();
 
    return true;
}

bool PatchKernel::sortCells()
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    // Sort internal cells
    if (m_nInternalCells > 0) {
        m_cells.sortBefore(m_lastInternalCellId, true);
        updateLastInternalCellId();
    }
 
#if BITPIT_ENABLE_MPI==1
    // Sort ghost cells
    if (m_nGhostCells > 0) {
        m_cells.sortAfter(m_firstGhostCellId, true);
        updateFirstGhostCellId();
    }
#endif
 
    // Synchronize storage
    m_cells.sync();
 
    return true;
}

bool PatchKernel::sortInterfaces()
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    m_interfaces.sort();
 
    m_interfaces.sync();
 
    return true;
}

bool PatchKernel::sort()
{
    bool status = sortVertices();
    status |= sortCells();
    status |= sortInterfaces();
 
    return status;
}

bool PatchKernel::squeezeVertices()
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    m_vertices.squeeze();
 
    m_vertices.sync();
 
    return true;
}

bool PatchKernel::squeezeCells()
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    m_cells.squeeze();
 
    m_cells.sync();
 
    return true;
}

bool PatchKernel::squeezeInterfaces()
{
    if (getAdaptionMode() != ADAPTION_MANUAL) {
        return false;
    }
 
    m_interfaces.squeeze();
 
    m_interfaces.sync();
 
    return true;
}

bool PatchKernel::squeeze()
{
    bool status = true;
 
    // Alteration flags
    if (m_alteredCells.empty()) {
        AlterationFlagsStorage().swap(m_alteredCells);
    }
 
    if (m_alteredInterfaces.empty()) {
        AlterationFlagsStorage().swap(m_alteredInterfaces);
    }
 
    // Vertices
    status |= squeezeVertices();
 
    // Cells
    status |= squeezeCells();
 
    // Interfaces
    status |= squeezeInterfaces();
 
    return status;
}

std::array<double, 3> PatchKernel::evalCellCentroid(long id) const
{
    const Cell &cell = getCell(id);
 
    return evalElementCentroid(cell);
}

void PatchKernel::evalCellBoundingBox(long id, std::array<double,3> *minPoint, std::array<double,3> *maxPoint) const
{
    const Cell &cell = getCell(id);
 
    evalElementBoundingBox(cell, minPoint, maxPoint);
}

ConstProxyVector<std::array<double, 3>> PatchKernel::getCellVertexCoordinates(long id) const
{
    // Get the vertices ids
    const Cell &cell = getCell(id);
    ConstProxyVector<long> cellVertexIds = cell.getVertexIds();
    const int nCellVertices = cellVertexIds.size();
 
    // Build the proxy vector with the coordinates
    ConstProxyVector<std::array<double, 3>> vertexCoordinates(ConstProxyVector<std::array<double, 3>>::INTERNAL_STORAGE, nCellVertices);
    ConstProxyVector<std::array<double, 3>>::storage_pointer storage = vertexCoordinates.storedData();
    getVertexCoords(cellVertexIds.size(), cellVertexIds.data(), storage);
 
    return vertexCoordinates;
}

void PatchKernel::getCellVertexCoordinates(long id, std::unique_ptr<std::array<double, 3>[]> *coordinates) const
{
    const Cell &cell = getCell(id);
 
    getElementVertexCoordinates(cell, coordinates);
}

void PatchKernel::getCellVertexCoordinates(long id, std::array<double, 3> *coordinates) const
{
    const Cell &cell = getCell(id);
 
    getElementVertexCoordinates(cell, coordinates);
}

std::array<double, 3> PatchKernel::evalInterfaceCentroid(long id) const
{
    const Interface &interface = getInterface(id);
 
    return evalElementCentroid(interface);
}

void PatchKernel::evalInterfaceBoundingBox(long id, std::array<double,3> *minPoint, std::array<double,3> *maxPoint) const
{
    const Interface &interface = getInterface(id);
 
    evalElementBoundingBox(interface, minPoint, maxPoint);
}

ConstProxyVector<std::array<double, 3>> PatchKernel::getInterfaceVertexCoordinates(long id) const
{
    // Get the vertices ids
    const Interface &interface = getInterface(id);
    ConstProxyVector<long> interfaceVertexIds = interface.getVertexIds();
    const int nInterfaceVertices = interfaceVertexIds.size();
 
    // Build the proxy vector with the coordinates
    ConstProxyVector<std::array<double, 3>> vertexCoordinates(ConstProxyVector<std::array<double, 3>>::INTERNAL_STORAGE, nInterfaceVertices);
    ConstProxyVector<std::array<double, 3>>::storage_pointer storage = vertexCoordinates.storedData();
    getVertexCoords(interfaceVertexIds.size(), interfaceVertexIds.data(), storage);
 
    return vertexCoordinates;
}

void PatchKernel::getInterfaceVertexCoordinates(long id, std::unique_ptr<std::array<double, 3>[]> *coordinates) const
{
    const Interface &interface = getInterface(id);
 
    getElementVertexCoordinates(interface, coordinates);
}

void PatchKernel::getInterfaceVertexCoordinates(long id, std::array<double, 3> *coordinates) const
{
    const Interface &interface = getInterface(id);
 
    getElementVertexCoordinates(interface, coordinates);
}

bool PatchKernel::intersectInterfacePlane(long id, std::array<double, 3> P, std::array<double, 3> nP,
                                          std::array<std::array<double, 3>, 2> *intersection,
                                          std::vector<std::array<double, 3>> *polygon) const
{
        double distanceTolerance = getTol();
 
    const Interface &interface = getInterface(id);
    ElementType interfaceType = interface.getType();
    switch (interfaceType)
    {
 
    case ElementType::LINE:
    {
        const std::array<double, 3> &P0 = getVertexCoords(interface.getVertexId(0));
        const std::array<double, 3> &P1 = getVertexCoords(interface.getVertexId(1));
 
        std::array<double, 3> x;
        bool intersectionFound = CGElem::intersectSegmentPlane(P0, P1, P, nP, x, distanceTolerance);
 
        if (norm2(x - P0) < distanceTolerance) return false;
        if (norm2(x - P1) < distanceTolerance) return false;
 
        if (!intersectionFound) return false;
 
        (*intersection)[0] = x;
        (*intersection)[1] = x;
                if (dotProduct(P - P0, nP) > 0.) {
            polygon->push_back(P0);
            polygon->push_back(x);
                        return true;
                }
        polygon->push_back(x);
        polygon->push_back(P1);
                return true;
    }
 
    case ElementType::TRIANGLE:
    {
        const std::array<double, 3> &P0 = getVertexCoords(interface.getVertexId(0));
        const std::array<double, 3> &P1 = getVertexCoords(interface.getVertexId(1));
        const std::array<double, 3> &P2 = getVertexCoords(interface.getVertexId(2));
            
        return CGElem::intersectPlaneTriangle(P0, P1, P2, P, nP,
                       *intersection, *polygon, distanceTolerance);
    }
 
    case ElementType::PIXEL:
    {
        // evaluate pixel normal
        std::unique_ptr<std::array<double, 3>[]> coordinates;
        getElementVertexCoordinates(interface, &coordinates);
 
        return CGElem::intersectPlanePixel(P, nP, coordinates.get(),
                       *intersection, *polygon, distanceTolerance);
    }
 
    default:
    {
 
                assert(interfaceType != ElementType::UNDEFINED);
                const Reference2DElementInfo &interfaceInfo = static_cast<const Reference2DElementInfo &>(ReferenceElementInfo::getInfo(interfaceType));
 
            ConstProxyVector<long> vertexIds = interface.getVertexIds();
                int nVertices = vertexIds.size();
 
                std::vector<std::array<double, 3>> coordinates(nVertices);
                for (int i = 0; i < nVertices; ++i) {
                    long vertexId = vertexIds[interfaceInfo.getCCWOrderedVertex(i)];
                    coordinates[i] = getVertex(vertexId).getCoords();
                }
 
                std::vector<std::array<std::array<double, 3>, 2>> intersections;
                std::vector< std::vector<std::array<double, 3>>> polygons;
        bool intersectionFound = CGElem::intersectPlanePolygon(P, nP, nVertices, coordinates.data(),
                                 intersections, polygons, distanceTolerance);
 
        if (intersections.size() > 1) {
            throw std::runtime_error("Intersection of concave interfaces is not supported.");
                }
 
        (*intersection) = intersections[0];
        (*polygon)      = polygons[0];
 
        return intersectionFound;
    }
 
    }
        return false;
}

std::array<double, 3> PatchKernel::evalElementCentroid(const Element &element) const
{
    ConstProxyVector<long> elementVertexIds = element.getVertexIds();
    std::size_t nCellVertices = elementVertexIds.size();
    BITPIT_CREATE_WORKSPACE(vertexCoordinates, std::array<double BITPIT_COMMA 3>, nCellVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    getVertexCoords(nCellVertices, elementVertexIds.data(), vertexCoordinates);
 
    return element.evalCentroid(vertexCoordinates);
}

void PatchKernel::evalElementBoundingBox(const Element &element, std::array<double,3> *minPoint, std::array<double,3> *maxPoint) const
{
    ElementType elementType = element.getType();
    switch (elementType)
    {
 
    case ElementType::VERTEX:
    {
        const long *elementConnect = element.getConnect();
        *minPoint = getVertexCoords(elementConnect[0]);
        *maxPoint = *minPoint;
        break;
    }
 
    case ElementType::PIXEL:
    {
        const long *elementConnect = element.getConnect();
 
        const std::array<double, 3> &vertexCoord_0 = getVertexCoords(elementConnect[0]);
        const std::array<double, 3> &vertexCoord_3 = getVertexCoords(elementConnect[3]);
 
        for (int d = 0; d < 3; ++d) {
            (*minPoint)[d] = std::min(vertexCoord_0[d], vertexCoord_3[d]);
            (*maxPoint)[d] = std::max(vertexCoord_0[d], vertexCoord_3[d]);
        }
 
        break;
    }
 
    case ElementType::VOXEL:
    {
        const long *elementConnect = element.getConnect();
 
        const std::array<double, 3> &vertexCoord_0 = getVertexCoords(elementConnect[0]);
        const std::array<double, 3> &vertexCoord_7 = getVertexCoords(elementConnect[7]);
 
        for (int d = 0; d < 3; ++d) {
            (*minPoint)[d] = std::min(vertexCoord_0[d], vertexCoord_7[d]);
            (*maxPoint)[d] = std::max(vertexCoord_0[d], vertexCoord_7[d]);
        }
 
        break;
    }
 
    default:
    {
        ConstProxyVector<long> elementVertexIds = element.getVertexIds();
        const int nElementVertices = elementVertexIds.size();
 
        *minPoint = getVertexCoords(elementVertexIds[0]);
        *maxPoint = *minPoint;
        for (int i = 1; i < nElementVertices; ++i) {
            const std::array<double, 3> &vertexCoord = getVertexCoords(elementVertexIds[i]);
            for (int d = 0; d < 3; ++d) {
                (*minPoint)[d] = std::min(vertexCoord[d], (*minPoint)[d]);
                (*maxPoint)[d] = std::max(vertexCoord[d], (*maxPoint)[d]);
            }
        }
 
        break;
    }
 
    }
}

long PatchKernel::locatePoint(double x, double y, double z) const
{
    return locatePoint({{x, y, z}});
}

bool PatchKernel::isSameFace(const Cell &cell_A, int face_A, const Cell &cell_B, int face_B) const
{
    if (cell_A.getFaceType(face_A) != cell_B.getFaceType(face_B)) {
        return false;
    }
 
    int nFaceVertices = cell_A.getFaceVertexCount(face_A);
 
    std::vector<long> faceVertexIds_A(nFaceVertices);
    std::vector<long> faceVertexIds_B(nFaceVertices);
    for (int k = 0; k < nFaceVertices; ++k) {
        faceVertexIds_A[k] = cell_A.getFaceVertexId(face_A, k);
        faceVertexIds_B[k] = cell_B.getFaceVertexId(face_B, k);
    }
 
    std::sort(faceVertexIds_A.begin(), faceVertexIds_A.end());
    std::sort(faceVertexIds_B.begin(), faceVertexIds_B.end());
 
    return (faceVertexIds_A == faceVertexIds_B);
}

PatchKernel::AdjacenciesBuildStrategy PatchKernel::getAdjacenciesBuildStrategy() const
{
    return m_adjacenciesBuildStrategy;
}

void PatchKernel::setAdjacenciesBuildStrategy(AdjacenciesBuildStrategy status)
{
    m_adjacenciesBuildStrategy = status;
}

bool PatchKernel::areAdjacenciesDirty(bool global) const
{
    if (getAdjacenciesBuildStrategy() == ADJACENCIES_NONE) {
        return false;
    }
 
    bool areDirty = false;
    for (const auto &entry : m_alteredCells) {
        AlterationFlags cellAlterationFlags = entry.second;
        if (testAlterationFlags(cellAlterationFlags, FLAG_ADJACENCIES_DIRTY)) {
            areDirty = true;
            break;
        }
    }
 
#if BITPIT_ENABLE_MPI==1
    if (global && isPartitioned()) {
        const auto &communicator = getCommunicator();
        MPI_Allreduce(MPI_IN_PLACE, &areDirty, 1, MPI_C_BOOL, MPI_LOR, communicator);
    }
#else
    BITPIT_UNUSED(global);
#endif
 
    return areDirty;
}

void PatchKernel::buildAdjacencies()
{
    initializeAdjacencies(ADJACENCIES_AUTOMATIC);
}

void PatchKernel::initializeAdjacencies(AdjacenciesBuildStrategy strategy)
{
    // Get current build stategy
    AdjacenciesBuildStrategy currentStrategy = getAdjacenciesBuildStrategy();
 
    // Early return if we don't need adjacencies
    if (strategy == ADJACENCIES_NONE) {
        if (currentStrategy != ADJACENCIES_NONE) {
            destroyAdjacencies();
        }
 
        return;
    }
 
    // Update the build strategy
    if (currentStrategy != strategy) {
        setAdjacenciesBuildStrategy(strategy);
    }
 
    // Reset adjacencies
    resetAdjacencies();
 
    // Update the adjacencies
    updateAdjacencies();
}

void PatchKernel::updateAdjacencies(bool forcedUpdated)
{
    // Early return if adjacencies are not built
    AdjacenciesBuildStrategy currentStrategy = getAdjacenciesBuildStrategy();
    if (currentStrategy == ADJACENCIES_NONE) {
        return;
    }
 
    // Check if the adjacencies are dirty
    bool adjacenciesDirty = areAdjacenciesDirty();
    if (!adjacenciesDirty && !forcedUpdated) {
        return;
    }
 
    // Update adjacencies
    if (adjacenciesDirty) {
        // Prune stale adjacencies
        pruneStaleAdjacencies();
 
        // Update adjacencies
        _updateAdjacencies();
 
        // Adjacencies are now updated
        unsetCellAlterationFlags(FLAG_ADJACENCIES_DIRTY);
    } else {
        initializeAdjacencies(currentStrategy);
    }
}

void PatchKernel::destroyAdjacencies()
{
    // Early return if no adjacencies have been built
    if (getAdjacenciesBuildStrategy() == ADJACENCIES_NONE) {
        return;
    }
 
    // Destroy the adjacencies
    _resetAdjacencies(true);
 
    // Clear list of cells with dirty adjacencies
    unsetCellAlterationFlags(FLAG_ADJACENCIES_DIRTY);
 
    // Set adjacencies status
    setAdjacenciesBuildStrategy(ADJACENCIES_NONE);
}

void PatchKernel::resetAdjacencies()
{
    // Early return if no adjacencies have been built
    if (getAdjacenciesBuildStrategy() == ADJACENCIES_NONE) {
        return;
    }
 
    // Reset adjacencies
    _resetAdjacencies(false);
 
    // All cells have now dirty adjacencies
    setCellAlterationFlags(FLAG_ADJACENCIES_DIRTY);
}

void PatchKernel::_resetAdjacencies(bool release)
{
    for (Cell &cell : m_cells) {
        cell.resetAdjacencies(!release);
    }
}

void PatchKernel::pruneStaleAdjacencies()
{
    // Early return if adjacencies are not built
    AdjacenciesBuildStrategy currentStrategy = getAdjacenciesBuildStrategy();
    if (currentStrategy == ADJACENCIES_NONE) {
        return;
    }
 
    // Remove stale adjacencies
    for (const auto &entry : m_alteredCells) {
        AlterationFlags cellAlterationFlags = entry.second;
        if (!testAlterationFlags(cellAlterationFlags, FLAG_ADJACENCIES_DIRTY)) {
            continue;
        }
 
        long cellId = entry.first;
        Cell &cell = m_cells.at(cellId);
        if (cell.getAdjacencyCount() == 0) {
            continue;
        }
 
        int nCellFaces = cell.getFaceCount();
        for (int face = nCellFaces - 1; face >= 0; --face) {
            long *faceAdjacencies = cell.getAdjacencies(face);
            int nFaceAdjacencies = cell.getAdjacencyCount(face);
            for (int i = nFaceAdjacencies - 1; i >= 0; --i) {
                long adjacency = faceAdjacencies[i];
                if (!testCellAlterationFlags(adjacency, FLAG_DELETED)) {
                    continue;
                }
 
                cell.deleteAdjacency(face, i);
            }
        }
    }
}

void PatchKernel::_updateAdjacencies()
{
    int dimension = getDimension();
 
    // Check if multiple adjacencies are allowed
    //
    // On a three-dimensional patch each internal face is shared between two,
    // and only two, half-faces. On lower dimension patches one internal face
    // may be shared among multiple half-faces (this happens with non-manifold
    // patches).
    bool multipleMatchesAllowed = (dimension < 3);
 
    // Define matching windings
    //
    // If faces can be shared only between two half-faces, there is a one-to-one
    // correspondence between the half-faces. Give one half-faces, its matching
    // can be found looking for an half-faces with the same vertices but with
    // reverse vertex winding order.
    //
    // If faces can be shared among more than two half-faces, it's not anymore
    // possible to define a unique winding order for the vertices of the faces.
    // There is no more a one-to-one correspondence between pairs of half-face.
    // In this case, when looking for half-face correspondence, it is necessary
    // to look for half-faces with reverse vertex winding order and also for
    // half-faces with the same vertex winding order.
    std::vector<CellHalfFace::Winding> matchingWindings;
    matchingWindings.push_back(CellHalfFace::WINDING_REVERSE);
    if (multipleMatchesAllowed) {
        matchingWindings.push_back(CellHalfFace::WINDING_NATURAL);
    }
 
    // Count cells with dirty adjacencies
    long nDirtyAdjacenciesCells = 0;
    for (const auto &entry : m_alteredCells) {
        AlterationFlags cellAlterationFlags = entry.second;
        if (!testAlterationFlags(cellAlterationFlags, FLAG_ADJACENCIES_DIRTY)) {
            continue;
        }
 
        ++nDirtyAdjacenciesCells;
    }
 
    // Get the vertices of cells with dirty adjacencies
    //
    // The list is only needed if multiple matches are allowed and there are
    // cells whose adjacencies don't need to be updated.
    std::unique_ptr<PiercedStorage<bool, long>> dirtyAdjacenciesVertices;
    if (multipleMatchesAllowed && (nDirtyAdjacenciesCells != getCellCount())) {
        dirtyAdjacenciesVertices = std::unique_ptr<PiercedStorage<bool, long>>(new PiercedStorage<bool, long>(1, &m_vertices));
        dirtyAdjacenciesVertices->fill(false);
        for (const auto &entry : m_alteredCells) {
            AlterationFlags cellAlterationFlags = entry.second;
            if (!testAlterationFlags(cellAlterationFlags, FLAG_ADJACENCIES_DIRTY)) {
                continue;
            }
 
            long cellId = entry.first;
            const Cell &cell = m_cells.at(cellId);
            ConstProxyVector<long> cellVertexIds = cell.getVertexIds();
            for (long vertexId : cellVertexIds) {
                dirtyAdjacenciesVertices->at(vertexId) = true;
            }
        }
    }
 
    // List the cells that needs to be processed
    //
    // The list should contain the cells whose adjacencies need to be updated
    // and their neighbours candidates.
    //
    // Cells whose adjacencies needs to be updated are cells having the dirty
    // adjaceny flag set.
    //
    // Neighbours candidates are all border cells and, if multiple matches are
    // allowed, cells that share vertices with a cell with dirty adjacencies.
    std::vector<Cell *> processList;
    processList.reserve(nDirtyAdjacenciesCells);
 
    std::size_t nMaxHalfFaces = 0;
    for (Cell &cell : m_cells) {
        // Get cell information
        long cellId = cell.getId();
        int nCellFaces = cell.getFaceCount();
 
        // Add cells with dirty adjacencies
        if (testCellAlterationFlags(cellId, FLAG_ADJACENCIES_DIRTY)) {
            nMaxHalfFaces += nCellFaces;
            processList.push_back(&cell);
            continue;
        }
 
        // Add border cells
        bool isBorderCell = false;
        for (int face = 0; face < nCellFaces; face++) {
            if (cell.isFaceBorder(face)) {
                isBorderCell = true;
                continue;
            }
        }
 
        if (isBorderCell) {
            nMaxHalfFaces += nCellFaces;
            processList.push_back(&cell);
            continue;
        }
 
        // Add cell that share vertices with a cell with dirty adjacencies
        //
        // For performace reasons, the check is a bit rough. All cells that
        // share a number of vertices at least equal to the dimension af the
        // patch will be added to the process list. This will add also cells
        // that are not neighbour candidates, but it's faster to add some false
        // positives, rather trying to filter them out.
        if (multipleMatchesAllowed) {
            ConstProxyVector<long> cellVertexIds = cell.getVertexIds();
 
            int nCelldirtyAdjacenciesVertices = 0;
            bool futureNeighbourCandidate = false;
            for (long vertexId : cellVertexIds) {
                if (dirtyAdjacenciesVertices->at(vertexId)) {
                    ++nCelldirtyAdjacenciesVertices;
                    if (nCelldirtyAdjacenciesVertices >= dimension) {
                        futureNeighbourCandidate = true;
                        break;
                    }
                }
            }
 
            if (futureNeighbourCandidate) {
                nMaxHalfFaces += nCellFaces;
                processList.push_back(&cell);
            }
        }
    }
 
    // Create the adjacencies
    std::unordered_set<CellHalfFace, CellHalfFace::Hasher> halfFaces;
    halfFaces.reserve(static_cast<std::size_t>(0.5 * nMaxHalfFaces));
 
    std::vector<std::pair<Cell *, int>> matchingAdjacencies;
    for (Cell *cell : processList) {
        long cellId = cell->getId();
        const int nCellFaces = cell->getFaceCount();
        bool areCellAdjacenciesDirty = testCellAlterationFlags(cellId, FLAG_ADJACENCIES_DIRTY);
        for (int face = 0; face < nCellFaces; face++) {
            // Generate the half-face
            CellHalfFace halfFace(*cell, face);
 
            // Find matching half-face
            auto matchingHalfFaceItr = halfFaces.end();
            for (CellHalfFace::Winding winding : matchingWindings) {
                // Set winding order
                halfFace.setWinding(winding);
 
                // Search for matching half-face
                matchingHalfFaceItr = halfFaces.find(halfFace);
                if (matchingHalfFaceItr != halfFaces.end()) {
                    break;
                }
            }
 
            // Early return if no matching half-face has been found
            //
            // If no matching half-face has been found, we need to add the
            // current half-face to the list because it may match the face
            // of a cell that will be processed later.
            if (matchingHalfFaceItr == halfFaces.end()) {
                halfFace.setWinding(CellHalfFace::WINDING_NATURAL);
                halfFaces.insert(std::move(halfFace));
                continue;
            }
 
            // Get matching half-face information
            const CellHalfFace &matchingHalfFace = *matchingHalfFaceItr;
            Cell &matchingCell  = matchingHalfFace.getCell();
            long matchingCellId = matchingCell.getId();
            long matchingFace   = matchingHalfFace.getFace();
 
            // Only process cells with dirty adjacencies
            //
            // In order to limit the half faces that need to be stored, we are
            // processing at the same time cells with dirty adjacencies and
            // their neighbour candidates. This means we will also find matches
            // between faces of cell whose adjacensies are already up-to-date.
            bool areMatchingCellAdjacenciesDirty = testCellAlterationFlags(matchingCellId, FLAG_ADJACENCIES_DIRTY);
            if (!areCellAdjacenciesDirty && !areMatchingCellAdjacenciesDirty) {
                continue;
            }
 
            // Identifty matching adjacencies
            //
            // If multiple matches are allowed, in addition to the entry
            // related to the matching half face, we also need to add the
            // entries associated to the adjacencies already identified
            // for the half face.
            matchingAdjacencies.clear();
            matchingAdjacencies.emplace_back(std::make_pair<Cell *, int>(&matchingCell, matchingFace));
            if (multipleMatchesAllowed) {
                const int nMachingFaceNeighs = matchingCell.getAdjacencyCount(matchingFace);
                const long *machingFaceNeighs = matchingCell.getAdjacencies(matchingFace);
                for (int k = 0; k < nMachingFaceNeighs; ++k) {
                    long neighId = machingFaceNeighs[k];
                    if (neighId == cellId) {
                        continue;
                    }
 
                    Cell &neigh = m_cells.at(neighId);
                    int neighFace = findAdjoinNeighFace(matchingCell, matchingFace, neigh);
                    matchingAdjacencies.emplace_back(std::make_pair<Cell *, int>(&neigh, std::move(neighFace)));
                }
            } else {
                // When cell adjacencies are makred dirty this means that some
                // adjacencies are dirty not that all adjacencies are dirty.
                // The matchin cell may have no adjaceny associated to the
                // matching face or a signle adjacency that should match
                // the current cell.
                assert(matchingCell.getAdjacencyCount(matchingFace) == 0 || (matchingCell.getAdjacencyCount(matchingFace) == 1 && (*(matchingCell.getAdjacencies(matchingFace)) == cellId)));
            }
 
            // Create adjacencies
            for (const std::pair<Cell *, int> &matchingAdjacency : matchingAdjacencies) {
                Cell *adjacentCell = matchingAdjacency.first;
                long adjacentCellId = adjacentCell->getId();
                int adjacentFace = matchingAdjacency.second;
 
                cell->pushAdjacency(face, adjacentCellId);
                adjacentCell->pushAdjacency(adjacentFace, cellId);
            }
 
            // Remove the matching half-face from the list
            //
            // If multiple matches are allowed, we need to keep the half
            // face, because that entry may match the face of a cell that
            // will be processed later.
            if (!multipleMatchesAllowed) {
                halfFaces.erase(matchingHalfFaceItr);
            }
        }
    }
}

PatchKernel::InterfacesBuildStrategy PatchKernel::getInterfacesBuildStrategy() const
{
    return m_interfacesBuildStrategy;
}

void PatchKernel::setInterfacesBuildStrategy(InterfacesBuildStrategy status)
{
    if (status == INTERFACES_AUTOMATIC) {
        if (!isInterfaceAutoIndexingEnabled()) {
            throw std::runtime_error("Automatic build strategy requires auto-indexing.");
        }
    }
 
    m_interfacesBuildStrategy = status;
}

bool PatchKernel::areInterfacesDirty(bool global) const
{
    if (getInterfacesBuildStrategy() == INTERFACES_NONE) {
        return false;
    }
 
    bool areDirty = !m_alteredInterfaces.empty();
    if (!areDirty) {
        for (const auto &entry : m_alteredCells) {
            AlterationFlags cellAlterationFlags = entry.second;
            if (testAlterationFlags(cellAlterationFlags, FLAG_INTERFACES_DIRTY)) {
                areDirty = true;
                break;
            }
        }
    }
 
#if BITPIT_ENABLE_MPI==1
    if (global && isPartitioned()) {
        const auto &communicator = getCommunicator();
        MPI_Allreduce(MPI_IN_PLACE, &areDirty, 1, MPI_C_BOOL, MPI_LOR, communicator);
    }
#else
    BITPIT_UNUSED(global);
#endif
 
    return areDirty;
}

void PatchKernel::buildInterfaces()
{
    initializeInterfaces(INTERFACES_AUTOMATIC);
}

void PatchKernel::initializeInterfaces(InterfacesBuildStrategy strategy)
{
    // Interfaces need adjacencies
    if (getAdjacenciesBuildStrategy() == ADJACENCIES_NONE) {
        throw std::runtime_error ("Adjacencies are mandatory for building the interfaces.");
    }
 
    // Get current build stategy
    InterfacesBuildStrategy currentStrategy = getInterfacesBuildStrategy();
 
    // Early return if we don't need interfaces
    if (strategy == INTERFACES_NONE) {
        if (currentStrategy != INTERFACES_NONE) {
            destroyInterfaces();
        }
 
        return;
    }
 
    // Update the build strategy
    if (currentStrategy != strategy) {
        setInterfacesBuildStrategy(strategy);
    }
 
    // Reset interfaces
    resetInterfaces();
 
    // Update the interfaces
    updateInterfaces();
}

void PatchKernel::updateInterfaces(bool forcedUpdated)
{
    // Early return if interfaces are not built
    InterfacesBuildStrategy currentStrategy = getInterfacesBuildStrategy();
    if (currentStrategy == INTERFACES_NONE) {
        return;
    }
 
    // Check if the interfaces are dirty
    bool interfacesDirty = areInterfacesDirty();
    if (!interfacesDirty && !forcedUpdated) {
        return;
    }
 
    // Interfaces need up-to-date adjacencies
    assert(getAdjacenciesBuildStrategy() != ADJACENCIES_NONE);
    updateAdjacencies();
 
    // Update interfaces
    if (interfacesDirty) {
        // Enable manual adaption
        AdaptionMode previousAdaptionMode = getAdaptionMode();
        setAdaptionMode(ADAPTION_MANUAL);
 
        // Prune stale interfaces
        pruneStaleInterfaces();
 
        // Update interfaces
        _updateInterfaces();
 
        // Interfaces are now updated
        unsetCellAlterationFlags(FLAG_INTERFACES_DIRTY);
        m_alteredInterfaces.clear();
 
        // Restore previous adaption mode
        setAdaptionMode(previousAdaptionMode);
    } else {
        initializeInterfaces(currentStrategy);
    }
}

void PatchKernel::destroyInterfaces()
{
    // Early return if no interfaces have been built
    if (getInterfacesBuildStrategy() == INTERFACES_NONE) {
        return;
    }
 
    // Destroy the interfaces
    _resetInterfaces(true);
 
    // Clear list of cells with dirty interfaces
    unsetCellAlterationFlags(FLAG_INTERFACES_DIRTY);
 
    // Clear list of altered interfaces
    m_alteredInterfaces.clear();
 
    // Set interface build strategy
    setInterfacesBuildStrategy(INTERFACES_NONE);
}

void PatchKernel::pruneStaleInterfaces()
{
    // Early return if no interfaces have been built
    if (getInterfacesBuildStrategy() == INTERFACES_NONE) {
        return;
    }
 
    // Remove dangling interfaces from cells
    for (const auto &entry : m_alteredCells) {
        AlterationFlags cellAlterationFlags = entry.second;
        if (!testAlterationFlags(cellAlterationFlags, FLAG_INTERFACES_DIRTY)) {
            continue;
        }
 
        long cellId = entry.first;
        Cell &cell = m_cells.at(cellId);
        if (cell.getInterfaceCount() == 0) {
            continue;
        }
 
        int nCellFaces = cell.getFaceCount();
        for (int face = nCellFaces - 1; face >= 0; --face) {
            long *faceInterfaces = cell.getInterfaces(face);
            int nFaceInterfaces = cell.getInterfaceCount(face);
            for (int i = nFaceInterfaces - 1; i >= 0; --i) {
                long interfaceId = faceInterfaces[i];
                if (!testInterfaceAlterationFlags(interfaceId, FLAG_DANGLING)) {
                    continue;
                }
 
                // Delete the interface from the cell
                cell.deleteInterface(face, i);
            }
        }
    }
 
    // Delete dangling interfaces
    std::vector<long> danglingInterfaces;
    danglingInterfaces.reserve(m_alteredInterfaces.size());
    for (const auto &entry : m_alteredInterfaces) {
        AlterationFlags interfaceAlterationFlags = entry.second;
        if (!testAlterationFlags(interfaceAlterationFlags, FLAG_DANGLING)) {
            continue;
        }
 
        long interfaceId = entry.first;
        danglingInterfaces.push_back(interfaceId);
    }
    deleteInterfaces(danglingInterfaces);
}

void PatchKernel::_updateInterfaces()
{
    // Update interfaces
    //
    // Adjacencies and interfaces of a face are paired: the i-th face adjacency
    // corresponds to the i-th face interface. Moreover if we loop through the
    // adjacencies of a face, the adjacencies that have an interface are always
    // listed first. This meas that, to update the interfaces, we can count the
    // interfaces already associated to a face and loop only on the adjacencies
    // which have an index past the one of the last interface.
    //
    // On border faces of internal cells we need to build an interface, also
    // if there are no adjacencies.
    for (const auto &entry : m_alteredCells) {
        AlterationFlags cellAlterationFlags = entry.second;
        if (!testAlterationFlags(cellAlterationFlags, FLAG_INTERFACES_DIRTY)) {
            continue;
        }
 
        long cellId = entry.first;
        Cell &cell = m_cells.at(cellId);
        const int nCellFaces = cell.getFaceCount();
        for (int face = 0; face < nCellFaces; face++) {
            int nFaceInterfaces= cell.getInterfaceCount(face);
 
            bool isFaceBorder = cell.isFaceBorder(face);
            if (!isFaceBorder) {
                // Find the range of adjacencies that need an interface
                int updateEnd   = cell.getAdjacencyCount(face);
                int updateBegin = nFaceInterfaces;
                if (updateBegin == updateEnd) {
                    continue;
                }
 
                // Build an interface for every adjacency
                const long *faceAdjacencies = cell.getAdjacencies(face);
                for (int k = updateBegin; k < updateEnd; ++k) {
                    long neighId = faceAdjacencies[k];
                    Cell *neigh  = &m_cells[neighId];
 
                    int neighFace = findAdjoinNeighFace(cell, face, *neigh);
 
                    buildCellInterface(&cell, face, neigh, neighFace);
                }
            } else if (nFaceInterfaces == 0) {
                // Internal borderes need an interface
                buildCellInterface(&cell, face, nullptr, -1);
            }
        }
    }
}

PatchKernel::InterfaceIterator PatchKernel::buildCellInterface(Cell *cell_1, int face_1, Cell *cell_2, int face_2, long interfaceId)
{
    // Validate the input
    assert(cell_1);
    assert(face_1 >= 0);
 
    if (cell_2) {
        assert(face_2 >= 0);
    } else {
        assert(face_2 < 0);
    }
 
    // Get the id of the first cell
    long id_1 = cell_1->getId();
 
    // Get the id of the second cell
    long id_2 = Cell::NULL_ID;
    if (cell_2) {
        id_2 = cell_2->getId();
    }
 
    // Owner and neighbour of the interface
    //
    // The interface is owned by the cell that has only one adjacency, i.e.,
    // by the cell that owns the smallest of the two faces. If the faces
    // of both cells have the same size, the interface is owned by the cell
    // with the "lower fuzzy positioning". It is not necessary to have a
    // precise comparison, it's only necessary to define a repeatable order
    // between the two cells. It is therefore possible to use the "fuzzy"
    // cell comparison.
    bool cellOwnsInterface = true;
    if (cell_2) {
        if (cell_1->getAdjacencyCount(face_1) > 1) {
            cellOwnsInterface = false;
        } else if (cell_2->getAdjacencyCount(face_2) == 1) {
            assert(cell_1->getAdjacencyCount(face_1) == 1);
            cellOwnsInterface = CellFuzzyPositionLess(*this)(id_1, id_2);
        }
    }
 
    Cell *intrOwner;
    Cell *intrNeigh;
    long intrOwnerId;
    long intrNeighId;
    int intrOwnerFace;
    int intrNeighFace;
    if (cellOwnsInterface) {
        intrOwnerId   = id_1;
        intrOwner     = cell_1;
        intrOwnerFace = face_1;
 
        if (cell_2) {
            intrNeighId   = id_2;
            intrNeigh     = cell_2;
            intrNeighFace = face_2;
        } else {
            intrNeighId   = Cell::NULL_ID;
            intrNeigh     = nullptr;
            intrNeighFace = -1;
        }
    } else {
        intrOwnerId   = id_2;
        intrOwner     = cell_2;
        intrOwnerFace = face_2;
 
        intrNeighId   = id_1;
        intrNeigh     = cell_1;
        intrNeighFace = face_1;
    }
 
    // Connectivity of the interface
    ConstProxyVector<long> faceConnect = intrOwner->getFaceConnect(intrOwnerFace);
 
    int nInterfaceVertices = faceConnect.size();
    std::unique_ptr<long[]> interfaceConnect = std::unique_ptr<long[]>(new long[nInterfaceVertices]);
    for (int k = 0; k < nInterfaceVertices; ++k) {
        interfaceConnect[k] = faceConnect[k];
    }
 
    // Create the interface
    Interface *interface;
    InterfaceIterator interfaceIterator;
 
    ElementType interfaceType = intrOwner->getFaceType(intrOwnerFace);
    if (interfaceId < 0) {
        interfaceIterator = addInterface(interfaceType, std::move(interfaceConnect), interfaceId);
        interface = &(*interfaceIterator);
        interfaceId = interface->getId();
    } else {
        interfaceIterator = restoreInterface(interfaceType, std::move(interfaceConnect), interfaceId);
        interface = &(*interfaceIterator);
    }
 
    // Set owner and neighbour
    interface->setOwner(intrOwnerId, intrOwnerFace);
    if (intrNeighId >= 0) {
        interface->setNeigh(intrNeighId, intrNeighFace);
    }
 
    // Update owner and neighbour cell data
    //
    // Adjacencies and interfaces are paired: the i-th adjacency correspondes
    // to the i-th interface. Moreover if we loop through the adjacencies of
    // a face, the adjacencies that have an interface are always listed first.
    // When the interfaces are fully built, each adjacency is associated to an
    // interface, however this will not be true during the generation of the
    // interfaces. If the adjacencies already associated to an interface are
    // listed first, it will be easy to detect which interfaces are still to
    // be built: the missing interfaces are associated to the adjacencies which
    // have an index past the last interface.
    //
    // The above only matters if the neighbour cell exists, if there is no
    // neighbour there will be only one interface on that face, so there are
    // no ordering issues.
    intrOwner->pushInterface(intrOwnerFace, interfaceId);
    if (intrNeigh) {
        intrNeigh->pushInterface(intrNeighFace, interfaceId);
 
        // Fix adjacency order on the owner cell
        int ownerInterfaceIndex   = intrOwner->getInterfaceCount(intrOwnerFace) - 1;
        long ownerPairedAdjacency = intrOwner->getAdjacency(intrOwnerFace, ownerInterfaceIndex);
        if (ownerPairedAdjacency != intrNeighId) {
            int ownerPairedAdjacencyIndex = intrOwner->findAdjacency(intrOwnerFace, intrNeighId);
            assert(ownerPairedAdjacencyIndex >= 0);
            intrOwner->setAdjacency(intrOwnerFace, ownerInterfaceIndex, intrNeighId);
            intrOwner->setAdjacency(intrOwnerFace, ownerPairedAdjacencyIndex, ownerPairedAdjacency);
        }
 
        // Fix adjacency order on the neighbour cell
        int neighInterfaceIndex   = intrNeigh->getInterfaceCount(intrNeighFace) - 1;
        long neighPairedAdjacency = intrNeigh->getAdjacency(intrNeighFace, neighInterfaceIndex);
        if (neighPairedAdjacency != intrOwnerId) {
            int neighPairedAdjacencyIndex = intrNeigh->findAdjacency(intrNeighFace, intrOwnerId);
            assert(neighPairedAdjacencyIndex >= 0);
            intrNeigh->setAdjacency(intrNeighFace, neighInterfaceIndex, intrOwnerId);
            intrNeigh->setAdjacency(intrNeighFace, neighPairedAdjacencyIndex, neighPairedAdjacency);
        }
    }
 
    return interfaceIterator;
}

int PatchKernel::findAdjoinNeighFace(const Cell &cell, int cellFace, const Cell &neigh) const
{
    // Evaluate list of candidate faces
    //
    // The cells may be neighbours through multiple faces. Identify which face
    // matches the target one is quite expensive. Moreover, if the cells are
    // neighbours through a single face (which is the common case), the check
    // is not needed at all. Therefore, first we identify all the faces through
    // which the two cells are neighbours and then, if and only if there are
    // multiple candidates, we identify the face that matches the target one.
    long cellId = cell.getId();
    const int nNeighFaces = neigh.getFaceCount();
 
    int nCandidates = 0;
    BITPIT_CREATE_WORKSPACE(candidates, std::array<int BITPIT_COMMA 3>, nNeighFaces, ReferenceElementInfo::MAX_ELEM_FACES);
    for (int neighFace = 0; neighFace < nNeighFaces; neighFace++) {
        int nFaceAdjacencies = neigh.getAdjacencyCount(neighFace);
        const long *faceAdjacencies = neigh.getAdjacencies(neighFace);
        for (int k = 0; k < nFaceAdjacencies; ++k) {
            long geussId = faceAdjacencies[k];
            if (geussId == cellId) {
                (*candidates)[nCandidates] = neighFace;
                ++nCandidates;
                break;
            }
        }
    }
 
    // Select the face of the neighbour which adjoin the given face
    if (nCandidates == 1) {
        return (*candidates)[0];
    } else {
        for (int i = 0; i < nCandidates; ++i) {
            int candidateFace = (*candidates)[i];
            if (isSameFace(cell, cellFace, neigh, candidateFace)) {
                return candidateFace;
            }
        }
    }
 
    return -1;
}

bool PatchKernel::testCellAlterationFlags(long id, AlterationFlags flags) const
{
    return testElementAlterationFlags(id, flags, m_alteredCells);
}

PatchKernel::AlterationFlags PatchKernel::getCellAlterationFlags(long id) const
{
    return getElementAlterationFlags(id, m_alteredCells);
}

void PatchKernel::resetCellAlterationFlags(long id, AlterationFlags flags)
{
    resetElementAlterationFlags(id, flags, &m_alteredCells);
}

void PatchKernel::setCellAlterationFlags(AlterationFlags flags)
{
    CellConstIterator endItr = cellConstEnd();
    for (CellConstIterator itr = cellConstBegin(); itr != endItr; ++itr) {
        setCellAlterationFlags(itr.getId(), flags);
    }
}

void PatchKernel::setCellAlterationFlags(long id, AlterationFlags flags)
{
    setElementAlterationFlags(id, flags, &m_alteredCells);
}

void PatchKernel::unsetCellAlterationFlags(AlterationFlags flags)
{
    unsetElementAlterationFlags(flags, &m_alteredCells);
}

void PatchKernel::unsetCellAlterationFlags(long id, AlterationFlags flags)
{
    unsetElementAlterationFlags(id, flags, &m_alteredCells);
}

bool PatchKernel::testInterfaceAlterationFlags(long id, AlterationFlags flags) const
{
    return testElementAlterationFlags(id, flags, m_alteredInterfaces);
}

PatchKernel::AlterationFlags PatchKernel::getInterfaceAlterationFlags(long id) const
{
    return getElementAlterationFlags(id, m_alteredInterfaces);
}

void PatchKernel::resetInterfaceAlterationFlags(long id, AlterationFlags flags)
{
    resetElementAlterationFlags(id, flags, &m_alteredInterfaces);
}

void PatchKernel::setInterfaceAlterationFlags(AlterationFlags flags)
{
    InterfaceConstIterator endItr = interfaceConstEnd();
    for (InterfaceConstIterator itr = interfaceConstBegin(); itr != endItr; ++itr) {
        setInterfaceAlterationFlags(itr.getId(), flags);
    }
}

void PatchKernel::setInterfaceAlterationFlags(long id, AlterationFlags flags)
{
    setElementAlterationFlags(id, flags, &m_alteredInterfaces);
}

void PatchKernel::unsetInterfaceAlterationFlags(AlterationFlags flags)
{
    unsetElementAlterationFlags(flags, &m_alteredInterfaces);
}

void PatchKernel::unsetInterfaceAlterationFlags(long id, AlterationFlags flags)
{
    unsetElementAlterationFlags(id, flags, &m_alteredInterfaces);
}

bool PatchKernel::testElementAlterationFlags(long id, AlterationFlags flags, const AlterationFlagsStorage &flagsStorage) const
{
    auto storedFlagsItr = flagsStorage.find(id);
    if (storedFlagsItr != flagsStorage.end()) {
        return testAlterationFlags(storedFlagsItr->second, flags);
    } else {
        return testAlterationFlags(FLAG_NONE, flags);
    }
}

bool PatchKernel::testAlterationFlags(AlterationFlags availableFlags, AlterationFlags requestedFlags) const
{
    return ((availableFlags & requestedFlags) == requestedFlags);
}

PatchKernel::AlterationFlags PatchKernel::getElementAlterationFlags(long id, const AlterationFlagsStorage &flagsStorage) const
{
    auto storedFlagsItr = flagsStorage.find(id);
    if (storedFlagsItr != flagsStorage.end()) {
        return storedFlagsItr->second;
    } else {
        return FLAG_NONE;
    }
}

void PatchKernel::resetElementAlterationFlags(long id, AlterationFlags flags, AlterationFlagsStorage *flagsStorage) const
{
    if (flags == FLAG_NONE) {
        flagsStorage->erase(id);
    } else {
        (*flagsStorage)[id] = flags;
    }
}

void PatchKernel::setElementAlterationFlags(long id, AlterationFlags flags, AlterationFlagsStorage *flagsStorage) const
{
    if (flags == FLAG_NONE) {
        return;
    }
 
    auto storedFlagsItr = flagsStorage->find(id);
    if (storedFlagsItr != flagsStorage->end()) {
        storedFlagsItr->second |= flags;
    } else {
        flagsStorage->insert({id, flags});
    }
}

void PatchKernel::unsetElementAlterationFlags(AlterationFlags flags, AlterationFlagsStorage *flagsStorage) const
{
    if (flags == FLAG_NONE) {
        return;
    }
 
    for (auto storedFlagsItr = flagsStorage->begin(); storedFlagsItr != flagsStorage->end();) {
        storedFlagsItr->second &= ~flags;
        if (storedFlagsItr->second == FLAG_NONE) {
            storedFlagsItr = flagsStorage->erase(storedFlagsItr);
        } else {
            ++storedFlagsItr;
        }
    }
}

void PatchKernel::unsetElementAlterationFlags(long id, AlterationFlags flags, AlterationFlagsStorage *flagsStorage) const
{
    if (flags == FLAG_NONE) {
        return;
    }
 
    auto storedFlagsItr = flagsStorage->find(id);
    if (storedFlagsItr == flagsStorage->end()) {
        return;
    }
 
    storedFlagsItr->second &= ~flags;
    if (storedFlagsItr->second == FLAG_NONE) {
        flagsStorage->erase(storedFlagsItr);
    }
}

void PatchKernel::clearBoundingBox()
{
    for (int k = 0; k < 3; ++k) {
        m_boxMinPoint[k]   =   std::numeric_limits<double>::max();
        m_boxMinCounter[k] = 0;
 
        m_boxMaxPoint[k]   = - std::numeric_limits<double>::max();
        m_boxMaxCounter[k] = 0;
    }
 
    setBoundingBoxDirty(getCellCount() > 0);
}

void PatchKernel::setBoundingBox(const std::array<double, 3> &minPoint, const std::array<double, 3> &maxPoint)
{
    m_boxMinPoint = minPoint;
    m_boxMaxPoint = maxPoint;
 
    setBoundingBoxDirty(false);
}

void PatchKernel::getBoundingBox(std::array<double, 3> &minPoint, std::array<double, 3> &maxPoint) const
{
    getBoundingBox(false, minPoint, maxPoint);
}

void PatchKernel::getBoundingBox(bool global, std::array<double, 3> &minPoint, std::array<double, 3> &maxPoint) const
{
    minPoint = m_boxMinPoint;
    maxPoint = m_boxMaxPoint;
 
#if BITPIT_ENABLE_MPI==1
    if (global && isPartitioned()) {
        MPI_Comm communicator = getCommunicator();
        MPI_Allreduce(MPI_IN_PLACE, minPoint.data(), 3, MPI_DOUBLE, MPI_MIN, communicator);
        MPI_Allreduce(MPI_IN_PLACE, maxPoint.data(), 3, MPI_DOUBLE, MPI_MAX, communicator);
    }
#else
    BITPIT_UNUSED(global);
#endif
}

bool PatchKernel::isBoundingBoxFrozen() const
{
    return m_boxFrozen;
}

void PatchKernel::setBoundingBoxFrozen(bool frozen)
{
    m_boxFrozen = frozen;
}

bool PatchKernel::isBoundingBoxDirty(bool global) const
{
    bool isDirty = m_boxDirty;
#if BITPIT_ENABLE_MPI==1
    if (global && isPartitioned()) {
        const auto &communicator = getCommunicator();
        MPI_Allreduce(const_cast<bool *>(&m_boxDirty), &isDirty, 1, MPI_C_BOOL, MPI_LOR, communicator);
    }
#else
    BITPIT_UNUSED(global);
#endif
 
    return isDirty;
}

void PatchKernel::setBoundingBoxDirty(bool dirty)
{
    m_boxDirty = dirty;
}

void PatchKernel::updateBoundingBox(bool forcedUpdated)
{
    if (isBoundingBoxFrozen()) {
        return;
    }
 
    // Check if the bounding box is dirty
    if (!isBoundingBoxDirty() && !forcedUpdated) {
        return;
    }
 
    // Initialize bounding box
    clearBoundingBox();
 
    // Unset the dirty flag in order to be able to update the bounding box
    setBoundingBoxDirty(false);
 
    // Compute bounding box
    for (const auto &vertex : m_vertices) {
        addPointToBoundingBox(vertex.getCoords());
    }
}

void PatchKernel::addPointToBoundingBox(const std::array<double, 3> &point)
{
    if (isBoundingBoxFrozen() || isBoundingBoxDirty()) {
        return;
    }
 
    for (size_t k = 0; k < point.size(); ++k) {
        double value = point[k];
 
        // Update maximum value
        if (utils::DoubleFloatingEqual()(value, m_boxMaxPoint[k], getTol())) {
            ++m_boxMaxCounter[k];
        } else if (value > m_boxMaxPoint[k]) {
            m_boxMaxPoint[k]   = value;
            m_boxMaxCounter[k] = 1;
        }
 
        // Update minimum value
        if (utils::DoubleFloatingEqual()(value, m_boxMinPoint[k], getTol())) {
            ++m_boxMinCounter[k];
        } else if (value < m_boxMinPoint[k]) {
            m_boxMinPoint[k]   = value;
            m_boxMinCounter[k] = 1;
        }
    }
}

void PatchKernel::removePointFromBoundingBox(const std::array<double, 3> &point)
{
    if (isBoundingBoxFrozen() || isBoundingBoxDirty()) {
        return;
    }
 
    double tolerance = getTol();
    for (size_t k = 0; k < point.size(); ++k) {
        double value = point[k];
 
        // Check if maximum value is still valid
        if (utils::DoubleFloatingEqual()(value, m_boxMaxPoint[k], tolerance)) {
            --m_boxMaxCounter[k];
            if (m_boxMaxCounter[k] == 0) {
                setBoundingBoxDirty(true);
                return;
            }
        } else if (value > m_boxMaxPoint[k]) {
            assert(false && "Bounding box is in inconsistent state.");
            setBoundingBoxDirty(true);
            return;
        }
 
        // Update minimum value
        if (utils::DoubleFloatingEqual()(value, m_boxMinPoint[k], tolerance)) {
            --m_boxMinCounter[k];
            if (m_boxMinCounter[k] == 0) {
                setBoundingBoxDirty(true);
                return;
            }
        } else if (value < m_boxMinPoint[k]) {
            assert(false && "Bounding box is in inconsistent state.");
            setBoundingBoxDirty(true);
            return;
        }
    }
}

ConstProxyVector<std::array<double, 3>> PatchKernel::getElementVertexCoordinates(const Element &element) const
{
    // Get the vertices ids
    ConstProxyVector<long> elementVertexIds = element.getVertexIds();
    const int nElementVertices = elementVertexIds.size();
 
    // Build the proxy vector with the coordinates
    ConstProxyVector<std::array<double, 3>> vertexCoordinates(ConstProxyVector<std::array<double, 3>>::INTERNAL_STORAGE, nElementVertices);
    ConstProxyVector<std::array<double, 3>>::storage_pointer storage = vertexCoordinates.storedData();
    getVertexCoords(elementVertexIds.size(), elementVertexIds.data(), storage);
 
    return vertexCoordinates;
}

void PatchKernel::getElementVertexCoordinates(const Element &element, std::unique_ptr<std::array<double, 3>[]> *coordinates) const
{
    ConstProxyVector<long> elementVertexIds = element.getVertexIds();
    getVertexCoords(elementVertexIds.size(), elementVertexIds.data(), coordinates);
}

void PatchKernel::getElementVertexCoordinates(const Element &element, std::array<double, 3> *coordinates) const
{
    ConstProxyVector<long> elementVertexIds = element.getVertexIds();
    getVertexCoords(elementVertexIds.size(), elementVertexIds.data(), coordinates);
}

std::unordered_map<long, std::vector<long>> PatchKernel::binGroupVertices(int nBins)
{
    return PatchKernel::binGroupVertices(m_vertices, nBins);
}

std::unordered_map<long, std::vector<long>> PatchKernel::binGroupVertices(const PiercedVector<Vertex> &vertices, int nBins)
{
    // Update bounding box
    updateBoundingBox();
 
    // Bin's spacing
    double dx = std::max(1.0e-12, m_boxMaxPoint[0] - m_boxMinPoint[0]) / ((double) nBins);
    double dy = std::max(1.0e-12, m_boxMaxPoint[1] - m_boxMinPoint[1]) / ((double) nBins);
    double dz = std::max(1.0e-12, m_boxMaxPoint[2] - m_boxMinPoint[2]) / ((double) nBins);
 
    // Identify bins of vertices
    std::unordered_map<long, std::vector<long>> bins;
    for (const Vertex &vertex : vertices) {
        const std::array<double, 3> &coordinates = vertex.getCoords();
 
        int i = std::min(nBins - 1L, (long) ((coordinates[0] - m_boxMinPoint[0]) / dx));
        int j = std::min(nBins - 1L, (long) ((coordinates[1] - m_boxMinPoint[1]) / dy));
        int k = std::min(nBins - 1L, (long) ((coordinates[2] - m_boxMinPoint[2]) / dz));
 
        long binId = nBins * nBins * k + nBins * j + i;
        bins[binId].emplace_back(vertex.getId());
    }
 
    return bins;
}

void PatchKernel::translate(const std::array<double, 3> &translation)
{
    // Translate the patch
    for (auto &vertex : m_vertices) {
        vertex.translate(translation);
    }
 
    // Update the bounding box
    if (!isBoundingBoxFrozen() || isBoundingBoxDirty()) {
        m_boxMinPoint += translation;
        m_boxMaxPoint += translation;
    }
}

void PatchKernel::translate(double sx, double sy, double sz)
{
    translate({{sx, sy, sz}});
}

void PatchKernel::rotate(const std::array<double, 3> &n0, const std::array<double, 3> &n1, double angle)
{
    // Translate the patch
    for (auto &vertex : m_vertices) {
        vertex.rotate(n0, n1, angle);
    }
 
    // Update the bounding box
    if (!isBoundingBoxFrozen() && !isBoundingBoxDirty()) {
        // Save current bounding box
        std::array<std::array<double, 3>, 3> originalBox = {{m_boxMinPoint, m_boxMaxPoint}};
 
        // Clear bounding box
        clearBoundingBox();
 
        // Unset the dirty flag in order to be able to update the bounding box
        setBoundingBoxDirty(false);
 
        // Evaluate rotated bounding box
        for (int i = 0; i < 2; ++i) {
            double xCorner = originalBox[i][0];
            for (int j = 0; j < 2; ++j) {
                double yCorner = originalBox[j][1];
                for (int k = 0; k < 2; ++k) {
                    double zCorner = originalBox[k][2];
 
                    std::array<double, 3> corner = {{xCorner, yCorner, zCorner}};
                    corner = CGElem::rotatePoint(corner, n0, n1, angle);
                    addPointToBoundingBox(corner);
                }
            }
        }
    }
}

void PatchKernel::rotate(double n0x, double n0y, double n0z, double n1x,
                         double n1y, double n1z, double angle)
{
    rotate({{n0x, n0y, n0z}}, {{n1x, n1y, n1z}}, angle);
}

void PatchKernel::scale(const std::array<double, 3> &scaling)
{
    scale(scaling, m_boxMinPoint);
}

void PatchKernel::scale(const std::array<double, 3> &scaling, const std::array<double, 3> &center)
{
    // Scale the patch
    for (auto &vertex : m_vertices) {
        vertex.scale(scaling, center);
    }
 
    // Update the bounding box
    if (!isBoundingBoxFrozen() || isBoundingBoxDirty()) {
        for (int k = 0; k < 3; ++k) {
            m_boxMinPoint[k] = center[k] + scaling[k] * (m_boxMinPoint[k] - center[k]);
            m_boxMaxPoint[k] = center[k] + scaling[k] * (m_boxMaxPoint[k] - center[k]);
        }
    }
}

void PatchKernel::scale(double scaling)
{
    scale({{scaling, scaling, scaling}}, m_boxMinPoint);
}

void PatchKernel::scale(double scaling, const std::array<double, 3> &center)
{
    scale({{scaling, scaling, scaling}}, center);
}

void PatchKernel::scale(double sx, double sy, double sz)
{
    scale({{sx, sy, sz}}, m_boxMinPoint);
}

void PatchKernel::scale(double sx, double sy, double sz, const std::array<double, 3> &center)
{
    scale({{sx, sy, sz}}, center);
}

void PatchKernel::setTol(double tolerance)
{
    _setTol(tolerance);
 
    m_toleranceCustom = true;
}

void PatchKernel::_setTol(double tolerance)
{
    m_tolerance = tolerance;
}

double PatchKernel::getTol() const
{
    return m_tolerance;
}

void PatchKernel::resetTol()
{
    _resetTol();
 
    m_toleranceCustom = false;
}

void PatchKernel::_resetTol()
{
    const double DEFAULT_TOLERANCE = 1e-14;
 
    m_tolerance = DEFAULT_TOLERANCE;
}

bool PatchKernel::isTolCustomized() const
{
    return m_toleranceCustom;
}

void PatchKernel::extractEnvelope(PatchKernel &envelope) const
{
 
    // ====================================================================== //
    // RESIZE DATA STRUCTURES                                                 //
    // ====================================================================== //
    envelope.reserveVertices(envelope.getVertexCount() + countBorderVertices());
    envelope.reserveCells(envelope.getCellCount() + countBorderFaces());
 
    // ====================================================================== //
    // LOOP OVER CELLS                                                        //
    // ====================================================================== //
    std::unordered_map<long, long> vertexMap;
    for (const Cell &cell : m_cells) {
        int nCellFaces = cell.getFaceCount();
        for (int i = 0; i < nCellFaces; ++i) {
            if (!cell.isFaceBorder(i)) {
                continue;
            }
 
            // Add face vertices to the envelope and get face
            // connectivity in the envelope
            ConstProxyVector<long> faceConnect = cell.getFaceConnect(i);
            int nFaceVertices = faceConnect.size();
 
            std::unique_ptr<long[]> faceEnvelopeConnect = std::unique_ptr<long[]>(new long[nFaceVertices]);
            for (int j = 0; j < nFaceVertices; ++j) {
                long vertexId = faceConnect[j];
 
                // If the vertex is not yet in the envelope
                // add it.
                if (vertexMap.count(vertexId) == 0) {
                    const Vertex &vertex = getVertex(vertexId);
                    VertexIterator envelopeVertex = envelope.addVertex(vertex);
                    vertexMap[vertexId] = envelopeVertex->getId();
                }
 
                // Update face ace connectivity in the envelope
                faceEnvelopeConnect[j] = vertexMap.at(vertexId);
            }
 
            // Add face to envelope
            ElementType faceType = cell.getFaceType(i);
            envelope.addCell(faceType, std::move(faceEnvelopeConnect));
        }
    }
}

void PatchKernel::displayTopologyStats(std::ostream &out, unsigned int padding) const
{
    std::string indent = std::string(padding, ' ');
 
    // ====================================================================== //
    // VERTEX STATS                                                           //
    // ====================================================================== //
    out << indent<< "Vertices --------------------------------"     << std::endl;
    out << indent<< "  # vertices        " << getVertexCount()      << std::endl;
    out << indent<< "  # orphan vertices " << countOrphanVertices() << std::endl;
    out << indent<< "  # border vertices " << countBorderVertices()   << std::endl;
        //out << indent<< "  # free vertices   " << countDoubleVertices()   << std::endl;
 
    // ====================================================================== //
    // FACE STATS                                                             //
    // ====================================================================== //
    out << indent<< "Faces -----------------------------------"     << std::endl;
    out << indent<< "  # faces           " << countFaces()          << std::endl;
    out << indent<< "  # border faces    " << countBorderFaces()    << std::endl;
 
    // ====================================================================== //
    // CELLS STATS                                                            //
    // ====================================================================== //
    out << indent<< "Cells -----------------------------------"     << std::endl;
    out << indent<< "  # cells           " << getCellCount()        << std::endl;
    out << indent<< "  # orphan cells    " << countOrphanCells()    << std::endl;
    out << indent<< "  # border cells    " << countBorderCells()    << std::endl;
        //out << indent<< "  # free vertices   " << countDoubleCells()   << std::endl;
}

void PatchKernel::displayVertices(std::ostream &out, unsigned int padding) const
{
    std::string indent = std::string(padding, ' ');
    for (const Vertex &vertex : m_vertices) {
        out << indent << "vertex: " << std::endl;
        vertex.display(out, padding + 2);
    }
}

void PatchKernel::displayCells(std::ostream &out, unsigned int padding) const
{
    std::string indent = std::string(padding, ' ');
    for (const Cell &cell : m_cells) {
        out << indent << "cell: " << std::endl;
        cell.display(out, padding + 2);
    }
}

void PatchKernel::displayInterfaces(std::ostream &out, unsigned int padding) const
{
    std::string indent = std::string(padding, ' ');
    for (const Interface &interface : m_interfaces) {
        out << indent << "interface: " << std::endl;
        interface.display(out, padding + 2);
    }
}

VTKUnstructuredGrid & PatchKernel::getVTK()
{
    return m_vtk;
}

const VTKUnstructuredGrid & PatchKernel::getVTK() const
{
    return m_vtk;
}

PatchKernel::WriteTarget PatchKernel::getVTKWriteTarget() const
{
    return m_vtkWriteTarget;
}

void PatchKernel::setVTKWriteTarget(WriteTarget writeTarget)
{
    m_vtkWriteTarget = writeTarget;
}

const PatchKernel::CellConstRange PatchKernel::getVTKCellWriteRange() const
{
    if (m_vtkWriteTarget == WRITE_TARGET_CELLS_ALL) {
        return CellConstRange(cellConstBegin(), cellConstEnd());
#if BITPIT_ENABLE_MPI==1
    } else if (m_vtkWriteTarget == WRITE_TARGET_CELLS_INTERNAL) {
        return CellConstRange(internalCellConstBegin(), internalCellConstEnd());
#endif
    } else {
        return CellConstRange(cellConstEnd(), cellConstEnd());
    }
}

void PatchKernel::replaceVTKStreamer(const VTKBaseStreamer *original, VTKBaseStreamer *updated)
{
    // Update the VTK streamer
    //
    // The pointer to VTK streamers are copied, if there are pointer to the
    // original object they have to be replace with a pointer to this object.
    // Geometry fields
    std::vector<std::string> streamedGeomFields;
    streamedGeomFields.reserve(m_vtk.getGeomDataCount());
    for (auto itr = m_vtk.getGeomDataBegin(); itr != m_vtk.getGeomDataEnd(); ++itr) {
        const VTKField &field = *itr;
        if (&field.getStreamer() != original) {
            continue;
        }
 
        streamedGeomFields.push_back(field.getName());
    }
 
    for (const std::string &name : streamedGeomFields) {
        const VTKField &field = *(m_vtk.findGeomData(name));
        VTKField updatedField(field);
        updatedField.setStreamer(*updated);
 
        m_vtk.setGeomData(std::move(updatedField));
    }
 
    std::vector<std::string> streamedDataFields;
    streamedDataFields.reserve(m_vtk.getDataCount());
    for (auto itr = m_vtk.getDataBegin(); itr != m_vtk.getDataEnd(); ++itr) {
        const VTKField &field = *itr;
        if (&field.getStreamer() != original) {
            continue;
        }
 
        streamedDataFields.push_back(field.getName());
    }
 
    for (const std::string &name : streamedDataFields) {
        const VTKField &field = *(m_vtk.findData(name));
        VTKField updatedField(field);
        updatedField.setStreamer(*updated);
 
        m_vtk.removeData(field.getName());
        m_vtk.addData(std::move(updatedField));
    }
}

void PatchKernel::flushData(std::fstream &stream, const std::string &name, VTKFormat format)
{
    if (name == "Points") {
        VertexConstIterator endItr = vertexConstEnd();
        for (VertexConstIterator itr = vertexConstBegin(); itr != endItr; ++itr) {
            std::size_t vertexRawId = itr.getRawIndex();
            long vertexVTKId = m_vtkVertexMap.rawAt(vertexRawId);
            if (vertexVTKId != Vertex::NULL_ID) {
                const Vertex &vertex = m_vertices.rawAt(vertexRawId);
                flushValue(stream, format,vertex.getCoords());
            }
        }
    } else if (name == "offsets") {
        long offset = 0;
        for (const Cell &cell : getVTKCellWriteRange()) {
            offset += cell.getVertexCount();
            flushValue(stream, format,offset);
        }
    } else if (name == "types") {
        for (const Cell &cell : getVTKCellWriteRange()) {
            VTKElementType VTKType;
            switch (cell.getType())  {
 
            case ElementType::VERTEX:
                VTKType = VTKElementType::VERTEX;
                break;
 
            case ElementType::LINE:
                VTKType = VTKElementType::LINE;
                break;
 
            case ElementType::TRIANGLE:
                VTKType = VTKElementType::TRIANGLE;
                break;
 
            case ElementType::PIXEL:
                VTKType = VTKElementType::PIXEL;
                break;
 
            case ElementType::QUAD:
                VTKType = VTKElementType::QUAD;
                break;
 
            case ElementType::POLYGON:
                VTKType = VTKElementType::POLYGON;
                break;
 
            case ElementType::TETRA:
                VTKType = VTKElementType::TETRA;
                break;
 
            case ElementType::VOXEL:
                VTKType = VTKElementType::VOXEL;
                break;
 
            case ElementType::HEXAHEDRON:
                VTKType = VTKElementType::HEXAHEDRON;
                break;
 
            case ElementType::WEDGE:
                VTKType = VTKElementType::WEDGE;
                break;
 
            case ElementType::PYRAMID:
                VTKType = VTKElementType::PYRAMID;
                break;
 
            case ElementType::POLYHEDRON:
                VTKType = VTKElementType::POLYHEDRON;
                break;
 
            default:
                VTKType = VTKElementType::UNDEFINED;
                break;
 
            }
 
            flushValue(stream, format,(int) VTKType);
        }
    } else if (name == "connectivity") {
        for (const Cell &cell : getVTKCellWriteRange()) {
            ConstProxyVector<long> cellVertexIds = cell.getVertexIds();
            const int nCellVertices = cellVertexIds.size();
            for (int k = 0; k < nCellVertices; ++k) {
                long vertexId = cellVertexIds[k];
                long vtkVertexId = m_vtkVertexMap.at(vertexId);
                flushValue(stream, format,vtkVertexId);
            }
        }
    } else if (name == "faces") {
        for (const Cell &cell : getVTKCellWriteRange()) {
            if (cell.getDimension() <= 2 || cell.hasInfo()) {
                flushValue(stream, format,(long) 0);
            } else {
                std::vector<long> faceStream = cell.getFaceStream();
                Cell::renumberFaceStream(m_vtkVertexMap, &faceStream);
                int faceStreamSize = faceStream.size();
                for (int k = 0; k < faceStreamSize; ++k) {
                    flushValue(stream, format,faceStream[k]);
                }
            }
        }
    } else if (name == "faceoffsets") {
        long offset = 0;
        for (const Cell &cell : getVTKCellWriteRange()) {
            if (cell.getDimension() <= 2 || cell.hasInfo()) {
                offset += 1;
            } else {
                offset += cell.getFaceStreamSize();
            }
 
            flushValue(stream, format,offset);
        }
    } else if (name == "cellIndex") {
        for (const Cell &cell : getVTKCellWriteRange()) {
            flushValue(stream, format,cell.getId());
        }
    } else if (name == "PID") {
        for (const Cell &cell : getVTKCellWriteRange()) {
            flushValue(stream, format,cell.getPID());
        }
    } else if (name == "vertexIndex") {
        VertexConstIterator endItr = vertexConstEnd();
        for (VertexConstIterator itr = vertexConstBegin(); itr != endItr; ++itr) {
            std::size_t vertexRawId = itr.getRawIndex();
            long vertexVTKId = m_vtkVertexMap.rawAt(vertexRawId);
            if (vertexVTKId != Vertex::NULL_ID) {
                long vertexId = itr.getId();
                flushValue(stream, format,vertexId);
            }
        }
#if BITPIT_ENABLE_MPI==1
    } else if (name == "cellGlobalIndex") {
        PatchNumberingInfo numberingInfo(this);
        for (const Cell &cell : getVTKCellWriteRange()) {
            flushValue(stream, format,numberingInfo.getCellGlobalId(cell.getId()));
        }
    } else if (name == "cellRank") {
        for (const Cell &cell : getVTKCellWriteRange()) {
            flushValue(stream, format,getCellOwner(cell.getId()));
        }
    } else if (name == "cellHaloLayer") {
        for (const Cell &cell : getVTKCellWriteRange()) {
            flushValue(stream, format,getCellHaloLayer(cell.getId()));
        }
    } else if (name == "vertexRank") {
        VertexConstIterator endItr = vertexConstEnd();
        for (VertexConstIterator itr = vertexConstBegin(); itr != endItr; ++itr) {
            std::size_t vertexRawId = itr.getRawIndex();
            long vertexVTKId = m_vtkVertexMap.rawAt(vertexRawId);
            if (vertexVTKId != Vertex::NULL_ID) {
                flushValue(stream, format,getVertexOwner(itr.getId()));
            }
        }
#endif
    }
}

void PatchKernel::consecutiveRenumberVertices(long offset)
{
    // Renumber vertices
    std::unordered_map<long, long > map = consecutiveItemRenumbering(m_vertices, offset);
    
    // Renumber cell connectivity
    for(Cell &cell : m_cells) {
        cell.renumberVertices(map);
    }
 
    // Renumber interface connectivity
    for(Interface &interface : getInterfaces()) {
        interface.renumberVertices(map);
    }
 
    // Reset index generator
    if (isVertexAutoIndexingEnabled()) {
        createVertexIndexGenerator(true);
    }
}   

void PatchKernel::consecutiveRenumberCells(long offset)
{
    // Renumber cells
    std::unordered_map<long, long > map = consecutiveItemRenumbering(m_cells, offset);
    
    // Renumber cell adjacencies
    for (auto &cell: m_cells) {
        long *adjacencies = cell.getAdjacencies();
        int nCellAdjacencies = cell.getAdjacencyCount();
        for (int i = 0; i < nCellAdjacencies; ++i) {
            long &neighId = adjacencies[i];
            neighId = map[neighId];
        }
    }
    
    // Renumber interface owner and neighbour
    for (Interface &interface: m_interfaces) {
        long ownerId = interface.getOwner();
        int ownerFace = interface.getOwnerFace();
        interface.setOwner(map.at(ownerId), ownerFace);
 
        long neighId = interface.getNeigh();
        if (neighId >= 0) {
            int neighFace = interface.getNeighFace();
            interface.setNeigh(map.at(neighId), neighFace);
        }
    }
 
    // Renumber last internal and first ghost markers
    if (m_lastInternalCellId >= 0) {
        m_lastInternalCellId = map.at(m_lastInternalCellId);
    }
 
#if BITPIT_ENABLE_MPI==1
    if (m_firstGhostCellId >= 0) {
        m_firstGhostCellId = map.at(m_firstGhostCellId);
    }
#endif
 
    // Reset index generator
    if (isCellAutoIndexingEnabled()) {
        createCellIndexGenerator(true);
    }
 
#if BITPIT_ENABLE_MPI==1
    // Update partitioning information
    if (isPartitioned()) {
        updatePartitioningInfo(true);
    }
#endif
}   

void PatchKernel::consecutiveRenumberInterfaces(long offset)
{
    // Renumber interfaces
    std::unordered_map<long, long > map = consecutiveItemRenumbering(m_interfaces, offset);
    
    // Renumber cell interfaces
    for (Cell &cell: m_cells) {
        long *interfaces = cell.getInterfaces();
        int nCellInterfaces = cell.getInterfaceCount();
        for (int i = 0; i < nCellInterfaces; ++i) {
            long &interfaceId = interfaces[i];
            interfaceId = map.at(interfaceId);
        }
    }
 
    // Reset index generator
    if (isInterfaceAutoIndexingEnabled()) {
        createInterfaceIndexGenerator(true);
    }
}

void PatchKernel::consecutiveRenumber(long vertexOffset, long cellOffset, long interfaceOffset)
{
    consecutiveRenumberVertices(vertexOffset);
    consecutiveRenumberCells(cellOffset);
    consecutiveRenumberInterfaces(interfaceOffset);
}

int PatchKernel::getDumpVersion() const
{
    const int KERNEL_DUMP_VERSION = 12;
 
    return (KERNEL_DUMP_VERSION + _getDumpVersion());
}

bool PatchKernel::dump(std::ostream &stream)
{
    // Update the patch
    update();
 
    // Dump the patch
    const PatchKernel *constPatch = this;
    return constPatch->dump(stream);
}

bool PatchKernel::dump(std::ostream &stream) const
{
    // Dumping a dirty patch is not supported.
    bool dirty = isDirty(true);
    assert(!dirty && "Dumping a patch that is not up-to-date is not supported.");
    if (dirty) {
        return false;
    }
 
    // Version
    utils::binary::write(stream, getDumpVersion());
 
    // Generic information
    utils::binary::write(stream, m_id);
    utils::binary::write(stream, m_dimension);
    utils::binary::write(stream, m_vtk.getName());
#if BITPIT_ENABLE_MPI==1
    utils::binary::write(stream, m_haloSize);
#else
    utils::binary::write(stream, 0);
#endif
 
    // Adaption information
    utils::binary::write(stream, m_adaptionMode);
    utils::binary::write(stream, m_adaptionStatus);
 
    // Partition information
#if BITPIT_ENABLE_MPI==1
    utils::binary::write(stream, m_partitioningMode);
    utils::binary::write(stream, m_partitioningStatus);
#else
    utils::binary::write(stream, PARTITIONING_DISABLED);
    utils::binary::write(stream, PARTITIONING_CLEAN);
#endif
 
    // Adjacencies build strategy
    utils::binary::write(stream, m_adjacenciesBuildStrategy);
 
    // Dump interfaces build strategy
    utils::binary::write(stream, m_interfacesBuildStrategy);
 
    // VTK data
    utils::binary::write(stream, m_vtkWriteTarget);
 
    // Specific dump
    _dump(stream);
 
    // Geometric tolerance
    utils::binary::write(stream, (int) m_toleranceCustom);
    if (m_toleranceCustom) {
        utils::binary::write(stream, m_tolerance);
    }
 
    // Index generators
    bool hasVertexAutoIndexing = isVertexAutoIndexingEnabled();
    utils::binary::write(stream, hasVertexAutoIndexing);
    if (hasVertexAutoIndexing) {
        dumpVertexAutoIndexing(stream);
    }
 
    bool hasInterfaceAutoIndexing = isInterfaceAutoIndexingEnabled();
    utils::binary::write(stream, hasInterfaceAutoIndexing);
    if (hasInterfaceAutoIndexing) {
        dumpInterfaceAutoIndexing(stream);
    }
 
    bool hasCellAutoIndexing = isCellAutoIndexingEnabled();
    utils::binary::write(stream, hasCellAutoIndexing);
    if (hasCellAutoIndexing) {
        dumpCellAutoIndexing(stream);
    }
 
    // The patch has been dumped successfully
    return true;
}

void PatchKernel::restore(std::istream &stream, bool reregister)
{
    // Reset the patch
    reset();
 
    // Version
    int version;
    utils::binary::read(stream, version);
    if (version != getDumpVersion()) {
        throw std::runtime_error ("The version of the file does not match the required version");
    }
 
    // Id
    int id;
    utils::binary::read(stream, id);
    if (reregister) {
        setId(id);
    }
 
    // Dimension
    int dimension;
    utils::binary::read(stream, dimension);
    setDimension(dimension);
 
    // Name
    std::string name;
    utils::binary::read(stream, name);
    m_vtk.setName(name);
 
    // Halo size
#if BITPIT_ENABLE_MPI==1
    utils::binary::read(stream, m_haloSize);
#else
    int dummyHaloSize;
    utils::binary::read(stream, dummyHaloSize);
#endif
 
    // Adaption information
    utils::binary::read(stream, m_adaptionMode);
    utils::binary::read(stream, m_adaptionStatus);
 
    // Partition information
#if BITPIT_ENABLE_MPI==1
    utils::binary::read(stream, m_partitioningMode);
    utils::binary::read(stream, m_partitioningStatus);
#else
    PartitioningStatus dummyPartitioningMode;
    utils::binary::read(stream, dummyPartitioningMode);
 
    PartitioningStatus dummyPartitioningStatus;
    utils::binary::read(stream, dummyPartitioningStatus);
#endif
 
    // Adjacencies build strategy
    utils::binary::read(stream, m_adjacenciesBuildStrategy);
 
    // Interfaces build strategy
    utils::binary::read(stream, m_interfacesBuildStrategy);
 
    // VTK data
    utils::binary::read(stream, m_vtkWriteTarget);
 
    // Specific restore
    _restore(stream);
 
#if BITPIT_ENABLE_MPI==1
    // Update partitioning information
    if (isPartitioned()) {
        updatePartitioningInfo(true);
    }
#endif
 
    // Geometric tolerance
    int hasCustomTolerance;
    utils::binary::read(stream, hasCustomTolerance);
    if (hasCustomTolerance) {
        double tolerance;
        utils::binary::read(stream, tolerance);
        setTol(tolerance);
    } else {
        resetTol();
    }
 
    // Index generators
    bool hasVertexAutoIndexing;
    utils::binary::read(stream, hasVertexAutoIndexing);
    if (hasVertexAutoIndexing) {
        restoreVertexAutoIndexing(stream);
    } else {
        setVertexAutoIndexing(false);
    }
 
    bool hasInterfaceAutoIndexing;
    utils::binary::read(stream, hasInterfaceAutoIndexing);
    if (hasInterfaceAutoIndexing) {
        restoreInterfaceAutoIndexing(stream);
    } else {
        setInterfaceAutoIndexing(false);
    }
 
    bool hasCellAutoIndexing;
    utils::binary::read(stream, hasCellAutoIndexing);
    if (hasCellAutoIndexing) {
        restoreCellAutoIndexing(stream);
    } else {
        setCellAutoIndexing(false);
    }
}

void PatchKernel::mergeAdaptionInfo(std::vector<adaption::Info> &&source, std::vector<adaption::Info> &destination)
{
    if (source.empty()) {
        return;
    } else if (destination.empty()) {
        destination.swap(source);
        return;
    }
 
    throw std::runtime_error ("Unable to merge the adaption info.");
}
 
}