Mesh.cpp

Go to the documentation of this file.

#include "Mesh.hpp"
#include "Mesh_PartitionHelpers.hpp"
#include "Solver/Direct.hpp"
 
#include <cstdlib>
#include <string>
#include <map>
#include <set>
#include <unordered_map>
#include <unordered_set>
#include <filesystem>
#include <fmt/core.h>
#include "DNDS/EigenPCH.hpp"
#include "Geom/Mesh_DeviceView.hpp"
#include "SerialAdjReordering.hpp"
 
#include "DNDS/DeviceStorage.hxx"
 
namespace DNDS
{
    // DNDS_DEVICE_STORAGE_BASE_DELETER_INST(Geom::ElemInfo, )
    // DNDS_DEVICE_STORAGE_INST(Geom::ElemInfo, DeviceBackend::Host, )
}
 
namespace DNDS::Geom
{
 
    /**
     * \brief Reserved skeleton for parallel topology interpolation.
     *
     * This commented-out method is preserved as a placeholder for a future generic
     * topology management API. Currently, mesh topology (cell2face, face2cell, etc.)
     * is constructed serially and distributed. However, a parallel interpolation
     * capability may be needed for:
     *   - Dynamic mesh adaptation (refinement/coarsening)
     *   - Load balancing with topology migration
     *   - Multi-physics coupling with different mesh partitions
     *
     * If implementing parallel topology construction, this skeleton provides the
     * basic structure for building face-based adjacencies from cell2node without
     * requiring a serial global mesh on rank 0.
     *
     * \todo Evaluate need for parallel topology API in Phase 5 refactoring.
     */
    // void UnstructuredMeshSerialRW::InterpolateTopology()
    // {
    //     // count node 2 face
    //     DNDS_MAKE_SSP(cell2faceSerial, mesh->getMPI());
    //     DNDS_MAKE_SSP(face2cellSerial, mesh->getMPI());
    //     DNDS_MAKE_SSP(face2nodeSerial, mesh->getMPI());
    //     DNDS_MAKE_SSP(faceElemInfoSerial, ElemInfo::CommType(), ElemInfo::CommMult(), mesh->getMPI());
 
    //     if (mRank != mesh->getMPI().rank)
    //         return;
 
    //     for (DNDS::index iCell = 0; iCell <cell2nodeSerial->Size(); iCell++)
    //     {
    //         // Parallel face construction logic would go here
    //     }
    // }
 
    /*******************************************************************************************************************/
    /*******************************************************************************************************************/
    /*******************************************************************************************************************/
 
    // Partition2LocalIdx, Partition2Serial2Global, ConvertAdjSerial2Global,
    // TransferDataSerial2Global are now in Mesh_PartitionHelpers.hpp
    // so they can be shared with Mesh_ReadSerializeDistributed.cpp.
 
    //! inefficient, use Partition2Serial2Global ! only used for convenient comparison
    void PushInfo2Serial2Global(std::vector<DNDS::index> &serial2Global,
                                DNDS::index localSize,
                                const std::vector<DNDS::index> &pushIndex,
                                const std::vector<DNDS::index> &pushIndexStart,
                                const DNDS::MPIInfo &mpi)
    {
        tIndPair Serial2Global;
        Serial2Global.InitPair("Serial2Global", mpi);
        Serial2Global.father->Resize(localSize);
        Serial2Global.TransAttach();
        Serial2Global.trans.createFatherGlobalMapping();
        Serial2Global.trans.createGhostMapping(pushIndex, pushIndexStart);
        Serial2Global.trans.createMPITypes();
        Serial2Global.son->createGlobalMapping();
        // Set son to son's global
        for (DNDS::index iSon = 0; iSon < Serial2Global.son->Size(); iSon++)
            (*Serial2Global.son)[iSon] = Serial2Global.son->pLGlobalMapping->operator()(mpi.rank, iSon);
        Serial2Global.trans.pushOnce();
        serial2Global.resize(localSize);
        for (DNDS::index iFat = 0; iFat < Serial2Global.father->Size(); iFat++)
            serial2Global[iFat] = Serial2Global.father->operator[](iFat);
    }
 
    // template <class TAdj = tAdj1>
    // void ConvertAdjSerial2Global(TAdj &arraySerialAdj,
    //                              const std::vector<DNDS::index> &partitionJSerial2Global,
    //                              const DNDS::MPIInfo &mpi)
    // {
    // }
    /*******************************************************************************************************************/
    /*******************************************************************************************************************/
    /*******************************************************************************************************************/
 
    void UnstructuredMeshSerialRW::
        PartitionReorderToMeshCell2Cell()
    {
        if (mesh->getMPI().rank == mRank)
            DNDS::log() << "UnstructuredMeshSerialRW === Doing  PartitionReorderToMeshCell2Cell" << std::endl;
        DNDS_assert(cnPart == mesh->getMPI().size);
        // * 1: get the nodal partition
        nodePartition.resize(coordSerial->Size(), static_cast<DNDS::MPI_int>(INT32_MAX));
        for (DNDS::index iCell = 0; iCell < cell2nodeSerial->Size(); iCell++)
            for (DNDS::rowsize ic2n = 0; ic2n < (*cell2nodeSerial).RowSize(iCell); ic2n++)
                nodePartition[(*cell2nodeSerial)(iCell, ic2n)] = std::min(nodePartition[(*cell2nodeSerial)(iCell, ic2n)], cellPartition.at(iCell));
        // * 1: get the bnd partition
        bndPartition.resize(bnd2cellSerial->Size());
        for (DNDS::index iBnd = 0; iBnd < bnd2cellSerial->Size(); iBnd++)
            bndPartition[iBnd] = cellPartition[(*bnd2cellSerial)(iBnd, 0)];
 
        std::vector<DNDS::index> cell_push, cell_pushStart, node_push, node_pushStart, bnd_push, bnd_pushStart;
        Partition2LocalIdx(cellPartition, cell_push, cell_pushStart, mesh->getMPI());
        Partition2LocalIdx(nodePartition, node_push, node_pushStart, mesh->getMPI());
        Partition2LocalIdx(bndPartition, bnd_push, bnd_pushStart, mesh->getMPI());
        std::vector<DNDS::index> cell_Serial2Global, node_Serial2Global, bnd_Serial2Global;
        Partition2Serial2Global(cellPartition, cell_Serial2Global, mesh->getMPI(), mesh->getMPI().size);
        Partition2Serial2Global(nodePartition, node_Serial2Global, mesh->getMPI(), mesh->getMPI().size);
        // Partition2Serial2Global(bndPartition, bnd_Serial2Global, mesh->getMPI(), mesh->getMPI().size);//seems not needed for now
        // PushInfo2Serial2Global(cell_Serial2Global, cellPartition.size(), cell_push, cell_pushStart, mesh->getMPI());//*safe validation version
        // PushInfo2Serial2Global(node_Serial2Global, nodePartition.size(), node_push, node_pushStart, mesh->getMPI());//*safe validation version
        // PushInfo2Serial2Global(bnd_Serial2Global, bndPartition.size(), bnd_push, bnd_pushStart, mesh->getMPI());    //*safe validation version
        if (mesh->getMPI().rank == mRank)
            DNDS::log() << "UnstructuredMeshSerialRW === Doing PartitionReorderToMeshCell2Cell ConvertAdjSerial2Global" << std::endl;
        ConvertAdjSerial2Global(cell2nodeSerial, node_Serial2Global, mesh->getMPI());
        // !cell2cell discarded
        // ConvertAdjSerial2Global(cell2cellSerial, cell_Serial2Global, mesh->getMPI());
        ConvertAdjSerial2Global(bnd2nodeSerial, node_Serial2Global, mesh->getMPI());
        ConvertAdjSerial2Global(bnd2cellSerial, cell_Serial2Global, mesh->getMPI());
 
        mesh->coords.InitPair("coords", mesh->getMPI());
        mesh->cellElemInfo.InitPair("cellElemInfo", ElemInfo::CommType(), ElemInfo::CommMult(), mesh->getMPI());
        mesh->bndElemInfo.InitPair("bndElemInfo", ElemInfo::CommType(), ElemInfo::CommMult(), mesh->getMPI());
        mesh->cell2node.InitPair("cell2node", mesh->getMPI());
        mesh->cell2cellOrig.InitPair("cell2cellOrig", mesh->getMPI());
        mesh->node2nodeOrig.InitPair("node2nodeOrig", mesh->getMPI());
        mesh->bnd2bndOrig.InitPair("bnd2bndOrig", mesh->getMPI());
        if (mesh->isPeriodic)
            mesh->cell2nodePbi.InitPair("cell2nodePbi", mesh->getMPI());
        // !cell2cell discarded
        mesh->bnd2node.InitPair("bnd2node", mesh->getMPI());
        mesh->bnd2cell.InitPair("bnd2cell", mesh->getMPI());
        if (mesh->isPeriodic)
            mesh->bnd2nodePbi.InitPair("bnd2nodePbi", mesh->getMPI());
 
        // coord transferring
        if (mesh->getMPI().rank == mRank)
            DNDS::log() << "UnstructuredMeshSerialRW === Doing PartitionReorderToMeshCell2Cell Trasfer Data Coord" << std::endl;
        TransferDataSerial2Global(coordSerial, mesh->coords.father, node_push, node_pushStart, mesh->getMPI());
        TransferDataSerial2Global(node2nodeOrigSerial, mesh->node2nodeOrig.father, node_push, node_pushStart, mesh->getMPI());
 
        if (mesh->getMPI().rank == mRank)
            DNDS::log() << "UnstructuredMeshSerialRW === Doing PartitionReorderToMeshCell2Cell Trasfer Data Cell" << std::endl;
        // cells transferring
        // !cell2cell discarded
        // TransferDataSerial2Global(cell2cellSerial, mesh->cell2cell.father, cell_push, cell_pushStart, mesh->getMPI());
        TransferDataSerial2Global(cell2nodeSerial, mesh->cell2node.father, cell_push, cell_pushStart, mesh->getMPI());
        TransferDataSerial2Global(cell2cellOrigSerial, mesh->cell2cellOrig.father, cell_push, cell_pushStart, mesh->getMPI());
        if (mesh->isPeriodic)
            TransferDataSerial2Global(cell2nodePbiSerial, mesh->cell2nodePbi.father, cell_push, cell_pushStart, mesh->getMPI());
        TransferDataSerial2Global(cellElemInfoSerial, mesh->cellElemInfo.father, cell_push, cell_pushStart, mesh->getMPI());
        if (mesh->getMPI().rank == mRank)
            DNDS::log() << "UnstructuredMeshSerialRW === Doing PartitionReorderToMeshCell2Cell Trasfer Data Bnd" << std::endl;
        // bnds transferring
        TransferDataSerial2Global(bnd2cellSerial, mesh->bnd2cell.father, bnd_push, bnd_pushStart, mesh->getMPI());
        TransferDataSerial2Global(bnd2nodeSerial, mesh->bnd2node.father, bnd_push, bnd_pushStart, mesh->getMPI());
        TransferDataSerial2Global(bndElemInfoSerial, mesh->bndElemInfo.father, bnd_push, bnd_pushStart, mesh->getMPI());
        TransferDataSerial2Global(bnd2bndOrigSerial, mesh->bnd2bndOrig.father, bnd_push, bnd_pushStart, mesh->getMPI());
        if (mesh->isPeriodic)
            TransferDataSerial2Global(bnd2nodePbiSerial, mesh->bnd2nodePbi.father, bnd_push, bnd_pushStart, mesh->getMPI());
 
        {
            DNDS::MPISerialDo(mesh->getMPI(), [&]()
                              { log() << "[" << mesh->getMPI().rank << ": nCell " << mesh->cell2node.father->Size() << "] " << (((mesh->getMPI().rank + 1) % 10) ? "" : "\n") << std::flush; });
            MPI::Barrier(mesh->getMPI().comm);
            if (mesh->getMPI().rank == 0)
                log() << std::endl;
            DNDS::MPISerialDo(mesh->getMPI(), [&]()
                              { log() << "[" << mesh->getMPI().rank << ": nNode " << mesh->coords.father->Size() << "] " << (((mesh->getMPI().rank + 1) % 10) ? "" : "\n") << std::flush; });
            MPI::Barrier(mesh->getMPI().comm);
            if (mesh->getMPI().rank == 0)
                log() << std::endl;
            DNDS::MPISerialDo(mesh->getMPI(), [&]()
                              { log() << "[" << mesh->getMPI().rank << ": nBnd " << mesh->bnd2node.father->Size() << "] " << (((mesh->getMPI().rank + 1) % 10) ? "" : "\n") << std::flush; });
            MPI::Barrier(mesh->getMPI().comm);
            if (mesh->getMPI().rank == 0)
                log() << std::endl;
        }
        mesh->adjPrimaryState = Adj_PointToGlobal;
        if (mesh->getMPI().rank == mRank)
            DNDS::log() << "UnstructuredMeshSerialRW === Done  PartitionReorderToMeshCell2Cell" << std::endl;
    }
 
    void UnstructuredMeshSerialRW::
        BuildSerialOut()
    {
        DNDS_assert(mesh->adjPrimaryState == Adj_PointToGlobal);
        mode = SerialOutput;
        dataIsSerialIn = false;
        dataIsSerialOut = true;
 
        std::vector<DNDS::index> serialPullCell;
        std::vector<DNDS::index> serialPullNode;
        // std::vector<DNDS::index> serialPullBnd;
 
        DNDS::index numCellGlobal = mesh->cellElemInfo.father->globalSize();
        // DNDS::index numBndGlobal = mesh->bndElemInfo.father->globalSize();
        DNDS::index numNodeGlobal = mesh->coords.father->globalSize();
 
        if (mesh->getMPI().rank == mRank)
        {
            serialPullCell.resize(numCellGlobal);
            serialPullNode.resize(numNodeGlobal);
            // serialPullBnd.reserve(numBndGlobal);
            for (DNDS::index i = 0; i < numCellGlobal; i++)
                serialPullCell[i] = i;
            for (DNDS::index i = 0; i < numNodeGlobal; i++)
                serialPullNode[i] = i;
            // for (DNDS::index i = 0; i < numBndGlobal; i++)
            //     serialPullBnd[i] = i;
        }
        cell2nodeSerial = make_ssp<decltype(cell2nodeSerial)::element_type>(ObjName{"cell2nodeSerial"}, mesh->getMPI());
        if (mesh->isPeriodic)
            cell2nodePbiSerial = make_ssp<decltype(cell2nodePbiSerial)::element_type>(ObjName{"cell2nodePbiSerial"}, NodePeriodicBits::CommType(), NodePeriodicBits::CommMult(), mesh->getMPI());
        coordSerial = make_ssp<decltype(coordSerial)::element_type>(ObjName{"coordSerial"}, mesh->getMPI());
        cellElemInfoSerial = make_ssp<decltype(cellElemInfoSerial)::element_type>(ObjName{"cellElemInfoSerial"}, ElemInfo::CommType(), ElemInfo::CommMult(), mesh->getMPI());
 
        coordSerialOutTrans.setFatherSon(mesh->coords.father, coordSerial);
        cell2nodeSerialOutTrans.setFatherSon(mesh->cell2node.father, cell2nodeSerial);
        if (mesh->isPeriodic)
            cell2nodePbiSerialOutTrans.setFatherSon(mesh->cell2nodePbi.father, cell2nodePbiSerial);
        // bnd2nodeSerialOutTrans.setFatherSon(mesh->bnd2node.father, bnd2nodeSerial);
        cellElemInfoSerialOutTrans.setFatherSon(mesh->cellElemInfo.father, cellElemInfoSerial);
        // bndElemInfoSerialOutTrans.setFatherSon(mesh->bndElemInfo.father, bndElemInfoSerial);
 
        // Father could already have global mapping, result should be the same
        coordSerialOutTrans.createFatherGlobalMapping();
        cell2nodeSerialOutTrans.createFatherGlobalMapping();
        if (mesh->isPeriodic)
            cell2nodePbiSerialOutTrans.createFatherGlobalMapping();
        // bnd2nodeSerialOutTrans.createFatherGlobalMapping();
        cellElemInfoSerialOutTrans.createFatherGlobalMapping();
        // bndElemInfoSerialOutTrans.createFatherGlobalMapping();
 
        coordSerialOutTrans.createGhostMapping(serialPullNode);
        cell2nodeSerialOutTrans.createGhostMapping(serialPullCell);
        if (mesh->isPeriodic)
            cell2nodePbiSerialOutTrans.createGhostMapping(serialPullCell);
        // bnd2nodeSerialOutTrans.createGhostMapping(serialPullBnd);
        cellElemInfoSerialOutTrans.BorrowGGIndexing(cell2nodeSerialOutTrans); // accidentally rewrites mesh->cellElemInfo.father's global mapping but ok
        // bndElemInfoSerialOutTrans.BorrowGGIndexing(bnd2nodeSerialOutTrans);
 
        coordSerialOutTrans.createMPITypes();
        cell2nodeSerialOutTrans.createMPITypes();
        if (mesh->isPeriodic)
            cell2nodePbiSerialOutTrans.createMPITypes();
        // bnd2nodeSerialOutTrans.createMPITypes();
        cellElemInfoSerialOutTrans.createMPITypes();
        // bndElemInfoSerialOutTrans.createMPITypes();
 
        coordSerialOutTrans.pullOnce();
        cell2nodeSerialOutTrans.pullOnce();
        if (mesh->isPeriodic)
            cell2nodePbiSerialOutTrans.pullOnce();
        // bnd2nodeSerialOutTrans.pullOnce();
        cellElemInfoSerialOutTrans.pullOnce();
        // bndElemInfoSerialOutTrans.pullOnce();
        if (mesh->getMPI().rank == mRank)
        {
            DNDS::log() << "UnstructuredMeshSerialRW === BuildSerialOut Done " << std::endl;
        }
    }
 
}
 
namespace DNDS::Geom
{
 
    bool UnstructuredMesh::IsO1()
    {
        using namespace Elem;
        int hasBad = 0;
        for (index iCell = 0; iCell < cellElemInfo.Size(); iCell++)
        {
            auto eType = cellElemInfo(iCell, 0).getElemType();
            if (eType == ElemType::Line2 ||
                eType == ElemType::Tri3 ||
                eType == ElemType::Quad4 ||
                eType == ElemType::Tet4 ||
                eType == ElemType::Hex8 ||
                eType == ElemType::Prism6 ||
                eType == ElemType::Pyramid5)
            {
                continue;
            }
            else
            {
                hasBad = 1;
                break;
            }
        }
        int hasBadAll;
        MPI::Allreduce(&hasBad, &hasBadAll, 1, MPI_INT, MPI_SUM, mpi.comm);
        return hasBad == 0;
    }
 
    bool UnstructuredMesh::IsO2()
    {
        using namespace Elem;
        int hasBad = 0;
        for (index iCell = 0; iCell < cellElemInfo.Size(); iCell++)
        {
            auto eType = cellElemInfo(iCell, 0).getElemType();
            if (eType == ElemType::Line3 ||
                eType == ElemType::Tri6 ||
                eType == ElemType::Quad9 ||
                eType == ElemType::Tet10 ||
                eType == ElemType::Hex27 ||
                eType == ElemType::Prism18 ||
                eType == ElemType::Pyramid14)
            {
                continue;
            }
            else
            {
                hasBad = 1;
                break;
            }
        }
        int hasBadAll;
        MPI::Allreduce(&hasBad, &hasBadAll, 1, MPI_INT, MPI_SUM, mpi.comm);
        return hasBad == 0;
    }
 
    using tIndexMapFunc = std::function<index(index)>;
 
    static void GeneralCell2NodeToNode2Cell(
        tCoordPair &coords, tAdjPair &cell2node, tAdjPair &node2cell,
        const tIndexMapFunc &CellIndexLocal2Global_NoSon,
        const tIndexMapFunc &NodeIndexLocal2Global_NoSon)
    {
        const auto &mpi = coords.father->getMPI();
        std::unordered_set<index> ghostNodesCompactSet;
        std::vector<index> ghostNodesCompact;
        std::unordered_map<index, std::unordered_set<index>> node2CellLocalRecord;
 
        for (index iCell = 0; iCell < cell2node.father->Size(); iCell++)
            for (auto iNode : cell2node[iCell])
            {
                auto [ret, rank, val] = coords.father->pLGlobalMapping->search(iNode);
                DNDS_assert_info(ret, "search failed");
                if (rank != mpi.rank)
                    ghostNodesCompact.push_back(iNode), ghostNodesCompactSet.insert(iNode);
                node2CellLocalRecord[iNode].insert(CellIndexLocal2Global_NoSon(iCell));
            }
 
        // MPI_Barrier(mpi.comm);
        // std::cout << "here2 " << std::endl;
 
        tAdj node2cellPast; // + node2cell * a triplet to deal with reverse inserting
        node2cell.InitPair("node2cell", mpi);
        node2cellPast = make_ssp<tAdj::element_type>(ObjName{"node2cellPast"}, mpi);
        //* fill into father
        node2cell.father->Resize(coords.father->Size());
        for (index iNode = 0; iNode < coords.father->Size(); iNode++)
        {
            index iNodeG = NodeIndexLocal2Global_NoSon(iNode);
            if (node2CellLocalRecord.count(iNodeG))
            {
                node2cell.ResizeRow(iNode, node2CellLocalRecord[iNodeG].size());
                rowsize in2c = 0;
                for (auto v : node2CellLocalRecord[iNodeG])
                    node2cell(iNode, in2c++) = v;
            }
        }
        node2cell.TransAttach();
        node2cell.trans.createFatherGlobalMapping();
        node2cell.trans.createGhostMapping(ghostNodesCompact);
        //* fill into son
        node2cell.son->Resize(node2cell.trans.pLGhostMapping->ghostIndex.size());
        // std::unordered_set<index> touched; // only used for checking
        for (auto &[k, s] : node2CellLocalRecord)
        {
            MPI_int rank{-1};
            index val{-1};
            if (!node2cell.trans.pLGhostMapping->search(k, rank, val))
                DNDS_assert_info(false, "search failed");
            if (rank >= 0)
            {
                node2cell.son->ResizeRow(val, s.size());
                rowsize in2c = 0;
                for (auto v : s)
                    node2cell.son->operator()(val, in2c++) = v;
                // touched.insert(val);
            }
        }
        // DNDS_assert(touched.size() == node2cell.son->Size());
 
        // node2cell.trans.pLGhostMapping->pushingIndexGlobal; // where to receive in a push
        DNDS::ArrayTransformerType<tAdj::element_type>::Type node2cellPastTrans;
        node2cellPastTrans.setFatherSon(node2cell.son, node2cellPast);
        node2cellPastTrans.createFatherGlobalMapping();
        std::vector<index> pushSonSeries(node2cell.son->Size());
        for (index i = 0; i < node2cell.son->Size(); i++)
            pushSonSeries[i] = i;
        node2cellPastTrans.createGhostMapping(pushSonSeries, node2cell.trans.pLGhostMapping->ghostStart);
        node2cellPastTrans.createMPITypes();
 
        node2cellPastTrans.pullOnce();
        DNDS_assert(DNDS::size_to_index(node2cell.trans.pLGhostMapping->ghostIndex.size()) == node2cell.son->Size());
        DNDS_assert(DNDS::size_to_index(node2cell.trans.pLGhostMapping->pushingIndexGlobal.size()) == node2cellPast->Size());
        // * this state of triplet: node2cell.father - node2cell.son - node2cellPast forms a "unique pushing" for the pair node2cell
        // * should be made into some standard
        for (index i = 0; i < node2cellPast->Size(); i++)
        {
            index iNodeG = node2cell.trans.pLGhostMapping->pushingIndexGlobal[i]; //?should be right
            for (auto iCell : (*node2cellPast)[i])
                node2CellLocalRecord[iNodeG].insert(iCell);
        }
        // MPISerialDo(
        //     mpi,
        //     [&]()
        //     {
        //         for (auto &[k, s] : node2CellLocalRecord)
        //         {
        //             if (NodeIndexGlobal2Local_NoSon(k) >= 0 && s.size() != 4)
        //                 std::cout << k << ", " << s.size() << "; " << std::flush;
        //         }
        //         std::cout << std::endl;
        //     });
 
        // reset pair
        node2cell.InitPair("node2cell", mpi);
        //* fill into father
        node2cell.father->Resize(coords.father->Size());
        for (index iNode = 0; iNode < coords.father->Size(); iNode++)
        {
            index iNodeG = NodeIndexLocal2Global_NoSon(iNode);
            if (node2CellLocalRecord.count(iNodeG))
            {
                node2cell.ResizeRow(iNode, node2CellLocalRecord[iNodeG].size());
                rowsize in2c = 0;
                for (auto v : node2CellLocalRecord[iNodeG])
                    node2cell(iNode, in2c++) = v;
            }
        }
    }
 
    void UnstructuredMesh::SetPeriodicGeometry(
        const tPoint &translation1,
        const tPoint &rotationCenter1,
        const tPoint &eulerAngles1,
        const tPoint &translation2,
        const tPoint &rotationCenter2,
        const tPoint &eulerAngles2,
        const tPoint &translation3,
        const tPoint &rotationCenter3,
        const tPoint &eulerAngles3)
    {
        periodicInfo.translation[1].map() = translation1;
        periodicInfo.translation[2].map() = translation2;
        periodicInfo.translation[3].map() = translation3;
        periodicInfo.rotationCenter[1].map() = rotationCenter1;
        periodicInfo.rotationCenter[2].map() = rotationCenter2;
        periodicInfo.rotationCenter[3].map() = rotationCenter3;
        periodicInfo.rotation[1].map() =
            Geom::RotZ(eulerAngles1[2]) *
            Geom::RotY(eulerAngles1[1]) *
            Geom::RotX(eulerAngles1[0]);
        periodicInfo.rotation[2].map() =
            Geom::RotZ(eulerAngles2[2]) *
            Geom::RotY(eulerAngles2[1]) *
            Geom::RotX(eulerAngles2[0]);
        periodicInfo.rotation[3].map() =
            Geom::RotZ(eulerAngles3[2]) *
            Geom::RotY(eulerAngles3[1]) *
            Geom::RotX(eulerAngles3[0]);
    }
 
    void UnstructuredMesh::
        RecoverNode2CellAndNode2Bnd()
    {
        DNDS_assert(adjPrimaryState == Adj_PointToGlobal);
        DNDS_assert(coords.father);
        DNDS_assert(cell2node.father);
        DNDS_assert(bnd2node.father);
 
        /*****************************************************/
        // * first recover node2cell
 
        if (!coords.father->pLGlobalMapping)
            coords.father->createGlobalMapping(); // for NodeIndexLocal2Global_NoSon
        if (!cell2node.father->pLGlobalMapping)
            cell2node.father->createGlobalMapping(); // for CellIndexLocal2Global_NoSon
 
        GeneralCell2NodeToNode2Cell(
            coords, cell2node, node2cell,
            [this](index v)
            { return this->CellIndexLocal2Global_NoSon(v); },
            [this](index v)
            { return this->NodeIndexLocal2Global_NoSon(v); });
 
        if (!bnd2node.father->pLGlobalMapping)
            bnd2node.father->createGlobalMapping(); // for BndIndexLocal2Global_NoSon
        GeneralCell2NodeToNode2Cell(
            coords, bnd2node, node2bnd,
            [this](index v)
            { return this->BndIndexLocal2Global_NoSon(v); },
            [this](index v)
            { return this->NodeIndexLocal2Global_NoSon(v); });
 
        this->adjN2CBState = Adj_PointToGlobal;
 
        // if (mpi.rank == 0)
        // {
        //     for (index i = 0; i < node2cell.father->Size(); i++)
        //         std::cout << node2cell.RowSize(i) - 4 << std::endl;
        //     for (index i = 0; i < node2bnd.father->Size(); i++)
        //         std::cout << node2bnd.RowSize(i) + 10 << std::endl;
        // }
 
        // node2cell.TransAttach();
        // node2cell.trans.createFatherGlobalMapping();
        // node2cell.trans.createGhostMapping(ghostNodesCompact);
        // node2cell.trans.createMPITypes();
        // node2cell.trans.pullOnce();
 
        // mesh->node2cell.TransAttach();
        // mesh->node2cell.trans.BorrowGGIndexing(mesh->coords);
        // mesh->node2cell.trans.createMPITypes();
        // mesh->node2cell.trans.pullOnce();
    }
 
    void UnstructuredMesh::RecoverCell2CellAndBnd2Cell()
    {
        DNDS_assert(adjPrimaryState == Adj_PointToGlobal);
        DNDS_assert(adjN2CBState == Adj_PointToGlobal);
        DNDS_assert(coords.father);
        DNDS_assert(cell2node.father);
        DNDS_assert(bnd2node.father);
        DNDS_assert(node2cell.father);
 
        coords.TransAttach();
        coords.trans.createFatherGlobalMapping(); // for NodeIndexLocal2Global_NoSon
        cell2node.TransAttach();
        cell2node.trans.createFatherGlobalMapping(); // for CellIndexLocal2Global_NoSon
        bnd2node.TransAttach();
        bnd2node.trans.createFatherGlobalMapping(); // for BndIndexLocal2Global_NoSon
 
        // std::cout << "RecoverCell2CellAndBnd2Cell here1" << std::endl;
 
        std::unordered_set<index> ghostNodesCompactSet;
        for (index i = 0; i < cell2node.father->Size(); i++)
            for (auto in : cell2node.father->operator[](i))
                if (NodeIndexGlobal2Local_NoSon(in) < 0)
                    ghostNodesCompactSet.insert(in);
        for (index i = 0; i < bnd2node.father->Size(); i++)
            for (auto in : bnd2node.father->operator[](i))
                if (NodeIndexGlobal2Local_NoSon(in) < 0)
                    ghostNodesCompactSet.insert(in);
        std::vector<index> ghostNodes;
        ghostNodes.reserve(ghostNodesCompactSet.size());
        for (auto v : ghostNodesCompactSet)
            ghostNodes.push_back(v);
        // std::cout << "RecoverCell2CellAndBnd2Cell here2" << std::endl;
        node2cell.son = make_ssp<decltype(node2cell.son)::element_type>(ObjName{"node2cell.son"}, mpi);
        node2cell.TransAttach();
        node2cell.trans.createFatherGlobalMapping();
        node2cell.trans.createGhostMapping(ghostNodes);
        node2cell.trans.createMPITypes(); //! warning, this is not actual final official trans, just needed temporarily
        node2cell.trans.pullOnce();
 
        std::unordered_map<index, index> iNodeGlobal2LocalAppendInNode2Cell;
 
        for (index i = 0; i < node2cell.Size(); i++)
            iNodeGlobal2LocalAppendInNode2Cell[node2cell.trans.pLGhostMapping->operator()(-1, i)] = i;
 
        for (index i = 0; i < cell2node.father->Size(); i++)
            for (auto in : cell2node.father->operator[](i))
                DNDS_assert(iNodeGlobal2LocalAppendInNode2Cell.count(in));
 
        cell2cell.InitPair("cell2cell", mpi); // actual outputs need empty but constructed son
        cell2cell.father->Resize(cell2node.father->Size());
        for (index i = 0; i < cell2node.father->Size(); i++)
        {
            std::set<index> cellRec;
            for (auto in : cell2node.father->operator[](i))
            {
                DNDS_assert(iNodeGlobal2LocalAppendInNode2Cell.count(in));
                for (auto ico : node2cell[iNodeGlobal2LocalAppendInNode2Cell.at(in)])
                    cellRec.insert(ico);
            }
            auto ret = cellRec.erase(CellIndexLocal2Global_NoSon(i));
            DNDS_assert(ret == 1);
            cell2cell.father->ResizeRow(i, cellRec.size());
            rowsize ic2c = 0;
            for (auto v : cellRec)
                cell2cell.father->operator()(i, ic2c++) = v;
        }
        // std::cout << "RecoverCell2CellAndBnd2Cell here2.5" << std::endl;
        bnd2cell.InitPair("bnd2cell", mpi); // actual outputs need empty but constructed son
        bnd2cell.father->Resize(bnd2node.father->Size());
 
        // For periodic meshes, store per-bnd candidate cell sets from the node
        // intersection pass, then resolve via pbi check in a second pass after
        // ghost-pulling cell2node/cell2nodePbi for the candidate cells.
        std::vector<std::set<index>> bndCellCandidates;
        if (isPeriodic)
            bndCellCandidates.resize(bnd2node.father->Size());
 
        for (index i = 0; i < bnd2node.father->Size(); i++)
        {
            std::set<index> cellRecCur;
            bool initDone{false};
            // std::cout << "RecoverCell2CellAndBnd2Cell here L1 1" << std::endl;
            for (auto in : bnd2node.father->operator[](i))
            {
                DNDS_assert(iNodeGlobal2LocalAppendInNode2Cell.count(in));
                std::set<index> cellRecCurNew;
                for (auto ico : node2cell[iNodeGlobal2LocalAppendInNode2Cell.at(in)])
                    if (!initDone || cellRecCur.count(ico))
                        cellRecCurNew.insert(ico);
                std::swap(cellRecCur, cellRecCurNew);
                initDone = true;
            }
 
            // Periodic pbi check and bnd2cell assignment are deferred to the
            // second pass below (after ghost-pulling cell2node/cell2nodePbi for
            // the candidate cells). Store candidates per bnd for now.
            if (isPeriodic)
                bndCellCandidates[i] = std::move(cellRecCur);
            else
            {
                DNDS_assert_info(cellRecCur.size() == 1, fmt::format("cellRecCur.size() is [{}]", cellRecCur.size()));
                bnd2cell.father->operator()(i, 0) = *cellRecCur.begin();
                bnd2cell.father->operator()(i, 1) = UnInitIndex;
            }
        }
 
        /************************************************************/
        // Periodic: ghost-pull cell2node and cell2nodePbi for all candidate
        // cells, then do the pbi filter and donor/receiver assignment.
        /************************************************************/
        if (isPeriodic)
        {
            // Collect all unique candidate cells across all periodic bnds
            std::vector<index> neededCells;
            for (index i = 0; i < bnd2cell.father->Size(); i++)
            {
                if (!Geom::FaceIDIsPeriodic(bndElemInfo.father->operator()(i, 0).zone))
                    continue;
                for (auto ic : bndCellCandidates[i])
                    neededCells.push_back(ic);
            }
 
            cell2nodePbi.son = make_ssp<decltype(cell2nodePbi.son)::element_type>(ObjName{"cell2nodePbi.son"}, NodePeriodicBits::CommType(), NodePeriodicBits::CommMult(), mpi);
            cell2nodePbi.TransAttach();
            cell2nodePbi.trans.createFatherGlobalMapping();
            cell2nodePbi.trans.createGhostMapping(neededCells); //! warning, this is not actual final official trans, just needed temporarily
            cell2nodePbi.trans.createMPITypes();
            cell2nodePbi.trans.pullOnce();
            cell2node.son = make_ssp<decltype(cell2node.son)::element_type>(ObjName{"cell2node.son"}, mpi);
            cell2node.BorrowAndPull(cell2nodePbi); //! warning, this is not actual final official trans, just needed temporarily
 
            // Now do the periodic pbi filter for each bnd
            for (index i = 0; i < bnd2cell.father->Size(); i++)
            {
                if (!Geom::FaceIDIsPeriodic(bndElemInfo.father->operator()(i, 0).zone))
                    continue;
 
                auto &cellRecCur = bndCellCandidates[i]; // the pre-pbi candidates (from node intersection)
                std::set<index> cellRecFiltered;
                for (auto ic : cellRecCur)
                {
                    auto [ret, rank, icAppend] = cell2nodePbi.trans.pLGhostMapping->search_indexAppend(ic);
                    DNDS_assert_info(ret, fmt::format("periodic bnd {} candidate cell {} not found in ghost mapping", i, ic));
 
                    bool cellContainsBnd = true;
                    for (int ib2n = 0; ib2n < bnd2node.father->operator[](i).size(); ib2n++)
                    {
                        index iNode = bnd2node.father->operator[](i)[ib2n];
                        auto iNodePbi = bnd2nodePbi.father->operator[](i)[ib2n];
                        int nIndexMatchNode{0}, nIndexPBIMatchNode{0};
                        for (int ic2n = 0; ic2n < cell2node[icAppend].size(); ic2n++)
                            if (iNode == cell2node(icAppend, ic2n))
                            {
                                nIndexMatchNode++;
                                if (iNodePbi == cell2nodePbi(icAppend, ic2n))
                                    nIndexPBIMatchNode++;
                            }
                        DNDS_assert(nIndexMatchNode >= 1);
                        if (nIndexPBIMatchNode == 0)
                            cellContainsBnd = false;
                    }
                    if (cellContainsBnd)
                        cellRecFiltered.insert(ic);
                }
 
                // Periodic bnd: 2 cells or 1 (self-periodic)
                DNDS_assert_info(cellRecFiltered.size() == 2 || cellRecFiltered.size() == 1,
                                 fmt::format("periodic bnd {} has {} cells after pbi filter", i, cellRecFiltered.size()));
                auto it = cellRecFiltered.begin();
                bnd2cell.father->operator()(i, 0) = *it;
                if (cellRecFiltered.size() == 2)
                    ++it;
                bnd2cell.father->operator()(i, 1) = *it;
            }
 
            // Swap check: ensure bnd2cell(i, 0) is the donor-side cell
            for (index i = 0; i < bnd2cell.father->Size(); i++)
                if (bnd2cell(i, 1) != UnInitIndex && bnd2cell(i, 0) != bnd2cell(i, 1) /* no need to check if both sides are same*/)
                {
                    index ic0 = bnd2cell(i, 0);
                    index ic1 = bnd2cell(i, 1);
                    auto [ret0, rank0, ic0L] = cell2nodePbi.trans.pLGhostMapping->search_indexAppend(ic0);
                    auto [ret1, rank1, ic1L] = cell2nodePbi.trans.pLGhostMapping->search_indexAppend(ic1);
                    // std::cout << fmt::format("ic0L ic1L, {},{}", ic0L, ic1L) << std::endl;
                    DNDS_assert_info(ret0 && ret1, "search failed");
                    std::vector<Geom::NodePeriodicBits> pbi0, pbi1;
 
                    for (auto in : bnd2node[i])
                    {
                        bool found0{false}, found1{false};
                        for (rowsize ic2n = 0; ic2n < cell2node[ic0L].size(); ic2n++)
                            if (cell2node[ic0L][ic2n] == in)
                                found0 = true, pbi0.push_back(cell2nodePbi(ic0L, ic2n));
                        for (rowsize ic2n = 0; ic2n < cell2node[ic1L].size(); ic2n++)
                            if (cell2node[ic1L][ic2n] == in)
                                found1 = true, pbi1.push_back(cell2nodePbi(ic1L, ic2n));
                        DNDS_assert(found0 && found1);
                    }
                    auto cleanOtherBits = [&](Geom::NodePeriodicBits &b)
                    {
                        if (Geom::FaceIDIsPeriodic1(bndElemInfo(i, 0).zone))
                            b = b & nodePB1;
                        else if (Geom::FaceIDIsPeriodic2(bndElemInfo(i, 0).zone))
                            b = b & nodePB2;
                        else if (Geom::FaceIDIsPeriodic3(bndElemInfo(i, 0).zone))
                            b = b & nodePB3;
                    };
                    for (auto &b : pbi0)
                        cleanOtherBits(b);
                    for (auto &b : pbi1)
                        cleanOtherBits(b);
 
                    Geom::NodePeriodicBits bndBit;
                    if (bndElemInfo(i, 0).zone == Geom::BC_ID_PERIODIC_1_DONOR)
                        bndBit.setP1True();
                    else if (bndElemInfo(i, 0).zone == Geom::BC_ID_PERIODIC_2_DONOR)
                        bndBit.setP2True();
                    else if (bndElemInfo(i, 0).zone == Geom::BC_ID_PERIODIC_3_DONOR)
                        bndBit.setP3True();
                    else if (!Geom::FaceIDIsPeriodic(bndElemInfo(i, 0).zone))
                        DNDS_assert_info(false, "this bnd with both sides has to be periodic");
 
                    bool match0{true}, match1{true};
                    for (auto b : pbi0)
                        if (b ^ bndBit)
                            match0 = false;
                    for (auto b : pbi1)
                        if (b ^ bndBit)
                            match1 = false;
                    if (match0 && !match1)
                        ; // keep
                    else if (match1 && !match0)
                        std::swap(bnd2cell(i, 0), bnd2cell(i, 1));
                    else
                    {
                        if (mpi.rank >= 0)
                        {
                            for (auto b : pbi0)
                                DNDS::log() << b << ", ";
                            DNDS::log() << " ---- ";
                            for (auto b : pbi1)
                                DNDS::log() << b << ", ";
                            DNDS::log() << " ---- " << bndBit;
                            DNDS::log() << std::endl;
                        }
                        DNDS_assert_info(false,
                                         match0
                                             ? "this periodic bnd matches both sides"
                                             : "this periodic bnd matches no sides");
                    }
                    // if (mpi.rank == 0)
                    // {
                    //     if (match1 && !match0)
                    //         std::cout << "swap " << i << std::endl;
                    // }
                }
        }
    }
 
    void UnstructuredMesh::
        BuildGhostPrimary()
    {
        DNDS_assert(adjPrimaryState == Adj_PointToGlobal);
        DNDS_assert(cell2cell.father && cell2cell.father->Size() == this->NumCell());
        DNDS_assert(bnd2cell.father && bnd2cell.father->Size() == this->NumBnd());
        /********************************/
        // cells
        {
            cell2cell.TransAttach();
            cell2node.TransAttach();
            cell2cellOrig.TransAttach();
            if (isPeriodic)
                cell2nodePbi.TransAttach();
            cellElemInfo.TransAttach();
 
            cell2cell.trans.createFatherGlobalMapping();
 
            std::vector<DNDS::index> ghostCells;
            for (DNDS::index iCell = 0; iCell < cell2cell.father->Size(); iCell++)
            {
                for (DNDS::rowsize ic2c = 0; ic2c < cell2cell.father->RowSize(iCell); ic2c++)
                {
                    auto iCellOther = (*cell2cell.father)(iCell, ic2c);
                    DNDS::MPI_int rank;
                    DNDS::index val;
                    if (!cell2cell.trans.pLGlobalMapping->search(iCellOther, rank, val))
                        DNDS_assert_info(false, "search failed");
                    if (rank != mpi.rank)
                        ghostCells.push_back(iCellOther);
                }
            }
            cell2cell.trans.createGhostMapping(ghostCells);
 
            cell2cell.trans.createMPITypes();
            cell2cell.trans.pullOnce();
            cell2node.BorrowAndPull(cell2cell);
            cell2cellOrig.BorrowAndPull(cell2cell);
            if (isPeriodic)
                cell2nodePbi.BorrowAndPull(cell2cell);
            cellElemInfo.BorrowAndPull(cell2cell);
        }
        // if(mpi.rank == 0)
        //     std::cout <<"XXXXXXXXXXXXXXXXXXXXXXXXX" <<std::endl;
        // if(mpi.rank == 0)
        //     for(index iC = 0; iC < coords.father->Size(); iC ++)
        //         std::cout << coords.father->operator[](iC).transpose() << std::endl;
 
        /********************************/
        // cells done, go on to nodes
        {
            coords.TransAttach();
            node2nodeOrig.TransAttach();
 
            coords.trans.createFatherGlobalMapping();
 
            std::vector<DNDS::index> ghostNodes;
            for (DNDS::index iCell = 0; iCell < cell2cell.Size(); iCell++) // note doing full (son + father) traverse
            {
                for (DNDS::rowsize ic2c = 0; ic2c < cell2node.RowSize(iCell); ic2c++)
                {
                    auto iNode = cell2node(iCell, ic2c);
                    DNDS::MPI_int rank;
                    DNDS::index val;
                    if (!coords.trans.pLGlobalMapping->search(iNode, rank, val))
                        DNDS_assert_info(false, "search failed");
                    if (rank != mpi.rank)
                        ghostNodes.push_back(iNode);
                }
            }
            coords.trans.createGhostMapping(ghostNodes);
            coords.trans.createMPITypes();
            coords.trans.pullOnce();
            node2nodeOrig.BorrowAndPull(coords);
        }
        // if(mpi.rank == 0)
        //     std::cout <<"XXXXXXXXXXXXXXXXXXXXXXXXX" <<std::endl;
        // if(mpi.rank == 0)
        //     for(index iC = 0; iC < coords.Size(); iC ++)
        //         std::cout << coords.operator[](iC).transpose() << std::endl;
 
        /********************************/
        // bnds: added via node2bnd's father part
        {
            DNDS_assert(node2bnd.father);
            DNDS_assert(this->adjN2CBState == Adj_PointToGlobal);
            bnd2cell.TransAttach();
            bnd2node.TransAttach();
            if (isPeriodic)
                bnd2nodePbi.TransAttach();
            bndElemInfo.TransAttach();
            bnd2bndOrig.TransAttach();
 
            bnd2cell.trans.createFatherGlobalMapping();
 
            std::vector<DNDS::index> ghostBnds;
            for (index iNode = 0; iNode < node2bnd.father->Size(); iNode++)
                for (auto iBnd : node2bnd[iNode])
                {
                    auto [ret, rank, val] = bnd2cell.trans.pLGlobalMapping->search(iBnd);
                    DNDS_assert_info(ret, "search failed");
                    if (rank != mpi.rank)
                        ghostBnds.push_back(iBnd);
                }
            bnd2cell.trans.createGhostMapping(ghostBnds);
            bnd2cell.trans.createMPITypes();
            bnd2cell.trans.pullOnce();
            bnd2node.BorrowAndPull(bnd2cell);
            if (isPeriodic)
                bnd2nodePbi.BorrowAndPull(bnd2cell);
            bndElemInfo.BorrowAndPull(bnd2cell);
            bnd2bndOrig.BorrowAndPull(bnd2cell);
 
            // Ghost bnds may reference nodes not yet in the coord ghost layer.
            // Add those nodes so that AdjGlobal2LocalPrimary can convert bnd2node.
            // This is collective: all ranks must participate even if some have no extras.
            {
                std::vector<DNDS::index> extraGhostNodes;
                for (DNDS::index iBnd = bnd2node.father->Size(); iBnd < bnd2node.Size(); iBnd++)
                    for (DNDS::rowsize j = 0; j < bnd2node.RowSize(iBnd); j++)
                    {
                        auto iNode = bnd2node(iBnd, j);
                        DNDS::MPI_int rank;
                        DNDS::index val;
                        if (!coords.trans.pLGhostMapping->search_indexAppend(iNode, rank, val))
                            extraGhostNodes.push_back(iNode);
                    }
                DNDS::index nExtraLocal = extraGhostNodes.size();
                DNDS::index nExtraGlobal{0};
                MPI::Allreduce(&nExtraLocal, &nExtraGlobal, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
                if (nExtraGlobal > 0)
                {
                    // Rebuild coord ghost mapping with the additional nodes
                    auto &existingGhost = coords.trans.pLGhostMapping->ghostIndex;
                    std::vector<DNDS::index> allGhostNodes(existingGhost.begin(), existingGhost.end());
                    allGhostNodes.insert(allGhostNodes.end(), extraGhostNodes.begin(), extraGhostNodes.end());
                    coords.trans.createGhostMapping(allGhostNodes);
                    node2nodeOrig.trans.BorrowGGIndexing(coords.trans);
                    coords.trans.createMPITypes();
                    node2nodeOrig.trans.createMPITypes();
                    coords.trans.pullOnce();
                    node2nodeOrig.trans.pullOnce();
                }
            }
        }
    }
 
    void UnstructuredMesh::
        AdjGlobal2LocalPrimary()
    {
        // needs results of BuildGhostPrimary()
        DNDS_assert(adjPrimaryState == Adj_PointToGlobal);
        DNDS_assert_info(cellElemInfo.trans.pLGhostMapping, "trans of cellElemInfo needed but not built");
        DNDS_assert_info(coords.trans.pLGhostMapping, "trans of coords needed but not built");
 
        /**********************************/
        // convert bnd2cell, bnd2node, cell2cell, cell2node ptrs global to local
        ConvertAdjEntries(cell2cell, cell2cell.Size(),
                          [&](index v) { return CellIndexGlobal2Local(v); });
 
        for (DNDS::index iBnd = 0; iBnd < bnd2cell.Size(); iBnd++)
            for (DNDS::rowsize j = 0; j < bnd2cell.RowSize(iBnd); j++)
                bnd2cell(iBnd, j) = CellIndexGlobal2Local(bnd2cell(iBnd, j)),
                               DNDS_assert(j == 0 ? (iBnd < bnd2cell.father->Size() ? bnd2cell(iBnd, j) >= 0 : true) : true); // father bnds' owner cell must be inside
 
        ConvertAdjEntries(cell2node, cell2node.Size(),
                          [&](index v) { auto r = NodeIndexGlobal2Local(v); DNDS_assert(r >= 0); return r; });
 
        ConvertAdjEntries(bnd2node, bnd2node.Size(),
                          [&](index v) { auto r = NodeIndexGlobal2Local(v); DNDS_assert(r >= 0); return r; });
        /**********************************/
        adjPrimaryState = Adj_PointToLocal;
    }
 
    void UnstructuredMesh::
        AdjLocal2GlobalPrimary()
    {
        DNDS_assert(adjPrimaryState == Adj_PointToLocal);
        DNDS_assert_info(cellElemInfo.trans.pLGhostMapping, "trans of cellElemInfo needed but not built");
        DNDS_assert_info(coords.trans.pLGhostMapping, "trans of coords needed but not built");
 
        /**********************************/
        // convert bnd2cell, bnd2node, cell2cell, cell2node ptrs local to global
        ConvertAdjEntries(cell2cell, cell2cell.Size(),
                          [&](index v) { return CellIndexLocal2Global(v); });
 
        for (DNDS::index iBnd = 0; iBnd < bnd2cell.Size(); iBnd++)
            for (DNDS::rowsize j = 0; j < bnd2cell.RowSize(iBnd); j++)
                bnd2cell(iBnd, j) = CellIndexLocal2Global(bnd2cell(iBnd, j)),
                               DNDS_assert(j == 0 ? (iBnd < bnd2cell.father->Size() ? bnd2cell(iBnd, j) >= 0 : true) : true); // father bnds' owner cell must be inside
 
        ConvertAdjEntries(cell2node, cell2node.Size(),
                          [&](index v) { auto r = NodeIndexLocal2Global(v); DNDS_assert(r >= 0); return r; });
 
        ConvertAdjEntries(bnd2node, bnd2node.Size(),
                          [&](index v) { auto r = NodeIndexLocal2Global(v); DNDS_assert(r >= 0); return r; });
        /**********************************/
        adjPrimaryState = Adj_PointToGlobal;
    }
 
    void UnstructuredMesh::
        AdjGlobal2LocalPrimaryForBnd() // a reduction of primary version
    {
        // needs results of BuildGhostPrimary()
        DNDS_assert(adjPrimaryState == Adj_PointToGlobal);
        /**********************************/
        // convert cell2node ptrs global to local
        ConvertAdjEntries(cell2node, cell2node.Size(),
                          [&](index v) { auto r = NodeIndexGlobal2Local(v); DNDS_assert(r >= 0); return r; });
        /**********************************/
        adjPrimaryState = Adj_PointToLocal;
    }
 
    void UnstructuredMesh::
        AdjLocal2GlobalPrimaryForBnd() // a reduction of primary version
    {
        DNDS_assert(adjPrimaryState == Adj_PointToLocal);
        /**********************************/
        // convert cell2node ptrs local to global
        ConvertAdjEntries(cell2node, cell2node.Size(),
                          [&](index v) { auto r = NodeIndexLocal2Global(v); DNDS_assert(r >= 0); return r; });
        /**********************************/
        adjPrimaryState = Adj_PointToGlobal;
    }
 
    void UnstructuredMesh::
        AdjGlobal2LocalFacial()
    {
        DNDS_assert(adjFacialState == Adj_PointToGlobal);
        /**********************************/
        // convert face2cell ptrs and face2node ptrs global to local
#ifdef DNDS_USE_OMP
#    pragma omp parallel for
#endif
        for (DNDS::index iFace = 0; iFace < face2cell.Size(); iFace++)
        {
            for (rowsize if2n = 0; if2n < face2node.RowSize(iFace); if2n++)
            {
                index &iNode = face2node(iFace, if2n);
                index val;
                MPI_int rank;
                auto ret = coords.trans.pLGhostMapping->search_indexAppend(iNode, rank, val);
                DNDS_assert(ret);
                iNode = val;
            }
            for (rowsize if2c = 0; if2c < 2; if2c++)
            {
                index &iCell = face2cell(iFace, if2c);
                if (iCell != UnInitIndex) // is not a bnd
                {
                    index val;
                    MPI_int rank;
                    auto ret = cell2node.trans.pLGhostMapping->search_indexAppend(iCell, rank, val);
                    DNDS_assert(ret);
                    iCell = val;
                }
            }
        }
        /**********************************/
        adjFacialState = Adj_PointToLocal;
    }
 
    void UnstructuredMesh::
        AdjLocal2GlobalFacial()
    {
        DNDS_assert(adjFacialState == Adj_PointToLocal);
        /**********************************/
        // convert face2cell ptrs and face2node ptrs to global
#ifdef DNDS_USE_OMP
#    pragma omp parallel for
#endif
        for (DNDS::index iFace = 0; iFace < face2cell.Size(); iFace++)
        {
            for (rowsize if2n = 0; if2n < face2node.RowSize(iFace); if2n++)
            {
                index &iNode = face2node(iFace, if2n);
                iNode = coords.trans.pLGhostMapping->operator()(-1, iNode);
            }
            for (rowsize if2c = 0; if2c < 2; if2c++)
            {
                index &iCell = face2cell(iFace, if2c);
                if (iCell != UnInitIndex) // is not a bnd
                    iCell = cell2node.trans.pLGhostMapping->operator()(-1, iCell);
            }
            index &iBnd = face2bnd(iFace);
            iBnd = this->BndIndexGlobal2Local(iBnd);
        }
        // MPI::Barrier(mpi.comm);
        /**********************************/
        adjFacialState = Adj_PointToGlobal;
    }
 
    void UnstructuredMesh::
        AdjLocal2GlobalC2F()
    {
        DNDS_assert(adjC2FState == Adj_PointToLocal);
        /**********************************/
        // convert cell2face
 
        auto FaceIndexLocal2Global = [&](DNDS::index &iF)
        {
            DNDS_assert(face2node.trans.pLGhostMapping);
            if (iF == UnInitIndex)
                return;
            if (iF < 0) // mapping to un-found in father-son
                iF = -1 - iF;
            else
                iF = face2node.trans.pLGhostMapping->operator()(-1, iF);
        };
 
#ifdef DNDS_USE_OMP
#    pragma omp parallel for
#endif
        for (DNDS::index iCell = 0; iCell < cell2face.Size(); iCell++)
        {
            for (rowsize ic2f = 0; ic2f < cell2face.RowSize(iCell); ic2f++)
            {
                index &iFace = cell2face(iCell, ic2f);
                FaceIndexLocal2Global(iFace);
            }
        }
        for (index iBnd = 0; iBnd < bnd2face.Size(); iBnd++)
            for (auto &iFace : bnd2face[iBnd])
                FaceIndexLocal2Global(iFace);
        // MPI::Barrier(mpi.comm);
        /**********************************/
        adjC2FState = Adj_PointToGlobal;
    }
 
    void UnstructuredMesh::
        AdjGlobal2LocalC2F()
    {
        DNDS_assert(adjC2FState == Adj_PointToGlobal);
        /**********************************/
        // convert cell2face
 
        auto FaceIndexGlobal2Local = [&](DNDS::index &iF)
        {
            DNDS_assert(face2node.trans.pLGhostMapping);
            if (iF == UnInitIndex)
                return;
            DNDS::MPI_int rank;
            DNDS::index val;
            auto result = face2node.trans.pLGhostMapping->search_indexAppend(iF, rank, val);
            if (result)
                iF = val;
            else
                iF = -1 - iF; // mapping to un-found in father-son
        };
 
#ifdef DNDS_USE_OMP
#    pragma omp parallel for
#endif
        for (DNDS::index iCell = 0; iCell < cell2face.Size(); iCell++)
        {
            for (rowsize ic2f = 0; ic2f < cell2face.RowSize(iCell); ic2f++)
            {
                index &iFace = cell2face(iCell, ic2f);
                FaceIndexGlobal2Local(iFace);
            }
        }
        for (index iBnd = 0; iBnd < bnd2face.Size(); iBnd++)
            for (auto &iFace : bnd2face[iBnd])
                FaceIndexGlobal2Local(iFace);
        /**********************************/
        adjC2FState = Adj_PointToLocal;
    }
 
    void UnstructuredMesh::
        BuildGhostN2CB()
    {
        DNDS_assert(adjN2CBState == Adj_PointToGlobal);
 
        DNDS_assert(coords.trans.father && coords.trans.pLGhostMapping);
 
        node2cell.TransAttach();
        node2cell.trans.BorrowGGIndexing(coords.trans);
        node2cell.trans.createMPITypes();
        node2cell.trans.pullOnce();
 
        node2bnd.TransAttach();
        node2bnd.trans.BorrowGGIndexing(coords.trans);
        node2bnd.trans.createMPITypes();
        node2bnd.trans.pullOnce();
    }
 
    void UnstructuredMesh::
        AdjLocal2GlobalN2CB()
    {
        DNDS_assert(adjN2CBState == Adj_PointToLocal);
        DNDS_assert_info(cellElemInfo.trans.pLGhostMapping, "trans of cellElemInfo needed but not built");
        DNDS_assert_info(bndElemInfo.trans.pLGhostMapping, "trans of bndElemInfo needed but not built");
        /**********************************/
#ifdef DNDS_USE_OMP
#    pragma omp parallel for
#endif
        for (index iNode = 0; iNode < node2cell.Size(); iNode++)
            for (index &iCell : node2cell[iNode])
                iCell = CellIndexLocal2Global(iCell);
 
#ifdef DNDS_USE_OMP
#    pragma omp parallel for
#endif
        for (index iNode = 0; iNode < node2bnd.Size(); iNode++)
            for (index &iBnd : node2bnd[iNode])
                iBnd = BndIndexLocal2Global(iBnd); // ! bnd has son now
 
        // MPI::Barrier(mpi.comm);
        /**********************************/
        adjN2CBState = Adj_PointToGlobal;
    }
 
    void UnstructuredMesh::
        AdjGlobal2LocalN2CB()
    {
        DNDS_assert(adjN2CBState == Adj_PointToGlobal);
        DNDS_assert_info(cellElemInfo.trans.pLGhostMapping, "trans of cellElemInfo needed but not built");
        DNDS_assert_info(bndElemInfo.trans.pLGhostMapping, "trans of bndElemInfo needed but not built");
        /**********************************/
// todo: ensure using dist omp settings not single!
#ifdef DNDS_USE_OMP
#    pragma omp parallel for
#endif
        for (index iNode = 0; iNode < node2cell.Size(); iNode++)
            for (index &iCell : node2cell[iNode])
                iCell = CellIndexGlobal2Local(iCell);
 
#ifdef DNDS_USE_OMP
#    pragma omp parallel for
#endif
        for (index iNode = 0; iNode < node2bnd.Size(); iNode++)
            for (index &iBnd : node2bnd[iNode])
                iBnd = BndIndexGlobal2Local(iBnd); // ! bnd has son now
 
        // MPI::Barrier(mpi.comm);
        /**********************************/
        adjN2CBState = Adj_PointToLocal;
    }
 
    void UnstructuredMesh::AssertOnN2CB()
    {
        for (index iNode = NumNode(); iNode < NumNodeProc(); iNode++)
        {
            int nCellAdjIn = 0;
            for (auto iCell : node2cell[iNode])
                if (iCell >= 0)
                {
                    DNDS_assert(iCell < NumCellProc());
                    nCellAdjIn++;
                    for (auto iNodeOther : cell2node[iCell])
                    {
                        DNDS_assert(iNodeOther < NumNodeProc() && iNodeOther >= 0);
                    }
                }
            DNDS_assert(nCellAdjIn);
        }
 
        std::set<index> bnd_main_nodes;
        for (index iNode = 0; iNode < NumNode(); iNode++)
        {
            for (auto iCell : node2cell[iNode])
            {
                DNDS_assert(iCell < NumCellProc() && iCell >= 0);
                for (auto iNodeOther : cell2node[iCell])
                {
                    DNDS_assert(iNodeOther < NumNodeProc() && iNodeOther >= 0);
                    if (iNodeOther >= NumNode())
                        bnd_main_nodes.insert(iNode);
                }
            }
        }
        std::map<index, int> bnd_main_nodes_adj_ghost_num;
        for (index iNode : bnd_main_nodes)
            bnd_main_nodes_adj_ghost_num[iNode] = 0;
        for (index iNode = NumNode(); iNode < NumNodeProc(); iNode++)
            for (auto iCell : node2cell[iNode])
                if (iCell >= 0)
                    for (auto iNodeOther : cell2node[iCell])
                        if (iNodeOther >= 0 && iNodeOther < NumNode())
                            bnd_main_nodes_adj_ghost_num.at(iNodeOther)++;
        for (auto [iNode, num_ghost_adj] : bnd_main_nodes_adj_ghost_num)
            DNDS_assert(num_ghost_adj);
    }
 
    void UnstructuredMesh::BuildCell2CellFace()
    {
        DNDS_assert(adjPrimaryState == Adj_PointToLocal);
        DNDS_assert(adjFacialState == Adj_PointToLocal);
        DNDS_assert_info(cellElemInfo.trans.pLGhostMapping, "trans of cellElemInfo needed but not built");
 
        cell2cellFace.InitPair("cell2cellFace", mpi);
        cell2cellFace.father->Resize(this->NumCell());
        for (index iCell = 0; iCell < this->NumCell(); iCell++)
        {
            cell2cellFace.ResizeRow(iCell, cell2face[iCell].size());
            for (rowsize ic2f = 0; ic2f < cell2face[iCell].size(); ic2f++)
            {
                index iFace = cell2face[iCell][ic2f];
                index iCellOther = this->CellFaceOther(iCell, iFace);
                DNDS_assert(iCellOther < this->NumCellProc());
                cell2cellFace[iCell][ic2f] = this->CellIndexLocal2Global(iCellOther);
            }
        }
        cell2cellFace.father->Compress();
        cell2cellFace.TransAttach();
        cell2cellFace.trans.BorrowGGIndexing(cell2node.trans);
        cell2cellFace.trans.createMPITypes();
        cell2cellFace.trans.pullOnce(); // warning! to pull the adj, must be in global state!
        adjC2CFaceState = Adj_PointToGlobal;
    }
 
    void UnstructuredMesh::AdjLocal2GlobalC2CFace()
    {
        // needs results of BuildGhostPrimary()
        DNDS_assert(adjC2CFaceState == Adj_PointToLocal);
        DNDS_assert_info(cellElemInfo.trans.pLGhostMapping, "trans of cellElemInfo needed but not built");
 
        /**********************************/
        ConvertAdjEntries(cell2cellFace, cell2cellFace.Size(),
                          [&](index v) { return CellIndexLocal2Global(v); });
        /**********************************/
 
        adjC2CFaceState = Adj_PointToGlobal;
    }
 
    void UnstructuredMesh::AdjGlobal2LocalC2CFace()
    {
        // needs results of BuildGhostPrimary()
        DNDS_assert(adjC2CFaceState == Adj_PointToGlobal);
        DNDS_assert_info(cellElemInfo.trans.pLGhostMapping, "trans of cellElemInfo needed but not built");
 
        /**********************************/
        ConvertAdjEntries(cell2cellFace, cell2cellFace.Size(),
                          [&](index v) { return CellIndexGlobal2Local(v); });
        /**********************************/
 
        adjC2CFaceState = Adj_PointToLocal;
    }
 
    // =================================================================
    // File-local helpers for InterpolateFace
    // =================================================================
    namespace
    {
        /**
         * \brief Result of face enumeration from cell connectivity.
         */
        struct FaceEnumerationResult
        {
            std::vector<std::vector<DNDS::index>> face2nodeV;
            std::vector<std::vector<NodePeriodicBits>> face2nodePbiV;
            std::vector<std::pair<DNDS::index, DNDS::index>> face2cellV;
            std::vector<ElemInfo> faceElemInfoV;
            DNDS::index nFaces = 0;
        };
 
        /**
         * \brief Enumerate unique faces from cell-to-node connectivity.
         *
         * Iterates over all cells (local + ghost), extracts topological faces,
         * deduplicates by sorted vertex comparison (with periodic-aware matching
         * when \p isPeriodic is true), and populates \p cell2face entries.
         *
         * \param[in,out] cell2face     Cell-to-face pair (rows resized and entries set).
         * \param[in]     cellElemInfo  Cell element info (for element type).
         * \param[in]     cell2node     Cell-to-node connectivity.
         * \param[in]     cell2nodePbi  Cell-to-node periodic bits (used when isPeriodic).
         * \param[in]     cell2cellSize Total number of cells (father + son).
         * \param[in]     coordsSize    Total number of nodes (father + son), for node2face sizing.
         * \param[in]     isPeriodic    Whether periodic mesh handling is enabled.
         */
        FaceEnumerationResult EnumerateFacesFromCells(
            tAdjPair &cell2face,
            const tElemInfoArrayPair &cellElemInfo,
            const tAdjPair &cell2node,
            const tPbiPair &cell2nodePbi,
            DNDS::index cell2cellSize,
            DNDS::index coordsSize,
            bool isPeriodic)
        {
            FaceEnumerationResult result;
            auto &face2nodeV = result.face2nodeV;
            auto &face2nodePbiV = result.face2nodePbiV;
            auto &face2cellV = result.face2cellV;
            auto &faceElemInfoV = result.faceElemInfoV;
            auto &nFaces = result.nFaces;
 
            std::vector<std::vector<DNDS::index>> node2face(coordsSize);
 
            using idx_pbi_pair = std::pair<index, NodePeriodicBits>;
            auto idx_pbi_pair_less = [](idx_pbi_pair &L, idx_pbi_pair &R)
            { return L.first == R.first ? uint8_t(L.second) < uint8_t(R.second) : L.first < R.first; };
            auto idx_pbi_pair_eq = [](idx_pbi_pair &L, idx_pbi_pair &R)
            { return L.first == R.first && uint8_t(L.second) == uint8_t(R.second); };
 
            for (DNDS::index iCell = 0; iCell < cell2cellSize; iCell++)
            {
                auto eCell = Elem::Element{cellElemInfo[iCell]->getElemType()};
                cell2face.ResizeRow(iCell, eCell.GetNumFaces());
                for (int ic2f = 0; ic2f < eCell.GetNumFaces(); ic2f++)
                {
                    auto eFace = eCell.ObtainFace(ic2f);
                    std::vector<DNDS::index> faceNodes(eFace.GetNumNodes());
                    eCell.ExtractFaceNodes(ic2f, cell2node[iCell], faceNodes);
                    DNDS::index iFound = -1;
                    std::vector<DNDS::index> faceVerts(faceNodes.begin(), faceNodes.begin() + eFace.GetNumVertices());
                    std::sort(faceVerts.begin(), faceVerts.end());
 
                    std::vector<NodePeriodicBits> faceNodePeriodicBits;
                    std::vector<idx_pbi_pair> faceNodeNodePeriodicBits;
                    if (isPeriodic)
                    {
                        faceNodePeriodicBits.resize(eFace.GetNumNodes());
                        faceNodeNodePeriodicBits.resize(eFace.GetNumNodes());
                        eCell.ExtractFaceNodes(ic2f, cell2nodePbi[iCell], faceNodePeriodicBits);
                        for (size_t i = 0; i < faceNodeNodePeriodicBits.size(); i++)
                            faceNodeNodePeriodicBits[i].first = faceNodes[i], faceNodeNodePeriodicBits[i].second = faceNodePeriodicBits[i];
                        std::sort(faceNodeNodePeriodicBits.begin(), faceNodeNodePeriodicBits.end(),
                                  idx_pbi_pair_less);
                        for (size_t i = 1; i < faceNodeNodePeriodicBits.size(); i++)
                        {
                            DNDS_assert_info(!idx_pbi_pair_eq(faceNodeNodePeriodicBits[i - 1], faceNodeNodePeriodicBits[i]),
                                             "the face has identical (periodic) nodes");
                        }
                    }
 
                    for (auto iV : faceVerts)
                        if (iFound < 0)
                            for (auto iFOther : node2face[iV])
                            {
                                auto eFaceOther = Elem::Element{faceElemInfoV[iFOther].getElemType()};
                                if (eFaceOther.type != eFace.type)
                                    continue;
                                std::vector<DNDS::index> faceVertsOther(
                                    face2nodeV[iFOther].begin(),
                                    face2nodeV[iFOther].begin() + eFace.GetNumVertices());
                                std::sort(faceVertsOther.begin(), faceVertsOther.end());
                                if (std::equal(faceVerts.begin(), faceVerts.end(), faceVertsOther.begin(), faceVertsOther.end()))
                                {
                                    if (isPeriodic)
                                    {
                                        std::vector<std::pair<index, NodePeriodicBits>> faceNodeNodePeriodicBitsOther(eFaceOther.GetNumNodes());
                                        for (size_t i = 0; i < faceNodeNodePeriodicBitsOther.size(); i++)
                                            faceNodeNodePeriodicBitsOther[i].first = face2nodeV[iFOther][i],
                                            faceNodeNodePeriodicBitsOther[i].second = face2nodePbiV[iFOther][i];
                                        std::sort(faceNodeNodePeriodicBitsOther.begin(), faceNodeNodePeriodicBitsOther.end(),
                                                  idx_pbi_pair_less);
                                        auto v0 = faceNodeNodePeriodicBits.at(0).second ^ faceNodeNodePeriodicBitsOther.at(0).second;
                                        DNDS_assert(faceNodeNodePeriodicBitsOther.size() == faceNodeNodePeriodicBits.size());
                                        bool collaborating = true;
                                        for (size_t i = 1; i < faceNodeNodePeriodicBitsOther.size(); i++)
                                            if ((faceNodeNodePeriodicBits[i].second ^ faceNodeNodePeriodicBitsOther[i].second) != v0)
                                                collaborating = false;
                                        if (collaborating)
                                            iFound = iFOther;
                                    }
                                    else
                                        iFound = iFOther;
                                }
                            }
                    if (iFound < 0)
                    {
                        // * face not existent yet
                        face2nodeV.emplace_back(faceNodes); // note: faceVerts invalid here!
                        if (isPeriodic)
                            face2nodePbiV.emplace_back(faceNodePeriodicBits);
                        face2cellV.emplace_back(std::make_pair(iCell, DNDS::UnInitIndex));
                        // important note: f2nPbi node pbi is always same as cell f2c[0]'s corresponding nodes
                        faceElemInfoV.emplace_back(ElemInfo{eFace.type, 0});
                        for (auto iV : faceVerts)
                            node2face[iV].push_back(nFaces);
                        cell2face(iCell, ic2f) = nFaces;
                        nFaces++;
                    }
                    else
                    {
                        DNDS_assert(face2cellV[iFound].second == DNDS::UnInitIndex);
                        face2cellV[iFound].second = iCell;
                        cell2face(iCell, ic2f) = iFound;
                    }
                }
            }
            return result;
        }
 
        /**
         * \brief Result of face collection/filtering pass.
         */
        struct FaceCollectionResult
        {
            std::vector<index> iFaceAllToCollected;  ///< mapping: old face idx -> collected idx (or UnInitIndex/-1)
            std::vector<std::vector<index>> faceSendLocals; ///< per-rank lists of local face indices to send as ghost
            index nFacesNew = 0;
        };
 
        /**
         * \brief Filter faces: discard ghost-only and duplicate cross-rank faces,
         *        assign ownership to the rank with the lower rank ID.
         *
         * \param[in] faceElemInfoV     Per-face element info (zone for internal/bnd distinction).
         * \param[in] face2cellV        Per-face cell pair (first, second).
         * \param[in] nFaces            Total enumerated faces.
         * \param[in] nLocalCells       Number of local (father) cells.
         * \param[in] pLGhostMapping    Ghost mapping from cell2node.trans (for global index lookup).
         * \param[in] pLGlobalMapping   Global mapping from cell2node.father (for rank search).
         * \param[in] nMPIRanks         Number of MPI ranks.
         * \param[in] myRank            This process's MPI rank.
         */
        FaceCollectionResult CollectFaces(
            const std::vector<ElemInfo> &faceElemInfoV,
            const std::vector<std::pair<DNDS::index, DNDS::index>> &face2cellV,
            DNDS::index nFaces,
            DNDS::index nLocalCells,
            const DNDS::OffsetAscendIndexMapping &ghostMapping,
            const DNDS::GlobalOffsetsMapping &globalMapping,
            DNDS::MPI_int nMPIRanks,
            DNDS::MPI_int myRank)
        {
            FaceCollectionResult result;
            result.iFaceAllToCollected.resize(nFaces);
            result.faceSendLocals.resize(nMPIRanks);
            result.nFacesNew = 0;
 
            for (index iFace = 0; iFace < nFaces; iFace++)
            {
                if (faceElemInfoV[iFace].zone <= 0) // if internal
                {
                    if (face2cellV[iFace].second == UnInitIndex && face2cellV[iFace].first >= nLocalCells) // has not other cell with ghost parent
                        result.iFaceAllToCollected[iFace] = UnInitIndex;                                   // * discard
                    else if (face2cellV[iFace].first >= nLocalCells &&
                             face2cellV[iFace].second >= nLocalCells) // both sides ghost
                        result.iFaceAllToCollected[iFace] = UnInitIndex; // * discard
                    else if (face2cellV[iFace].first >= nLocalCells ||
                             face2cellV[iFace].second >= nLocalCells)
                    {
                        DNDS_assert(face2cellV[iFace].second >= nLocalCells); // should only be the internal as first
                        // * check both sided's info //TODO: optimize so that pLGhostMapping returns rank directly ?
                        index cellGlobL = ghostMapping(-1, face2cellV[iFace].first);
                        index cellGlobR = ghostMapping(-1, face2cellV[iFace].second);
                        MPI_int rankL, rankR;
                        index valL, valR;
                        auto retL = globalMapping.search(cellGlobL, rankL, valL);
                        auto retR = globalMapping.search(cellGlobR, rankR, valR);
                        DNDS_assert(retL && retR && (rankL != rankR));
                        if (rankL > rankR)
                        {
                            result.iFaceAllToCollected[iFace] = -1; // * discard but with ghost
                        }
                        else
                        {
                            DNDS_assert(rankL == myRank);
                            result.faceSendLocals[rankR].push_back(result.nFacesNew);
                            result.iFaceAllToCollected[iFace] = result.nFacesNew++; //*use
                        }
                    }
                    else
                    {
                        result.iFaceAllToCollected[iFace] = result.nFacesNew++; //*use
                    }
                }
                else // all bnds would be non duplicate
                {
                    result.iFaceAllToCollected[iFace] = result.nFacesNew++; //*use
                }
            }
            return result;
        }
 
        /**
         * \brief Copy collected (non-discarded) face data from temporary vectors
         *        into the resized face member arrays, then remap cell2face entries.
         *
         * \param[in]  faceEnum           Face enumeration result (temp vectors).
         * \param[in]  faceCollect        Face collection result (mapping + nFacesNew).
         * \param[out] face2cell          Face-to-cell pair (father resized and filled).
         * \param[out] face2node          Face-to-node pair (father resized and filled).
         * \param[out] face2nodePbi       Face-to-node periodic bits pair (if periodic).
         * \param[out] faceElemInfo        Face element info pair (father resized and filled).
         * \param[in,out] cell2face       Cell-to-face pair (entries remapped via iFaceAllToCollected).
         * \param[in]  isPeriodic         Whether periodic mesh handling is enabled.
         * \param[in]  mpiComm            MPI communicator for barrier.
         */
        void CompactFacesAndRemapCell2Face(
            const FaceEnumerationResult &faceEnum,
            const FaceCollectionResult &faceCollect,
            tAdj2Pair &face2cell,
            tAdjPair &face2node,
            tPbiPair &face2nodePbi,
            tElemInfoArrayPair &faceElemInfo,
            tAdjPair &cell2face,
            bool isPeriodic,
            MPI_Comm mpiComm)
        {
            const auto &face2nodeV = faceEnum.face2nodeV;
            const auto &face2nodePbiV = faceEnum.face2nodePbiV;
            const auto &face2cellV = faceEnum.face2cellV;
            const auto &faceElemInfoV = faceEnum.faceElemInfoV;
            const auto &iFaceAllToCollected = faceCollect.iFaceAllToCollected;
            index nFacesNew = faceCollect.nFacesNew;
            index nFaces = faceEnum.nFaces;
 
            face2cell.father->Resize(nFacesNew);
            face2node.father->Resize(nFacesNew);
            if (isPeriodic)
                face2nodePbi.father->Resize(nFacesNew);
            faceElemInfo.father->Resize(nFacesNew); //! considering globally duplicate faces
            nFacesNew = 0;
            for (DNDS::index iFace = 0; iFace < nFaces; iFace++)
            {
                if (iFaceAllToCollected[iFace] >= 0) // ! -1 is also ignored!
                {
                    face2node.ResizeRow(nFacesNew, face2nodeV[iFace].size());
                    for (DNDS::rowsize if2n = 0; if2n < face2node.RowSize(nFacesNew); if2n++)
                        face2node(nFacesNew, if2n) = face2nodeV[iFace][if2n];
                    if (isPeriodic)
                    {
                        DNDS_assert(face2nodeV[iFace].size() == face2nodePbiV[iFace].size());
                        face2nodePbi.ResizeRow(nFacesNew, face2nodePbiV[iFace].size());
                        for (DNDS::rowsize if2n = 0; if2n < face2nodePbi.RowSize(nFacesNew); if2n++)
                            face2nodePbi(nFacesNew, if2n) = face2nodePbiV[iFace][if2n];
                    }
                    face2cell(nFacesNew, 0) = face2cellV[iFace].first;
                    face2cell(nFacesNew, 1) = face2cellV[iFace].second;
                    faceElemInfo(nFacesNew, 0) = faceElemInfoV[iFace];
                    nFacesNew++;
                }
            }
 
            MPI::Barrier(mpiComm);
#ifdef DNDS_USE_OMP
#    pragma omp parallel for
#endif
            for (DNDS::index iCell = 0; iCell < cell2face.Size(); iCell++) // convert face indices pointers
            {
                for (rowsize ic2f = 0; ic2f < cell2face.RowSize(iCell); ic2f++)
                {
                    cell2face(iCell, ic2f) = iFaceAllToCollected[cell2face(iCell, ic2f)]; // Uninit if to discard
                }
            }
        }
 
        /**
         * \brief Match boundary elements to faces, populating faceElemInfo zone IDs
         *        and building bnd2face / face2bnd mappings.
         *
         * For each boundary, finds the face that shares the same node set via
         * the parent cell's cell2face connectivity.
         */
        void MatchBoundariesToFaces(
            const tElemInfoArrayPair &bndElemInfo,
            const tAdj2Pair &bnd2cell,
            const tAdjPair &bnd2node,
            const tAdjPair &cell2face,
            const tAdjPair &face2node,
            const tAdjPair &cell2node,
            tElemInfoArrayPair &faceElemInfo,
            std::vector<index> &bnd2faceV,
            std::unordered_map<index, index> &face2bndM,
            tAdj1Pair &face2bnd,
            tAdj1Pair &bnd2face)
        {
            bnd2faceV.resize(bndElemInfo.father->Size(), -1); // this mapping only uses main (father) part
            face2bndM.reserve(bndElemInfo.father->Size());
            std::unordered_map<index, index> iFace2iBnd;
            for (DNDS::index iBnd = 0; iBnd < bndElemInfo.father->Size(); iBnd++)
            {
                DNDS::index pCell = bnd2cell(iBnd, 0);
                std::vector<DNDS::index> b2nRow = bnd2node[iBnd];
                std::sort(b2nRow.begin(), b2nRow.end());
                int nFound = 0;
                auto faceID = bndElemInfo[iBnd]->zone;
                for (int ic2f = 0; ic2f < cell2face.RowSize(pCell); ic2f++)
                {
                    auto iFace = cell2face(pCell, ic2f);
                    if (iFace < 0) //==-1, pointing to ghost face
                        continue;
                    std::vector<DNDS::index> f2nRow = face2node[iFace];
                    std::sort(f2nRow.begin(), f2nRow.end());
                    if (std::equal(b2nRow.begin(), b2nRow.end(), f2nRow.begin(), f2nRow.end()))
                    {
                        if (iFace2iBnd.count(iFace))
                        {
                            DNDS_assert(FaceIDIsPeriodic(faceID)); // only periodic gets to be duplicated
                            index iBndOther = iFace2iBnd[iFace];
                            index iCellA = bnd2cell[iBnd][0];
                            index iCellB = bnd2cell[iBndOther][0];
                            DNDS_assert(iCellA < cell2node.father->Size()); // both points to local non-ghost cells
                            DNDS_assert(iCellB < cell2node.father->Size()); // both points to local non-ghost cells
                            // std::cout << iCellA << " vs " << iCellB << std::endl;
                            if (iCellA > iCellB) // make face corresponds with f2c[0]'s bnd if possible
                                continue;
                        }
                        iFace2iBnd[iFace] = iBnd;
                        nFound++; // two things:
                        // if is periodic, then only gets the bnd info of the main cell's bnd;
                        // if is external bc, then must be non-ghost face
                        faceElemInfo(iFace, 0) = bndElemInfo(iBnd, 0);
                        bnd2faceV[iBnd] = iFace;
                        face2bndM[iFace] = iBnd;
                        DNDS_assert_info(FaceIDIsExternalBC(faceID) ||
                                             FaceIDIsPeriodic(faceID),
                                         "bnd elem should have a BC id not interior");
                    }
                }
                DNDS_assert(nFound > 0 || (FaceIDIsPeriodic(faceID) && nFound == 0)); // periodic could miss the face
            }
 
            face2bnd.father->Resize(face2node.father->Size());
            for (index iFace = 0; iFace < face2bnd.father->Size(); iFace++)
                face2bnd.father->operator()(iFace) = face2bndM.count(iFace) ? face2bndM[iFace] : UnInitIndex;
            bnd2face.father->Resize(bnd2node.father->Size());
            for (index iBnd = 0; iBnd < bnd2node.father->Size(); iBnd++)
                bnd2face.father->operator()(iBnd) = bnd2faceV.at(iBnd);
        }
 
        /**
         * \brief Match received ghost faces to cells and assign cell2face entries
         *        that were previously marked -1 (pointing to ghost).
         *
         * For each ghost face, iterates over both its cells, finds the matching
         * topological face slot in cell2face, and assigns the ghost face index.
         */
        void AssignGhostFacesToCells(
            const tAdj2 &face2cellSon,
            const tAdj &face2nodeSon,
            const tElemInfoArray &faceElemInfoSon,
            tAdjPair &cell2face,
            const tAdjPair &cell2node,
            const tElemInfoArrayPair &cellElemInfo,
            DNDS::index nFatherFaces)
        {
#ifdef DNDS_USE_OMP
#    pragma omp parallel for
#endif
            for (DNDS::index iFace = 0; iFace < face2cellSon->Size(); iFace++) // face2cell points to local now
            {
                // before: first points to inner, //!relies on the order of setting face2cell
                DNDS_assert((*face2cellSon)(iFace, 0) >= cell2node.father->Size());
                auto eFace = Elem::Element{(*faceElemInfoSon)(iFace, 0).getElemType()};
                auto faceVerts = std::vector<index>((*face2nodeSon)[iFace].begin(), (*face2nodeSon)[iFace].begin() + eFace.GetNumVertices());
                std::sort(faceVerts.begin(), faceVerts.end()); //* do not forget to do set operation sort first
                for (rowsize if2c = 0; if2c < 2; if2c++)
                {
                    index iCell = (*face2cellSon)(iFace, if2c);
                    auto cell2faceRow = cell2face[iCell];
                    auto cellNodes = cell2node[iCell];
                    auto eCell = Elem::Element{cellElemInfo(iCell, 0).getElemType()};
                    bool found = false;
                    for (rowsize ic2f = 0; ic2f < cell2face.RowSize(iCell); ic2f++)
                    {
                        auto eFace = eCell.ObtainFace(ic2f);
                        std::vector<index> faceNodesC(eFace.GetNumNodes());
                        eCell.ExtractFaceNodes(ic2f, cellNodes, faceNodesC);
                        std::sort(faceNodesC.begin(), faceNodesC.end());
                        if (std::includes(faceNodesC.begin(), faceNodesC.end(), faceVerts.begin(), faceVerts.end()))
                        {
                            DNDS_assert(cell2face(iCell, ic2f) == -1);
                            cell2face(iCell, ic2f) = iFace + nFatherFaces; // remember is ghost
                            found = true;
                        }
                    }
                    DNDS_assert(found);
                }
            }
        }
    } // anonymous namespace
 
    /// @todo //TODO: handle periodic cases
    void UnstructuredMesh::
        InterpolateFace()
    {
        DNDS_assert(adjPrimaryState == Adj_PointToLocal); // And also should have primary ghost comm
 
        // Allocate face-related array pairs
        cell2face.InitPair("cell2face", mpi);
        face2cell.InitPair("face2cell", mpi);
        face2node.InitPair("face2node", mpi);
        if (isPeriodic)
            face2nodePbi.InitPair("face2nodePbi", mpi);
        faceElemInfo.InitPair("faceElemInfo", mpi);
        face2bnd.InitPair("face2bnd", mpi);
        bnd2face.InitPair("bnd2face", mpi);
 
        cell2face.father->Resize(cell2cell.father->Size());
        cell2face.son->Resize(cell2cell.son->Size());
 
        // Section B: Enumerate unique faces from cell connectivity
        auto faceEnum = EnumerateFacesFromCells(
            cell2face, cellElemInfo, cell2node, cell2nodePbi,
            cell2cell.Size(), coords.Size(), isPeriodic);
 
        // Section C: Filter faces — discard ghost-only and duplicate cross-rank
        auto faceCollect = CollectFaces(
            faceEnum.faceElemInfoV, faceEnum.face2cellV, faceEnum.nFaces,
            cell2face.father->Size(),
            *cell2node.trans.pLGhostMapping, *cell2node.father->pLGlobalMapping,
            mpi.size, mpi.rank);
 
        // Section D: Compact collected faces into member arrays, remap cell2face
        CompactFacesAndRemapCell2Face(
            faceEnum, faceCollect,
            face2cell, face2node, face2nodePbi, faceElemInfo,
            cell2face, isPeriodic, mpi.comm);
        adjFacialState = Adj_PointToLocal;
 
        // Section E: Match boundary elements to faces
        MatchBoundariesToFaces(
            bndElemInfo, bnd2cell, bnd2node, cell2face, face2node, cell2node,
            faceElemInfo, bnd2faceV, face2bndM, face2bnd, bnd2face);
 
        // Section F: Ghost face communication
        this->AdjLocal2GlobalFacial();
 
        // Flatten faceSendLocals into CSR format
        std::vector<index> faceSendLocalsIdx;
        std::vector<index> faceSendLocalsStarts(mpi.size + 1);
        faceSendLocalsStarts[0] = 0;
        for (MPI_int r = 0; r < mpi.size; r++)
            faceSendLocalsStarts[r + 1] = faceSendLocalsStarts[r] + faceCollect.faceSendLocals[r].size();
        faceSendLocalsIdx.resize(faceSendLocalsStarts.back());
        for (MPI_int r = 0; r < mpi.size; r++)
            std::copy(faceCollect.faceSendLocals[r].begin(), faceCollect.faceSendLocals[r].end(),
                      faceSendLocalsIdx.begin() + faceSendLocalsStarts[r]);
 
        face2node.father->Compress(); // before comm
        if (isPeriodic)
            face2nodePbi.father->Compress();
        face2cell.TransAttach();
        face2cell.trans.createFatherGlobalMapping();
        face2cell.trans.createGhostMapping(faceSendLocalsIdx, faceSendLocalsStarts);
        face2cell.trans.createMPITypes();
        face2cell.trans.pullOnce();
        face2node.BorrowAndPull(face2cell);
        if (isPeriodic)
            face2nodePbi.BorrowAndPull(face2cell);
        faceElemInfo.BorrowAndPull(face2cell);
        face2bnd.BorrowAndPull(face2cell);
 
        this->AdjGlobal2LocalFacial();
 
        // Section G: Assign ghost faces to cell2face entries marked -1
        AssignGhostFacesToCells(
            face2cell.son, face2node.son, faceElemInfo.son,
            cell2face, cell2node, cellElemInfo,
            face2cell.father->Size());
 
        cell2face.father->Compress();
        cell2face.son->Compress();
        adjC2FState = Adj_PointToLocal;
 
        // Section H: Communicate cell2face and bnd2face ghost data
        this->AdjLocal2GlobalC2F();
        cell2face.TransAttach();
        cell2face.trans.BorrowGGIndexing(cell2node.trans);
        cell2face.trans.createMPITypes();
        cell2face.trans.pullOnce();
        bnd2face.TransAttach();
        bnd2face.trans.BorrowGGIndexing(bnd2node.trans);
        bnd2face.trans.createMPITypes();
        bnd2face.trans.pullOnce();
        this->AdjGlobal2LocalC2F();
 
        for (DNDS::index iFace = 0; iFace < faceElemInfo.Size(); iFace++)
        {
            if (FaceIDIsPeriodicMain(faceElemInfo(iFace, 0).zone))
            {
                // DNDS_assert(false);
            }
        }
        for (DNDS::index iFace = 0; iFace < faceElemInfo.Size(); iFace++)
        {
            if (FaceIDIsPeriodicDonor(faceElemInfo(iFace, 0).zone))
            {
                // DNDS_assert(false);
            }
        }
 
        auto gSize = face2node.father->globalSize(); //! sync call!!!
        if (mpi.rank == 0)
            log() << "UnstructuredMesh === InterpolateFace: total faces " << gSize << std::endl;
    }
 
    void UnstructuredMesh::
        AssertOnFaces()
    {
 
        //* some assertions on faces
        std::vector<uint16_t> cCont(cell2cell.Size(), 0); // simulate flux
        for (DNDS::index iFace = 0; iFace < faceElemInfo.Size(); iFace++)
        {
            auto faceID = faceElemInfo(iFace, 0).zone;
            if (FaceIDIsInternal(faceID))
            {
                // if (FaceIDIsPeriodic(faceID))
                // {
                //     // TODO: tend to the case of face is PeriodicDonor with Main in same proc
                //     continue;
                // }
                // if (face2cell[iFace][0] < cell2cell.father->Size()) // other side prime cell, periodic also
                DNDS_assert_info(face2cell[iFace][1] != DNDS::UnInitIndex,
                                 fmt::format(
                                     "Face {} is internal, but f2c[1] is null, at {},{},{} - {},{},{}", iFace,
                                     coords[face2node[iFace][0]](0),
                                     coords[face2node[iFace][0]](1),
                                     coords[face2node[iFace][0]](2),
                                     face2node[iFace].size() > 1 ? coords[face2node[iFace][1]](0) : 0.,
                                     face2node[iFace].size() > 1 ? coords[face2node[iFace][1]](1) : 0.,
                                     face2node[iFace].size() > 1 ? coords[face2node[iFace][1]](2) : 0.)); // Assert has enough cell donors
                DNDS_assert(face2cell[iFace][0] >= 0 && face2cell[iFace][0] < cell2cell.Size());
                DNDS_assert(face2cell[iFace][1] >= 0 && face2cell[iFace][1] < cell2cell.Size());
                cCont[face2cell[iFace][0]]++;
                cCont[face2cell[iFace][1]]++;
            }
            else // a external BC
            {
                DNDS_assert(face2cell[iFace][1] == DNDS::UnInitIndex);
                DNDS_assert(face2cell[iFace][0] >= 0 && face2cell[iFace][0] < cell2cell.father->Size());
                cCont[face2cell[iFace][0]]++;
            }
        }
        for (DNDS::index iCell = 0; iCell < cellElemInfo.father->Size(); iCell++) // for every non-ghost
        {
            for (auto iFace : cell2face[iCell])
            {
                DNDS_assert(iFace >= 0 && iFace < face2cell.Size());
                DNDS_assert(face2cell[iFace][0] == iCell || face2cell[iFace][1] == iCell);
            }
            DNDS_assert(cCont[iCell] == cell2face.RowSize(iCell));
        }
    }
 
    void UnstructuredMesh::
        WriteSerialize(Serializer::SerializerBaseSSP serializerP, const std::string &name)
    {
        // TODO: Make these read/write of mesh independent of ghost part!
        DNDS_assert(adjPrimaryState == Adj_PointToGlobal);
 
        auto cwd = serializerP->GetCurrentPath();
        serializerP->CreatePath(name);
        serializerP->GoToPath(name);
 
        serializerP->WriteString("mesh", "UnstructuredMesh");
        serializerP->WriteIndex("dim", dim);
        if (serializerP->IsPerRank())
            serializerP->WriteIndex("MPIRank", mpi.rank);
        serializerP->WriteIndex("MPISize", mpi.size);
        serializerP->WriteInt("isPeriodic", isPeriodic);
 
        coords.WriteSerialize(serializerP, "coords");
        cell2node.WriteSerialize(serializerP, "cell2node");
        // cell2cell.WriteSerialize(serializerP, "cell2cell");
        cellElemInfo.WriteSerialize(serializerP, "cellElemInfo");
        bnd2node.WriteSerialize(serializerP, "bnd2node");
        // bnd2cell.WriteSerialize(serializerP, "bnd2cell");
        bndElemInfo.WriteSerialize(serializerP, "bndElemInfo");
        if (isPeriodic)
        {
            cell2nodePbi.WriteSerialize(serializerP, "cell2nodePbi");
            bnd2nodePbi.WriteSerialize(serializerP, "bnd2nodePbi");
            periodicInfo.WriteSerializer(serializerP, "periodicInfo");
        }
        bnd2bndOrig.WriteSerialize(serializerP, "bnd2bndOrig");
        cell2cellOrig.WriteSerialize(serializerP, "cell2cellOrig");
        node2nodeOrig.WriteSerialize(serializerP, "node2nodeOrig");
 
        serializerP->GoToPath(cwd);
    }
 
    void UnstructuredMesh::
        ReadSerialize(Serializer::SerializerBaseSSP serializerP, const std::string &name)
    {
        auto cwd = serializerP->GetCurrentPath();
        // serializerP->CreatePath(name);//! remember no create!
        serializerP->GoToPath(name);
 
        std::string meshRead;
        index dimRead{0}, rankRead{0}, sizeRead{0};
        int isPeriodicRead;
        serializerP->ReadString("mesh", meshRead);
        serializerP->ReadIndex("dim", dimRead);
        if (serializerP->IsPerRank())
            serializerP->ReadIndex("MPIRank", rankRead);
        serializerP->ReadIndex("MPISize", sizeRead);
        serializerP->ReadInt("isPeriodic", isPeriodicRead);
        isPeriodic = bool(isPeriodicRead);
        DNDS_assert(meshRead == "UnstructuredMesh");
        DNDS_assert(dimRead == dim);
        DNDS_assert((!serializerP->IsPerRank() || rankRead == mpi.rank) && sizeRead == mpi.size);
 
        // make the empty arrays
        coords.InitPair("coords", getMPI());
        cellElemInfo.InitPair("cellElemInfo", getMPI());
        bndElemInfo.InitPair("bndElemInfo", getMPI());
        cell2node.InitPair("cell2node", getMPI());
        if (isPeriodic)
        {
            cell2nodePbi.InitPair("cell2nodePbi", getMPI());
            bnd2nodePbi.InitPair("bnd2nodePbi", getMPI());
        }
        bnd2node.InitPair("bnd2node", getMPI());
 
        bnd2bndOrig.InitPair("bnd2bndOrig", getMPI());
        cell2cellOrig.InitPair("cell2cellOrig", getMPI());
        node2nodeOrig.InitPair("node2nodeOrig", getMPI());
 
        coords.ReadSerialize(serializerP, "coords");
        cell2node.ReadSerialize(serializerP, "cell2node");
        // cell2cell.ReadSerialize(serializerP, "cell2cell");
        cellElemInfo.ReadSerialize(serializerP, "cellElemInfo");
        bnd2node.ReadSerialize(serializerP, "bnd2node");
        // bnd2cell.ReadSerialize(serializerP, "bnd2cell");
        bndElemInfo.ReadSerialize(serializerP, "bndElemInfo");
        if (isPeriodic)
        {
            cell2nodePbi.ReadSerialize(serializerP, "cell2nodePbi");
            bnd2nodePbi.ReadSerialize(serializerP, "bnd2nodePbi");
            periodicInfo.ReadSerializer(serializerP, "periodicInfo");
        }
        bnd2bndOrig.ReadSerialize(serializerP, "bnd2bndOrig");
        cell2cellOrig.ReadSerialize(serializerP, "cell2cellOrig");
        node2nodeOrig.ReadSerialize(serializerP, "node2nodeOrig");
 
        // after matters:
        coords.trans.createMPITypes();
        cell2node.trans.createMPITypes();
        // cell2cell.trans.createMPITypes();
        cellElemInfo.trans.createMPITypes();
        bnd2node.trans.createMPITypes();
        // bnd2cell.trans.createMPITypes();
        bndElemInfo.trans.createMPITypes();
        if (isPeriodic)
        {
            cell2nodePbi.trans.createMPITypes();
            bnd2nodePbi.trans.createMPITypes();
        }
        bnd2bndOrig.trans.createMPITypes();
        cell2cellOrig.trans.createMPITypes();
        node2nodeOrig.trans.createMPITypes();
        adjPrimaryState = Adj_PointToGlobal; // the file is pointing to local
 
        index nCellG = this->NumCellGlobal(); // collective call!
        index nNodeG = this->NumNodeGlobal(); // collective call!
        index nNodeB = this->NumBndGlobal();  // collective call!
        if (mpi.rank == mRank)
        {
            log() << "UnstructuredMesh === ReadSerialize "
                  << "Global NumCell [ " << nCellG << " ]" << std::endl;
            log() << "UnstructuredMesh === ReadSerialize "
                  << "Global NumNode [ " << nNodeG << " ]" << std::endl;
            log() << "UnstructuredMesh === ReadSerialize "
                  << "Global NumBnd  [ " << nNodeB << " ]" << std::endl;
        }
        MPISerialDo(mpi, [&]()
                    { log() << "    Rank: " << mpi.rank << " nCell " << this->NumCell() << " nCellGhost " << this->NumCellGhost() << std::endl; });
        MPISerialDo(mpi, [&]()
                    { log() << "    Rank: " << mpi.rank << " nNode " << this->NumNode() << " nNodeGhost " << this->NumNodeGhost() << std::endl; });
 
        serializerP->GoToPath(cwd);
    }
 
    void UnstructuredMesh::ConstructBndMesh(UnstructuredMesh &bMesh)
    {
        DNDS_assert(bMesh.dim == dim - 1 && bMesh.mpi == mpi);
        bMesh.cell2node.InitPair("bMesh.cell2node", mpi);
        bMesh.coords.InitPair("bMesh.coords", mpi);
        if (isPeriodic)
        {
            bMesh.isPeriodic = true;
            bMesh.periodicInfo = this->periodicInfo;
            bMesh.cell2nodePbi.InitPair("bMesh.cell2nodePbi", mpi);
        }
        bMesh.cellElemInfo.InitPair("bMesh.cellElemInfo", mpi);
        bMesh.cell2cellOrig.InitPair("bMesh.cell2cellOrig", mpi);
        bMesh.node2nodeOrig.InitPair("bMesh.node2nodeOrig", mpi);
 
        bMesh.bnd2cell.InitPair("bMesh.bnd2cell", mpi);       // which will remain 0 sized
        bMesh.bnd2node.InitPair("bMesh.bnd2node", mpi);       // which will remain 0 sized
        bMesh.bndElemInfo.InitPair("bMesh.bndElemInfo", mpi);  // which will remain 0 sized
        bMesh.bnd2bndOrig.InitPair("bMesh.bnd2bndOrig", mpi); // which will remain 0 sized
        if (isPeriodic)
            bMesh.bnd2nodePbi.InitPair("bMesh.bnd2nodePbi", mpi); // which will remain 0 sized
 
        tIndPair node2bndNodeGlobal;
        node2bndNodeGlobal.InitPair("node2bndNodeGlobal", mpi);
 
        node2bndNodeGlobal.father->Resize(this->NumNode());
        node2bndNodeGlobal.TransAttach();
        node2bndNodeGlobal.trans.BorrowGGIndexing(coords.trans);
        node2bndNodeGlobal.trans.createMPITypes();
        index bndNodeCount{0}, bndNodeStart{0};
        {
            for (index i = 0; i < node2bndNodeGlobal.Size(); i++)
                node2bndNodeGlobal[i] = 0;
            for (index iBnd = 0; iBnd < this->NumBnd(); iBnd++)
                for (auto iNode : bnd2node.father->operator[](iBnd))
                    node2bndNodeGlobal[iNode]++; // now stores num-reference
 
            {
                auto node2bndNodeGlobalPast = make_ssp<tInd::element_type>(ObjName{"node2bndNodeGlobalPast"}, mpi);
                DNDS::ArrayTransformerType<tInd::element_type>::Type node2bndNodeGlobalPastTrans;
                node2bndNodeGlobalPastTrans.setFatherSon(node2bndNodeGlobal.son, node2bndNodeGlobalPast);
                node2bndNodeGlobalPastTrans.createFatherGlobalMapping();
                std::vector<index> pushSonSeries(node2bndNodeGlobal.son->Size());
                for (size_t i = 0; i < pushSonSeries.size(); i++)
                    pushSonSeries[i] = DNDS::size_to_index(i);
                node2bndNodeGlobalPastTrans.createGhostMapping(pushSonSeries, node2bndNodeGlobal.trans.pLGhostMapping->ghostStart);
                node2bndNodeGlobalPastTrans.createMPITypes();
                node2bndNodeGlobalPastTrans.pullOnce();
                DNDS_assert(DNDS::size_to_index(node2bndNodeGlobal.trans.pLGhostMapping->ghostIndex.size()) == node2bndNodeGlobal.son->Size());
                DNDS_assert(DNDS::size_to_index(node2bndNodeGlobal.trans.pLGhostMapping->pushingIndexGlobal.size()) == node2bndNodeGlobalPast->Size());
                for (index i = 0; i < node2bndNodeGlobalPast->Size(); i++)
                {
                    index iNodeG = node2bndNodeGlobal.trans.pLGhostMapping->pushingIndexGlobal[i]; //?should be right
                    auto [ret, rank, iNodeL] = node2bndNodeGlobal.trans.pLGlobalMapping->search(iNodeG);
                    DNDS_assert(ret && rank == node2bndNodeGlobal.father->getMPI().rank); // has to be local main
                    node2bndNodeGlobal[iNodeL] += (*node2bndNodeGlobalPast)[i];
                }
            }
            {
                for (index i = 0; i < node2bndNodeGlobal.father->Size(); i++)
                    if (node2bndNodeGlobal[i] == 0)
                        node2bndNodeGlobal[i] = UnInitIndex;
                    else
                        node2bndNodeGlobal[i] = bndNodeCount++;
            }
            MPI::Scan(&bndNodeCount, &bndNodeStart, 1, DNDS_MPI_INDEX, MPI_SUM, mpi.comm);
            bndNodeStart -= bndNodeCount;
            DNDS_assert(bndNodeStart >= 0 && bndNodeCount >= 0);
            for (index i = 0; i < node2bndNodeGlobal.Size(); i++)
                if (node2bndNodeGlobal[i] >= 0)
                    node2bndNodeGlobal[i] += bndNodeStart;
        }
        node2bndNodeGlobal.trans.pullOnce();
        bMesh.coords.father->Resize(bndNodeCount);
        bMesh.coords.TransAttach();
        bMesh.coords.trans.createFatherGlobalMapping();
        std::vector<index> bndNodePullingSet;
        for (index iNode = 0; iNode < this->NumNodeProc(); iNode++)
        {
            if (node2bndNodeGlobal[iNode] < 0)
                continue;
            // then node2bndNodeGlobal[iNode] >= 0, which covers all bnd2node entries
            auto [ret, rank, iNodeL_in_bnd] = bMesh.coords.trans.pLGlobalMapping->search(node2bndNodeGlobal[iNode]);
            DNDS_assert(ret);
            if (rank != node2bndNodeGlobal.father->getMPI().rank)
                bndNodePullingSet.push_back(node2bndNodeGlobal[iNode]); // bndNodePullingSet now also covers bnd2node needs
        }
        bMesh.coords.trans.createGhostMapping(bndNodePullingSet);
        bMesh.coords.trans.createMPITypes();
        bMesh.node2nodeOrig.father->Resize(bndNodeCount);
        bMesh.node2nodeOrig.TransAttach();
        bMesh.node2nodeOrig.trans.BorrowGGIndexing(bMesh.coords.trans);
        bMesh.node2nodeOrig.trans.createMPITypes();
 
        // std::cout << mpi.rank << ", " << bndNodePullingSet.size() << ", " << bndNodeCount << std::endl;
        node2bndNode.resize(this->NumNodeProc(), -1);
        bMesh.node2parentNode.resize(bMesh.coords.Size());
        for (index iN = 0; iN < node2bndNodeGlobal.Size(); iN++)
            node2bndNode.at(iN) = bMesh.NodeIndexGlobal2Local(node2bndNodeGlobal[iN]),
            DNDS_assert(node2bndNodeGlobal[iN] != UnInitIndex ? (node2bndNode.at(iN) >= 0) : true);
        for (size_t iNode = 0; iNode < node2bndNode.size(); iNode++)
            if (node2bndNode[iNode] >= 0)
                bMesh.node2parentNode.at(node2bndNode[iNode]) = iNode;
 
        // std::cout << mpi.rank << ", " << bndNodeCount << std::endl;
        for (index iBNode = 0; iBNode < bMesh.coords.Size(); iBNode++)
            bMesh.coords[iBNode] = coords[bMesh.node2parentNode[iBNode]];
        bMesh.coords.trans.pullOnce(); // excessive, should be identical
        for (index iBNode = 0; iBNode < bMesh.node2nodeOrig.Size(); iBNode++)
            bMesh.node2nodeOrig[iBNode] = node2nodeOrig[bMesh.node2parentNode[iBNode]];
 
        index nBndCellUse{0};
        for (index iB = 0; iB < this->NumBnd(); iB++)
            if (!FaceIDIsPeriodic(this->bndElemInfo(iB, 0).zone))
                nBndCellUse++;
        bMesh.cell2node.father->Resize(nBndCellUse);
        if (isPeriodic)
            bMesh.cell2nodePbi.father->Resize(nBndCellUse);
        bMesh.cell2cellOrig.father->Resize(nBndCellUse);
        bMesh.cellElemInfo.father->Resize(nBndCellUse);
        bMesh.cell2parentCell.resize(nBndCellUse, -1);
        nBndCellUse = 0;
        for (index iB = 0; iB < this->NumBnd(); iB++)
        {
            if (FaceIDIsPeriodic(this->bndElemInfo(iB, 0).zone))
                continue;
            bMesh.cell2parentCell.at(nBndCellUse) = iB;
            bMesh.cell2node.ResizeRow(nBndCellUse, bnd2node.RowSize(iB));
            if (isPeriodic)
                bMesh.cell2nodePbi.ResizeRow(nBndCellUse, bnd2node.RowSize(iB));
            bMesh.cellElemInfo(nBndCellUse, 0) = bndElemInfo(iB, 0);
            bMesh.cell2cellOrig(nBndCellUse, 0) = bnd2bndOrig(iB, 0);
 
            for (rowsize ib2n = 0; ib2n < bnd2node.RowSize(iB); ib2n++)
            {
                if (bnd2faceV.at(iB) < 0) // where bnd has not a face!
                    bMesh.cell2node[nBndCellUse][ib2n] = node2bndNode.at(bnd2node[iB][ib2n]);
                else
                    bMesh.cell2node[nBndCellUse][ib2n] = node2bndNode.at(face2node[bnd2faceV.at(iB)][ib2n]); //* respect the face ordering if possible // this can be omitted if all bnds used are not periodic
                DNDS_assert(node2bndNode.at(bnd2node[iB][ib2n]) >= 0);
                if (isPeriodic)
                {
                    if (bnd2faceV.at(iB) < 0)                                             // where bnd has not a face!
                        bMesh.cell2nodePbi[nBndCellUse][ib2n] = Geom::NodePeriodicBits{}; // a invalid value
                    else
                    {
                        bMesh.cell2nodePbi[nBndCellUse][ib2n] = face2nodePbi[bnd2faceV.at(iB)][ib2n];
                    }
                }
            }
            nBndCellUse++;
        }
 
        bMesh.cell2node.father->Compress();
        bMesh.cell2node.father->createGlobalMapping();
        bMesh.cell2node.TransAttach();
        bMesh.cell2node.trans.createGhostMapping(std::vector<int>{}); // now bnd mesh has no valid cell2cell and ghost cells
 
        // Ensure global mappings exist on all arrays that BuildSerialOut
        // (and other consumers) may query via globalSize().  The ghost
        // mappings were already set up for coords/node2nodeOrig above;
        // cellElemInfo and cell2cellOrig only need the global mapping.
        bMesh.cellElemInfo.father->createGlobalMapping();
        bMesh.cell2cellOrig.father->createGlobalMapping();
 
        bMesh.adjPrimaryState = Adj_PointToLocal;
        if (mpi.rank == mRank)
            log() << "UnstructuredMesh === ConstructBndMesh Done" << std::endl;
    }
 
    void UnstructuredMesh::ObtainSymmetricSymbolicFactorization(Direct::SerialSymLUStructure &symLU, Direct::DirectPrecControl control) const
    {
        if (control.useDirectPrec)
            symLU.ObtainSymmetricSymbolicFactorization(
                cell2cellFaceVLocalParts,
                this->localPartitionStarts,
                control.getILUCode());
    }
 
    // =================================================================
    // File-local helpers for ReorderLocalCells
    // =================================================================
    namespace
    {
        /**
         * \brief Result of local cell permutation computation.
         */
        struct CellPermutationResult
        {
            std::vector<index> cellOld2New;
            std::vector<index> cellNew2Old;
            std::vector<index> localPartitionStarts;
            index bwOld = 0;
            index bwNew = 0;
        };
 
        /**
         * \brief Compute a cell reordering permutation using Metis partitioning
         *        with optional inner partitioning and contiguous sorting.
         *
         * 1. Partition via Metis + RCM.
         * 2. Optionally sub-partition each first-level partition.
         * 3. Within each partition, sort cells so interior (private) cells come
         *    before cells that touch ghost neighbors.
         * 4. Build inverse permutation.
         *
         * \param[in]  cell2cellFaceV  Local face-adjacency graph (no ghost edges).
         * \param[in]  cell2cell       Full cell-to-cell adjacency (with ghost).
         * \param[in]  nCell           Number of local (father) cells.
         * \param[in]  nParts          Number of first-level partitions.
         * \param[in]  nPartsInner     Number of inner partitions per first-level part.
         * \param[in]  localPartStartsIn  Current partition starts (for contiguous sort).
         */
        CellPermutationResult ComputeCellPermutation(
            tLocalMatStruct &cell2cellFaceV,
            const tAdjPair &cell2cell,
            index nCell,
            int nParts,
            int nPartsInner)
        {
            CellPermutationResult result;
            result.cellOld2New.resize(nCell, -1);
            result.cellNew2Old.resize(nCell);
            for (index i = 0; i < nCell; i++)
                result.cellNew2Old[i] = i;
 
            result.localPartitionStarts = ReorderSerialAdj_PartitionMetisC(
                cell2cellFaceV.begin(),
                cell2cellFaceV.end(),
                result.cellNew2Old.begin(),
                result.cellNew2Old.end(), nParts, 0, nPartsInner <= 1, result.bwOld, result.bwNew);
 
            if (nPartsInner > 1)
            {
                //! Debug assertions: verify graph invariants around inner partitioning.
                auto dbgCheckSubGraphRanges = [&](const char *tag)
                {
                    for (int p = 0; p < static_cast<int>(result.localPartitionStarts.size()) - 1; p++)
                    {
                        index pStart = result.localPartitionStarts[p];
                        index pEnd = result.localPartitionStarts[p + 1];
                        for (index iC = pStart; iC < pEnd; iC++)
                            for (auto jC : cell2cellFaceV[iC])
                                DNDS_assert_infof(
                                    jC >= 0 && jC < nCell,
                                    "%s: partition %d [%lld,%lld): cell %lld has neighbor %lld outside [0,%lld)",
                                    tag, p, (long long)pStart, (long long)pEnd,
                                    (long long)iC, (long long)jC, (long long)nCell);
                    }
                };
                auto dbgCheckBidir = [&](const char *tag)
                {
                    for (index iC = 0; iC < nCell; iC++)
                        for (auto jC : cell2cellFaceV[iC])
                        {
                            bool found = false;
                            for (auto kC : cell2cellFaceV[jC])
                                if (kC == iC)
                                {
                                    found = true;
                                    break;
                                }
                            DNDS_assert_infof(found,
                                              "%s: edge %lld->%lld exists but reverse %lld->%lld missing",
                                              tag, (long long)iC, (long long)jC, (long long)jC, (long long)iC);
                        }
                };
 
                dbgCheckSubGraphRanges("before inner partitioning");
                dbgCheckBidir("before inner partitioning");
 
                //! Each inner call partitions + RCM-reorders a first-level partition's sub-range,
                //! passing the full graph so cross-sub-graph references are updated in-place.
                for (int iPart = 0; iPart < static_cast<int>(result.localPartitionStarts.size()) - 1; iPart++)
                {
                    index bwOldC{0}, bwNewC{0};
                    index offset = result.localPartitionStarts[iPart];
                    index offsetN = result.localPartitionStarts[iPart + 1];
                    auto inner_parts_start = ReorderSerialAdj_PartitionMetisC(
                        cell2cellFaceV.begin() + offset,
                        cell2cellFaceV.begin() + offsetN,
                        result.cellNew2Old.begin() + offset,
                        result.cellNew2Old.begin() + offsetN, nPartsInner, offset, true, bwOldC, bwNewC,
                        cell2cellFaceV.begin(), nCell);
                    result.bwOld = std::max(result.bwOld, bwOldC);
                    result.bwNew = std::max(result.bwNew, bwNewC);
 
                    dbgCheckSubGraphRanges(fmt::format("after inner part {}", iPart).c_str());
                }
 
                dbgCheckBidir("after all inner partitioning");
            }
 
            // Contiguous sorting: within each partition, put interior cells before
            // cells that touch ghost neighbors.
            {
                auto cellIsNotPrivate = [&](index iCell)
                {
                    for (auto iCellOther : cell2cell[iCell])
                    {
                        if (iCellOther >= nCell)
                            return 1;
                    }
                    return 0;
                };
                int nLocalParts = result.localPartitionStarts.size() ? static_cast<int>(result.localPartitionStarts.size()) - 1 : 1;
                auto localPartStart = [&](int iPart) -> index
                { return result.localPartitionStarts.size() ? result.localPartitionStarts.at(iPart) : 0; };
                auto localPartEnd = [&](int iPart) -> index
                { return result.localPartitionStarts.size() ? result.localPartitionStarts.at(iPart + 1) : nCell; };
 
                std::vector<index> cellNew2Old_new;
                cellNew2Old_new.reserve(nCell);
                for (index i = 0; i < nCell; i++)
                    result.cellOld2New[i] /*tmp storage*/ = cellIsNotPrivate(result.cellNew2Old[i]);
                for (int iPart = 0; iPart < nLocalParts; iPart++) // we have to keep the local partitions intact
                {
                    for (index i = localPartStart(iPart); i < localPartEnd(iPart); i++)
                        if (!result.cellOld2New[i])
                            cellNew2Old_new.push_back(result.cellNew2Old[i]);
                    for (index i = localPartStart(iPart); i < localPartEnd(iPart); i++)
                        if (result.cellOld2New[i])
                            cellNew2Old_new.push_back(result.cellNew2Old[i]);
                }
 
                result.cellNew2Old = std::move(cellNew2Old_new);
                DNDS_assert(static_cast<index>(result.cellNew2Old.size()) == nCell);
                for (auto v : result.cellNew2Old)
                    DNDS_assert(v < nCell && v >= 0);
            }
 
            // Build inverse permutation
            std::unordered_set<index> set;
            set.reserve(result.cellNew2Old.size());
            for (index i = 0; i < nCell; i++)
            {
                DNDS_assert(set.count(result.cellNew2Old[i]) == 0);
                set.insert(result.cellNew2Old[i]);
                result.cellOld2New.at(result.cellNew2Old[i]) = i;
            }
            return result;
        }
    } // anonymous namespace
 
    void UnstructuredMesh::ReorderLocalCells(int nParts, int nPartsInner)
    {
        DNDS_assert(this->adjPrimaryState == Adj_PointToLocal);
        DNDS_check_throw(int64_t(nParts) * int64_t(nPartsInner) < std::numeric_limits<int>::max());
        nParts = std::max(nParts, 1);
        nPartsInner = std::max(nPartsInner, 1);
 
        // Section A: Compute cell permutation (Metis partition + contiguous sort)
        auto cell2cellFaceV = this->GetCell2CellFaceVLocal();
        auto perm = ComputeCellPermutation(
            cell2cellFaceV, cell2cell, NumCell(), nParts, nPartsInner);
        this->localPartitionStarts = std::move(perm.localPartitionStarts);
 
        MPI::AllreduceOneIndex(perm.bwOld, MPI_MAX, mpi);
        MPI::AllreduceOneIndex(perm.bwNew, MPI_MAX, mpi);
        if (mpi.rank == mRank)
            log()
                << fmt::format("UnstructuredMesh === ReorderLocalCells, nPart0 [{}], got reordering, bw [{}] to [{}]", nParts, perm.bwOld, perm.bwNew) << std::endl;
 
        // Section B: Build cellOld2NewArr for MPI communication of new cell indices
        tAdj1Pair cellOld2NewArr;
        cellOld2NewArr.InitPair("cellOld2NewArr", mpi);
        cellOld2NewArr.father->Resize(this->NumCell());
        for (index iCell = 0; iCell < this->NumCell(); iCell++)
            (*cellOld2NewArr.father)[iCell][0] = this->CellIndexLocal2Global_NoSon(perm.cellOld2New.at(iCell)); // this is right but bad syntax
        cellOld2NewArr.TransAttach();
 
        std::set<index> cellsNonLocalFull;
 
        auto record_iCellNonLocal = [&](index iCell)
        {
            if (iCell >= 0 && iCell < cell2node.father->pLGlobalMapping->globalSize())
            {
                if (std::get<1>(cell2node.father->pLGlobalMapping->search(iCell)) != mpi.rank)
                    cellsNonLocalFull.insert(iCell);
            }
            else
                DNDS_assert(iCell == UnInitIndex);
        };
 
        if (this->adjFacialState != Adj_Unknown)
        {
            DNDS_assert(this->adjFacialState == Adj_PointToLocal);
            this->AdjLocal2GlobalFacial();
            // see face2cell
            for (index iF = 0; iF < this->NumFace(); iF++)
                for (rowsize if2c = 0; if2c < 2; if2c++)
                    record_iCellNonLocal(face2cell(iF, if2c));
        }
        if (this->adjC2FState != Adj_Unknown)
        {
            DNDS_assert(this->adjC2FState == Adj_PointToLocal);
            this->AdjLocal2GlobalC2F();
        }
        if (this->adjN2CBState != Adj_Unknown)
        {
            DNDS_assert(this->adjN2CBState == Adj_PointToLocal);
            this->AdjLocal2GlobalN2CB();
            // see node2cell
            for (index iN = 0; iN < this->NumNode(); iN++)
                for (rowsize in2c = 0; in2c < node2cell[iN].size(); in2c++)
                    record_iCellNonLocal(node2cell(iN, in2c));
        }
        this->AdjLocal2GlobalPrimary();
        // see cell2cell
        for (index iC = 0; iC < this->NumCell(); iC++)
            for (rowsize ic2c = 0; ic2c < cell2cell[iC].size(); ic2c++)
                record_iCellNonLocal(cell2cell(iC, ic2c));
        // see bnd2cell
        for (index iB = 0; iB < this->NumBnd(); iB++)
            for (rowsize ib2c = 0; ib2c < 2; ib2c++)
                record_iCellNonLocal(bnd2cell(iB, ib2c));
 
        for (auto iCellGhost : cell2node.trans.pLGhostMapping->ghostIndex)
            cellsNonLocalFull.insert(iCellGhost); // ! should have no effect normally
 
        cellOld2NewArr.trans.createFatherGlobalMapping();
        cellOld2NewArr.trans.createGhostMapping(std::vector<index>(cellsNonLocalFull.begin(), cellsNonLocalFull.end()));
        cellOld2NewArr.trans.createMPITypes();
        cellOld2NewArr.trans.pullOnce();
 
        // Section C: Replace cell indices in RHS of xxx2cell arrays
        auto replaceCellIndexToNew = [&](index iCellOld) -> index
        {
            if (iCellOld == UnInitIndex)
                return UnInitIndex;
            auto [ret, rank, val] = cellOld2NewArr.trans.pLGhostMapping->search_indexAppend(iCellOld);
            DNDS_assert(ret);
            return cellOld2NewArr(val, 0);
        };
 
        if (this->adjFacialState != Adj_Unknown)
            for (index iFace = 0; iFace < this->NumFace(); iFace++)
                for (index &iCell : face2cell[iFace])
                    iCell = replaceCellIndexToNew(iCell);
        if (this->adjN2CBState != Adj_Unknown)
            for (index iNode = 0; iNode < this->NumNode(); iNode++)
                for (index &iCell : node2cell[iNode])
                    iCell = replaceCellIndexToNew(iCell);
        for (index iCell = 0; iCell < this->NumCell(); iCell++)
            for (index &iCellOther : cell2cell[iCell])
                iCellOther = replaceCellIndexToNew(iCellOther);
        for (index iBnd = 0; iBnd < this->NumBnd(); iBnd++)
            for (index &iCell : bnd2cell[iBnd])
                iCell = replaceCellIndexToNew(iCell);
 
        // Section D: Pull ghost data for xxx2cell
        if (this->adjFacialState != Adj_Unknown)
            face2cell.trans.pullOnce();
        if (this->adjN2CBState != Adj_Unknown)
            node2cell.trans.pullOnce();
        cell2cell.trans.pullOnce(); // should be unnecessary
        bnd2cell.trans.pullOnce();
 
        // Section E: Permute LHS of cell2xxx arrays
        auto cellOld2NewLocal = [&](index iCell) -> index
        { return this->CellIndexGlobal2Local_NoSon(cellOld2NewArr(iCell, 0)); };
 
        PermuteRows(cell2node, this->NumCell(), cellOld2NewLocal);
        PermuteRows(cell2cell, this->NumCell(), cellOld2NewLocal);
        PermuteRows(cell2cellOrig, this->NumCell(), cellOld2NewLocal);
        if (this->adjC2FState != Adj_Unknown)
            PermuteRows(cell2face, this->NumCell(), cellOld2NewLocal);
        if (this->isPeriodic)
            PermuteRows(cell2nodePbi, this->NumCell(), cellOld2NewLocal);
        PermuteRows(cellElemInfo, this->NumCell(), cellOld2NewLocal);
 
        // Section F: Rebuild ghost mappings with new cell indices
        { // cell2node
            std::vector<index> ghostCellGlobalsNew = cell2node.trans.pLGhostMapping->ghostIndex;
            for (index &iCell : ghostCellGlobalsNew)
                iCell = replaceCellIndexToNew(iCell);
            cell2node.trans.createGhostMapping(ghostCellGlobalsNew);
            cell2node.trans.createMPITypes();
            cell2node.trans.pullOnce();
        }
        { // cell2cell
            cell2cell.trans.BorrowGGIndexing(cell2node.trans);
            cell2cell.trans.createMPITypes();
            cell2cell.trans.pullOnce();
        }
        { // cell2cellOrig
            cell2cellOrig.trans.BorrowGGIndexing(cell2node.trans);
            cell2cellOrig.trans.createMPITypes();
            cell2cellOrig.trans.pullOnce();
        }
        if (this->adjC2FState != Adj_Unknown)
        { // cell2face
            cell2face.trans.BorrowGGIndexing(cell2node.trans);
            cell2face.trans.createMPITypes();
            cell2face.trans.pullOnce();
        }
        if (this->isPeriodic)
        { // cell2nodePbi
            cell2nodePbi.trans.BorrowGGIndexing(cell2node.trans);
            cell2nodePbi.trans.createMPITypes();
            cell2nodePbi.trans.pullOnce();
        }
        { // cellElemInfo
            cellElemInfo.trans.BorrowGGIndexing(cell2node.trans);
            cellElemInfo.trans.createMPITypes();
            cellElemInfo.trans.pullOnce();
        }
 
        // Section G: Restore all adjacencies to local indices
        if (this->adjFacialState != Adj_Unknown)
        {
            this->AdjGlobal2LocalFacial();
        }
        if (this->adjC2FState != Adj_Unknown)
        {
            this->AdjGlobal2LocalC2F();
        }
        if (this->adjN2CBState != Adj_Unknown)
        {
            this->AdjGlobal2LocalN2CB();
        }
        this->AdjGlobal2LocalPrimary();
        if (mpi.rank == mRank)
            log() << fmt::format("UnstructuredMesh === ReorderLocalCells finished") << std::endl;
    }
}