CoupledLocalToGlobalC0ContMap.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: CoupledLocalToGlobalC0ContMap.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Wrapper class around the library
// LocalToGlobalC0ContMap class for use in the Couplied Linearised NS
// solver.
//
///////////////////////////////////////////////////////////////////////////////
 
#include <IncNavierStokesSolver/EquationSystems/CoupledLocalToGlobalC0ContMap.h>
#include <LocalRegions/Expansion1D.h>
#include <LocalRegions/Expansion2D.h>
#include <LocalRegions/SegExp.h>
#include <MultiRegions/GlobalLinSysDirectStaticCond.h>
#include <SpatialDomains/MeshGraph.h>
 
namespace Nektar
{
/**
 * This is an vector extension of
 * MultiRegions::AssemblyMapCG::SetUp2DExpansionC0ContMap related to the
 * Linearised Navier Stokes problem
 */
CoupledLocalToGlobalC0ContMap::CoupledLocalToGlobalC0ContMap(
    const LibUtilities::SessionReaderSharedPtr &pSession,
    [[maybe_unused]] const SpatialDomains::MeshGraphSharedPtr &graph,
    [[maybe_unused]] const SpatialDomains::BoundaryConditionsSharedPtr
        &boundaryConditions,
    const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
    const MultiRegions::ExpListSharedPtr &pressure, const size_t nz_loc,
    const bool CheckforSingularSys)
    : AssemblyMapCG(pSession, fields[0]->GetComm())
{
    size_t i, j, k, n;
    size_t cnt = 0, offset = 0;
    size_t meshVertId;
    size_t meshEdgeId;
    size_t bndEdgeCnt;
    size_t globalId;
    size_t nEdgeCoeffs;
    size_t nEdgeInteriorCoeffs;
    int firstNonDirGraphVertId;
    size_t nLocBndCondDofs    = 0;
    size_t nLocDirBndCondDofs = 0;
    int nExtraDirichlet       = 0;
    LocalRegions::Expansion2DSharedPtr locExpansion;
    LocalRegions::SegExpSharedPtr bndSegExp;
    LibUtilities::BasisType bType;
    StdRegions::Orientation edgeOrient;
    Array<OneD, unsigned int> edgeInteriorMap;
    Array<OneD, int> edgeInteriorSign;
    size_t nvel = fields.size();
 
    const LocalRegions::ExpansionVector &locExpVector = *(fields[0]->GetExp());
    int id, diff;
    size_t nel = fields[0]->GetNumElmts();
 
    MultiRegions::PeriodicMap periodicVerts;
    MultiRegions::PeriodicMap periodicEdges;
    MultiRegions::PeriodicMap periodicFaces;
    std::vector<std::map<int, int>> ReorderedGraphVertId(3);
    MultiRegions::BottomUpSubStructuredGraphSharedPtr bottomUpGraph;
    int staticCondLevel = 0;
 
    if (CheckforSingularSys) // all singularity checking by setting flag to true
    {
        m_systemSingular = true;
    }
    else // Turn off singular checking by setting flag to false
    {
        m_systemSingular = false;
    }
 
    /**
     * STEP 1: Wrap boundary conditions vector in an array
     * (since routine is set up for multiple fields) and call
     * the graph re-odering subroutine to obtain the reordered
     * values
     */
 
    // Obtain any periodic information and allocate default mapping array
    fields[0]->GetPeriodicEntities(periodicVerts, periodicEdges, periodicFaces);
 
    const Array<OneD, const MultiRegions::ExpListSharedPtr> bndCondExp =
        fields[0]->GetBndCondExpansions();
 
    Array<OneD, Array<OneD, const SpatialDomains::BoundaryConditionShPtr>>
        bndConditionsVec(nvel);
    for (i = 0; i < nvel; ++i)
    {
        bndConditionsVec[i] = fields[i]->GetBndConditions();
    }
 
    std::map<int, int> IsDirVertDof;
    std::map<int, int> IsDirEdgeDof;
 
    SpatialDomains::Geometry1D *g;
    for (j = 0; j < bndCondExp.size(); ++j)
    {
        std::map<int, int> BndExpVids;
        // collect unique list of vertex ids for this expansion
        for (k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
        {
            g = bndCondExp[j]
                    ->GetExp(k)
                    ->as<LocalRegions::Expansion1D>()
                    ->GetGeom1D();
            BndExpVids[g->GetVid(0)] = g->GetVid(0);
            BndExpVids[g->GetVid(1)] = g->GetVid(1);
        }
 
        for (i = 0; i < nvel; ++i)
        {
            if (bndConditionsVec[i][j]->GetBoundaryConditionType() ==
                SpatialDomains::eDirichlet)
            {
                // set number of Dirichlet conditions along edge
                for (k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
                {
                    IsDirEdgeDof[bndCondExp[j]
                                     ->GetExp(k)
                                     ->as<LocalRegions::Expansion1D>()
                                     ->GetGeom1D()
                                     ->GetGlobalID()] += 1;
                }
 
                // Set number of Dirichlet conditions at vertices
                // with a clamp on its maximum value being nvel to
                // handle corners between expansions
                for (auto &mapIt : BndExpVids)
                {
                    id                         = IsDirVertDof[mapIt.second] + 1;
                    IsDirVertDof[mapIt.second] = (id > (int)nvel) ? nvel : id;
                }
            }
            else
            {
                // Check to see that edge normals have non-zero
                // component in this direction since otherwise
                // also can be singular.
                /// @TODO: Fix this so that we can extract normals from edges
                for (k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
                {
                    Array<OneD, Array<OneD, NekDouble>> locnorm;
                    LocalRegions::Expansion1DSharedPtr loc_exp =
                        bndCondExp[j]
                            ->GetExp(k)
                            ->as<LocalRegions::Expansion1D>();
                    locnorm =
                        loc_exp->GetLeftAdjacentElementExp()->GetTraceNormal(
                            loc_exp->GetLeftAdjacentElementTrace());
 
                    size_t ndir = locnorm.size();
                    if (i < ndir) // account for Fourier version where n can be
                                  // larger then ndir
                    {
                        for (size_t l = 0; l < locnorm[0].size(); ++l)
                        {
                            if (fabs(locnorm[i][l]) > NekConstants::kNekZeroTol)
                            {
                                m_systemSingular = false;
                                break;
                            }
                        }
                    }
                    if (m_systemSingular == false)
                    {
                        break;
                    }
                }
            }
        }
    }
 
    Array<OneD, std::map<int, int>> Dofs(2);
 
    Array<OneD, int> AddMeanPressureToEdgeId(nel, -1);
    size_t edgeId, vertId;
 
    // special case of singular problem - need to fix one pressure
    // dof to a dirichlet edge. Since we attached pressure dof to
    // last velocity component of edge need to make sure this
    // component is Dirichlet
    if (m_systemSingular)
    {
        id = -1;
        for (i = 0; i < bndConditionsVec[0].size(); ++i)
        {
            if (bndConditionsVec[nvel - 1][i]->GetBoundaryConditionType() ==
                SpatialDomains::eDirichlet)
            {
                id = bndCondExp[i]
                         ->GetExp(0)
                         ->as<LocalRegions::Expansion1D>()
                         ->GetGeom1D()
                         ->GetGlobalID();
                break;
            }
        }
 
        ASSERTL0(id != -1, " Did not find an edge to attach singular pressure "
                           "degree of freedom");
 
        // determine element with this edge id. There may be a
        // more direct way of getting element from spatialDomains
        for (i = 0; i < nel; ++i)
        {
            for (j = 0; j < (size_t)locExpVector[i]->GetNverts(); ++j)
            {
                edgeId = (locExpVector[i]
                              ->as<LocalRegions::Expansion2D>()
                              ->GetGeom2D())
                             ->GetEid(j);
 
                if ((int)edgeId == id)
                {
                    AddMeanPressureToEdgeId[i] = id;
                    break;
                }
            }
 
            if (AddMeanPressureToEdgeId[i] != -1)
            {
                break;
            }
        }
    }
 
    for (i = 0; i < nel; ++i)
    {
        for (j = 0; j < (size_t)locExpVector[i]->GetNverts(); ++j)
        {
            vertId =
                (locExpVector[i]->as<LocalRegions::Expansion2D>()->GetGeom2D())
                    ->GetVid(j);
            if (Dofs[0].count(vertId) == 0)
            {
                Dofs[0][vertId] = nvel * nz_loc;
 
                // Adjust for a Dirichlet boundary condition to give number to
                // be solved
                if (IsDirVertDof.count(vertId) != 0)
                {
                    Dofs[0][vertId] -= IsDirVertDof[vertId] * nz_loc;
                }
            }
 
            edgeId =
                (locExpVector[i]->as<LocalRegions::Expansion2D>()->GetGeom2D())
                    ->GetEid(j);
            if (Dofs[1].count(edgeId) == 0)
            {
                Dofs[1][edgeId] =
                    nvel * (locExpVector[i]->GetTraceNcoeffs(j) - 2) * nz_loc;
            }
 
            // Adjust for Dirichlet boundary conditions to give number to be
            // solved
            if (IsDirEdgeDof.count(edgeId) != 0)
            {
                Dofs[1][edgeId] -= IsDirEdgeDof[edgeId] * nz_loc *
                                   (locExpVector[i]->GetTraceNcoeffs(j) - 2);
            }
        }
    }
 
    std::set<int> extraDirVerts, extraDirEdges;
 
    CreateGraph(*fields[0], bndCondExp, bndConditionsVec, false, periodicVerts,
                periodicEdges, periodicFaces, ReorderedGraphVertId,
                bottomUpGraph, extraDirVerts, extraDirEdges,
                firstNonDirGraphVertId, nExtraDirichlet, 4);
    /*
    SetUp2DGraphC0ContMap(*fields[0],
                          bndCondExp,
                          bndConditionsVec,
                          periodicVerts,          periodicEdges,
                          Dofs,                   ReorderedGraphVertId,
                          firstNonDirGraphVertId, nExtraDirichlet,
                          bottomUpGraph, extraDir,  false,  4);
    */
 
    /**
     * STEP 2a: Set the mean pressure modes to edges depending on
     * type of direct solver technique;
     */
 
    // determine which edge to add mean pressure dof based on
    // ensuring that at least one pressure dof from an internal
    // patch is associated with its boundary system
    if (m_session->MatchSolverInfoAsEnum(
            "GlobalSysSoln", MultiRegions::eDirectMultiLevelStaticCond))
    {
 
        FindEdgeIdToAddMeanPressure(ReorderedGraphVertId, nel, locExpVector,
                                    edgeId, vertId, firstNonDirGraphVertId,
                                    IsDirEdgeDof, bottomUpGraph,
                                    AddMeanPressureToEdgeId);
    }
 
    // Set unset elmts to non-Dirichlet edges.
    // special case of singular problem - need to fix one
    // pressure dof to a dirichlet edge
    for (i = 0; i < nel; ++i)
    {
        for (j = 0; j < (size_t)locExpVector[i]->GetNverts(); ++j)
        {
            edgeId =
                (locExpVector[i]->as<LocalRegions::Expansion2D>()->GetGeom2D())
                    ->GetEid(j);
 
            if (IsDirEdgeDof.count(edgeId) == 0) // interior edge
            {
                // setup AddMeanPressureToEdgeId to decide where to
                // put pressure
                if (AddMeanPressureToEdgeId[i] == -1)
                {
                    AddMeanPressureToEdgeId[i] = edgeId;
                }
            }
        }
        ASSERTL0((AddMeanPressureToEdgeId[i] != -1),
                 "Did not determine "
                 "an edge to attach mean pressure dof");
        // Add the mean pressure degree of freedom to this edge
        Dofs[1][AddMeanPressureToEdgeId[i]] += nz_loc;
    }
 
    std::map<int, int> pressureEdgeOffset;
 
    /**
     * STEP 2: Count out the number of Dirichlet vertices and edges first
     */
    for (i = 0; i < bndCondExp.size(); i++)
    {
        for (j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
        {
            bndSegExp = bndCondExp[i]->GetExp(j)->as<LocalRegions::SegExp>();
            for (k = 0; k < nvel; ++k)
            {
                if (bndConditionsVec[k][i]->GetBoundaryConditionType() ==
                    SpatialDomains::eDirichlet)
                {
                    nLocDirBndCondDofs += bndSegExp->GetNcoeffs() * nz_loc;
                }
 
                if (bndConditionsVec[k][i]->GetBoundaryConditionType() !=
                    SpatialDomains::ePeriodic)
                {
                    nLocBndCondDofs += bndSegExp->GetNcoeffs() * nz_loc;
                }
            }
        }
    }
 
    if (m_systemSingular)
    {
        m_numLocalDirBndCoeffs = nLocDirBndCondDofs + nExtraDirichlet + nz_loc;
    }
    else
    {
        m_numLocalDirBndCoeffs = nLocDirBndCondDofs + nExtraDirichlet;
    }
 
    /**
     * STEP 3: Set up an array which contains the offset information of
     * the different graph vertices.
     *
     * This basically means to identify how many global degrees of
     * freedom the individual graph vertices correspond. Obviously,
     * the graph vertices corresponding to the mesh-vertices account
     * for a single global DOF. However, the graph vertices
     * corresponding to the element edges correspond to 2*(N-2) global DOF
     * where N is equal to the number of boundary modes on this edge.
     */
    Array<OneD, int> graphVertOffset(
        nvel * nz_loc *
            (ReorderedGraphVertId[0].size() + ReorderedGraphVertId[1].size()),
        0);
    graphVertOffset[0] = 0;
 
    m_signChange = false;
 
    for (i = 0; i < nel; ++i)
    {
        locExpansion = locExpVector[i]->as<LocalRegions::Expansion2D>();
 
        for (j = 0; j < (size_t)locExpansion->GetNtraces(); ++j)
        {
            nEdgeCoeffs = locExpansion->GetTraceNcoeffs(j);
            meshEdgeId  = (locExpansion->GetGeom2D())->GetEid(j);
            meshVertId  = (locExpansion->GetGeom2D())->GetVid(j);
 
            for (k = 0; k < nvel * nz_loc; ++k)
            {
                graphVertOffset[ReorderedGraphVertId[0][meshVertId] * nvel *
                                    nz_loc +
                                k] = 1;
                graphVertOffset[ReorderedGraphVertId[1][meshEdgeId] * nvel *
                                    nz_loc +
                                k] = (nEdgeCoeffs - 2);
            }
 
            bType = locExpansion->GetBasisType(0);
            // need a sign vector for modal expansions if nEdgeCoeffs >=4
            if ((nEdgeCoeffs >= 4) && (bType == LibUtilities::eModified_A))
            {
                m_signChange = true;
            }
        }
    }
 
    // Add mean pressure modes;
    for (i = 0; i < nel; ++i)
    {
        graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]] +
                         1) *
                            nvel * nz_loc -
                        1] += nz_loc;
        // graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]])*nvel*nz_loc]
        // += nz_loc;
    }
 
    // Negate the vertices and edges with only a partial
    // Dirichlet conditon. Essentially we check to see if an edge
    // has a mixed Dirichlet with Neumann/Robin Condition and if
    // so negate the offset associated with this vertex.
 
    std::map<int, int> DirVertChk;
 
    for (i = 0; i < bndConditionsVec[0].size(); ++i)
    {
        cnt = 0;
        for (j = 0; j < nvel; ++j)
        {
            if (bndConditionsVec[j][i]->GetBoundaryConditionType() ==
                SpatialDomains::eDirichlet)
            {
                cnt++;
            }
        }
 
        // Case where partial Dirichlet boundary condition
        if ((cnt > 0) && (cnt < nvel))
        {
            for (j = 0; j < nvel; ++j)
            {
                if (bndConditionsVec[j][i]->GetBoundaryConditionType() ==
                    SpatialDomains::eDirichlet)
                {
                    // negate graph offsets which should be
                    // Dirichlet conditions
                    for (k = 0; k < bndCondExp[i]->GetNumElmts(); ++k)
                    {
                        // vertices with mix condition;
                        id = bndCondExp[i]
                                 ->GetExp(k)
                                 ->as<LocalRegions::Expansion1D>()
                                 ->GetGeom1D()
                                 ->GetVid(0);
                        if (DirVertChk.count(id * nvel + j) == 0)
                        {
                            DirVertChk[id * nvel + j] = 1;
                            for (n = 0; n < nz_loc; ++n)
                            {
                                graphVertOffset[ReorderedGraphVertId[0][id] *
                                                    nvel * nz_loc +
                                                j * nz_loc + n] *= -1;
                            }
                        }
 
                        id = bndCondExp[i]
                                 ->GetExp(k)
                                 ->as<LocalRegions::Expansion1D>()
                                 ->GetGeom1D()
                                 ->GetVid(1);
                        if (DirVertChk.count(id * nvel + j) == 0)
                        {
                            DirVertChk[id * nvel + j] = 1;
                            for (n = 0; n < nz_loc; ++n)
                            {
                                graphVertOffset[ReorderedGraphVertId[0][id] *
                                                    nvel * nz_loc +
                                                j * nz_loc + n] *= -1;
                            }
                        }
 
                        // edges with mixed id;
                        id = bndCondExp[i]
                                 ->GetExp(k)
                                 ->as<LocalRegions::Expansion1D>()
                                 ->GetGeom1D()
                                 ->GetGlobalID();
                        for (n = 0; n < nz_loc; ++n)
                        {
                            graphVertOffset[ReorderedGraphVertId[1][id] * nvel *
                                                nz_loc +
                                            j * nz_loc + n] *= -1;
                        }
                    }
                }
            }
        }
    }
 
    cnt = 0;
    // assemble accumulative list of full Dirichlet values.
    for (i = 0; i < firstNonDirGraphVertId * nvel * nz_loc; ++i)
    {
        diff               = abs(graphVertOffset[i]);
        graphVertOffset[i] = cnt;
        cnt += diff;
    }
 
    // set Dirichlet values with negative values to Dirichlet value
    for (i = firstNonDirGraphVertId * nvel * nz_loc; i < graphVertOffset.size();
         ++i)
    {
        if (graphVertOffset[i] < 0)
        {
            diff               = -graphVertOffset[i];
            graphVertOffset[i] = -cnt;
            cnt += diff;
        }
    }
 
    // Accumulate all interior degrees of freedom with positive values
    m_numGlobalDirBndCoeffs = cnt;
 
    // offset values
    for (i = firstNonDirGraphVertId * nvel * nz_loc; i < graphVertOffset.size();
         ++i)
    {
        if (graphVertOffset[i] >= 0)
        {
            diff               = graphVertOffset[i];
            graphVertOffset[i] = cnt;
            cnt += diff;
        }
    }
 
    // Finally set negative entries (corresponding to Dirichlet
    // values ) to be positive
    for (i = firstNonDirGraphVertId * nvel * nz_loc; i < graphVertOffset.size();
         ++i)
    {
        if (graphVertOffset[i] < 0)
        {
            graphVertOffset[i] = -graphVertOffset[i];
        }
    }
 
    // Allocate the proper amount of space for the class-data and fill
    // information that is already known
    cnt                 = 0;
    m_numLocalBndCoeffs = 0;
    m_numLocalCoeffs    = 0;
 
    for (i = 0; i < nel; ++i)
    {
        m_numLocalBndCoeffs +=
            nz_loc * (nvel * locExpVector[i]->NumBndryCoeffs() + 1);
        // add these coeffs up separately since
        // pressure->GetNcoeffs can include the coefficient in
        // multiple planes.
        m_numLocalCoeffs += (pressure->GetExp(i)->GetNcoeffs() - 1) * nz_loc;
    }
 
    m_numLocalCoeffs += m_numLocalBndCoeffs;
 
    m_localToGlobalMap         = Array<OneD, int>(m_numLocalCoeffs, -1);
    m_localToGlobalBndMap      = Array<OneD, int>(m_numLocalBndCoeffs, -1);
    m_bndCondIDToGlobalTraceID = Array<OneD, int>(nLocBndCondDofs, -1);
 
    // Set default sign array.
    m_localToGlobalSign    = Array<OneD, NekDouble>(m_numLocalCoeffs, 1.0);
    m_localToGlobalBndSign = Array<OneD, NekDouble>(m_numLocalBndCoeffs, 1.0);
 
    m_staticCondLevel = staticCondLevel;
    m_numPatches      = nel;
 
    m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(nel);
    m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(nel);
 
    for (i = 0; i < nel; ++i)
    {
        m_numLocalBndCoeffsPerPatch[i] =
            (unsigned int)nz_loc *
            (nvel * locExpVector[i]->NumBndryCoeffs() + 1);
        m_numLocalIntCoeffsPerPatch[i] =
            (unsigned int)nz_loc * (pressure->GetExp(i)->GetNcoeffs() - 1);
    }
 
    // Set up local to local bnd and local int maps
    m_localToLocalBndMap = Array<OneD, int>(m_numLocalBndCoeffs, -1);
    m_localToLocalIntMap =
        Array<OneD, int>(m_numLocalCoeffs - m_numLocalBndCoeffs, -1);
 
    size_t bndcnt = 0;
    size_t intcnt = 0;
    cnt           = 0;
    for (i = 0; i < nel; ++i)
    {
        for (j = 0; j < nz_loc * (nvel * locExpVector[i]->NumBndryCoeffs());
             ++j)
        {
            m_localToLocalBndMap[bndcnt++] = cnt++;
        }
 
        for (n = 0; n < nz_loc; ++n)
        {
            m_localToLocalBndMap[bndcnt++] = cnt++;
            for (j = 1; j < (size_t)pressure->GetExp(i)->GetNcoeffs(); ++j)
            {
                m_localToLocalIntMap[intcnt++] = cnt++;
            }
        }
    }
 
    /**
     * STEP 4: Now, all ingredients are ready to set up the actual
     * local to global mapping.
     *
     * The remainder of the map consists of the element-interior
     * degrees of freedom. This leads to the block-diagonal submatrix
     * as each element-interior mode is globally orthogonal to modes
     * in all other elements.
     */
    cnt = 0;
    size_t nv, velnbndry;
    Array<OneD, unsigned int> bmap;
 
    // Loop over all the elements in the domain in shuffled
    // ordering (element type consistency)
    for (i = 0; i < nel; ++i)
    {
        locExpansion = locExpVector[i]->as<LocalRegions::Expansion2D>();
 
        velnbndry = locExpansion->NumBndryCoeffs();
 
        // Require an inverse ordering of the bmap system to store
        // local numbering system. Therefore get hold of elemental
        // bmap and set up an inverse map
        std::map<int, int> inv_bmap;
        locExpansion->GetBoundaryMap(bmap);
        for (j = 0; j < bmap.size(); ++j)
        {
            inv_bmap[bmap[j]] = j;
        }
 
        // Loop over all edges (and vertices) of element i
        for (j = 0; j < (size_t)locExpansion->GetNtraces(); ++j)
        {
            nEdgeInteriorCoeffs = locExpansion->GetTraceNcoeffs(j) - 2;
            edgeOrient          = (locExpansion->GetGeom())->GetEorient(j);
            meshEdgeId          = (locExpansion->GetGeom())->GetEid(j);
            meshVertId          = (locExpansion->GetGeom())->GetVid(j);
 
            auto pIt = periodicEdges.find(meshEdgeId);
 
            // See if this edge is periodic. If it is, then we map all
            // edges to the one with lowest ID, and align all
            // coefficients to this edge orientation.
            if (pIt != periodicEdges.end())
            {
                std::pair<int, StdRegions::Orientation> idOrient =
                    DeterminePeriodicEdgeOrientId(meshEdgeId, edgeOrient,
                                                  pIt->second);
                edgeOrient = idOrient.second;
            }
 
            locExpansion->GetTraceInteriorToElementMap(
                j, edgeInteriorMap, edgeInteriorSign, edgeOrient);
 
            // Set the global DOF for vertex j of element i
            for (nv = 0; nv < nvel * nz_loc; ++nv)
            {
                m_localToGlobalMap[cnt + nv * velnbndry +
                                   inv_bmap[locExpansion->GetVertexMap(j)]] =
                    graphVertOffset[ReorderedGraphVertId[0][meshVertId] * nvel *
                                        nz_loc +
                                    nv];
 
                // Set the global DOF's for the interior modes of edge j
                for (k = 0; k < nEdgeInteriorCoeffs; ++k)
                {
                    m_localToGlobalMap[cnt + nv * velnbndry +
                                       inv_bmap[edgeInteriorMap[k]]] =
                        graphVertOffset[ReorderedGraphVertId[1][meshEdgeId] *
                                            nvel * nz_loc +
                                        nv] +
                        k;
                }
            }
 
            // Fill the sign vector if required
            if (m_signChange)
            {
                for (nv = 0; nv < nvel * nz_loc; ++nv)
                {
                    for (k = 0; k < nEdgeInteriorCoeffs; ++k)
                    {
                        m_localToGlobalSign[cnt + nv * velnbndry +
                                            inv_bmap[edgeInteriorMap[k]]] =
                            (NekDouble)edgeInteriorSign[k];
                    }
                }
            }
        }
 
        // use difference between two edges of the
        // AddMeanPressureEdgeId to det nEdgeInteriorCoeffs.
        nEdgeInteriorCoeffs =
            graphVertOffset[(ReorderedGraphVertId[1]
                                                 [AddMeanPressureToEdgeId[i]]) *
                                nvel * nz_loc +
                            1] -
            graphVertOffset[(ReorderedGraphVertId[1]
                                                 [AddMeanPressureToEdgeId[i]]) *
                            nvel * nz_loc];
 
        size_t psize = pressure->GetExp(i)->GetNcoeffs();
        for (n = 0; n < nz_loc; ++n)
        {
            m_localToGlobalMap[cnt + nz_loc * nvel * velnbndry + n * psize] =
                graphVertOffset
                    [(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]] + 1) *
                         nvel * nz_loc -
                     1] +
                nEdgeInteriorCoeffs +
                pressureEdgeOffset[AddMeanPressureToEdgeId[i]];
 
            pressureEdgeOffset[AddMeanPressureToEdgeId[i]] += 1;
        }
 
        cnt += (velnbndry * nvel + psize) * nz_loc;
    }
 
    // Set up the mapping for the boundary conditions
    offset = cnt = 0;
    for (nv = 0; nv < nvel; ++nv)
    {
        for (i = 0; i < bndCondExp.size(); i++)
        {
            if (bndConditionsVec[nv][i]->GetBoundaryConditionType() ==
                SpatialDomains::ePeriodic)
            {
                continue;
            }
 
            for (n = 0; n < nz_loc; ++n)
            {
                int ncoeffcnt = 0;
                for (j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
                {
                    bndSegExp =
                        bndCondExp[i]->GetExp(j)->as<LocalRegions::SegExp>();
 
                    cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
                    for (k = 0; k < 2; k++)
                    {
                        meshVertId = (bndSegExp->GetGeom1D())->GetVid(k);
                        m_bndCondIDToGlobalTraceID[cnt +
                                                   bndSegExp->GetVertexMap(k)] =
                            graphVertOffset[ReorderedGraphVertId[0]
                                                                [meshVertId] *
                                                nvel * nz_loc +
                                            nv * nz_loc + n];
                    }
 
                    meshEdgeId  = (bndSegExp->GetGeom1D())->GetGlobalID();
                    bndEdgeCnt  = 0;
                    nEdgeCoeffs = bndSegExp->GetNcoeffs();
                    for (k = 0; k < nEdgeCoeffs; k++)
                    {
                        if (m_bndCondIDToGlobalTraceID[cnt + k] == -1)
                        {
                            m_bndCondIDToGlobalTraceID[cnt + k] =
                                graphVertOffset
                                    [ReorderedGraphVertId[1][meshEdgeId] *
                                         nvel * nz_loc +
                                     nv * nz_loc + n] +
                                bndEdgeCnt;
                            bndEdgeCnt++;
                        }
                    }
                    ncoeffcnt += nEdgeCoeffs;
                }
                // Note: Can not use bndCondExp[i]->GetNcoeffs()
                // due to homogeneous extension not returning just
                // the value per plane
                offset += ncoeffcnt;
            }
        }
    }
 
    globalId = Vmath::Vmax(m_numLocalCoeffs, &m_localToGlobalMap[0], 1) + 1;
    m_numGlobalBndCoeffs = globalId;
 
    /**
     * STEP 5: The boundary condition mapping is generated from the
     * same vertex renumbering and fill in a unique interior map.
     */
    cnt = 0;
    for (i = 0; i < (size_t)m_numLocalCoeffs; ++i)
    {
        if (m_localToGlobalMap[i] == -1)
        {
            m_localToGlobalMap[i] = globalId++;
        }
        else
        {
            if (m_signChange)
            {
                m_localToGlobalBndSign[cnt] = m_localToGlobalSign[i];
            }
            m_localToGlobalBndMap[cnt++] = m_localToGlobalMap[i];
        }
    }
    m_numGlobalCoeffs = globalId;
 
    // Set up the local to global map for the next level when using
    // multi-level static condensation
    if (m_session->MatchSolverInfoAsEnum(
            "GlobalSysSoln", MultiRegions::eDirectMultiLevelStaticCond))
    {
        if (m_staticCondLevel < (bottomUpGraph->GetNlevels() - 1))
        {
            Array<OneD, int> vwgts_perm(Dofs[0].size() + Dofs[1].size() -
                                        firstNonDirGraphVertId);
            for (i = 0; i < locExpVector.size(); ++i)
            {
                locExpansion = locExpVector[i]->as<LocalRegions::Expansion2D>();
                size_t nVert = locExpansion->GetNverts();
                for (j = 0; j < nVert; ++j)
                {
                    meshEdgeId = (locExpansion->GetGeom2D())->GetEid(j);
                    meshVertId = (locExpansion->GetGeom2D())->GetVid(j);
 
                    if (ReorderedGraphVertId[0][meshVertId] >=
                        firstNonDirGraphVertId)
                    {
                        vwgts_perm[ReorderedGraphVertId[0][meshVertId] -
                                   firstNonDirGraphVertId] =
                            Dofs[0][meshVertId];
                    }
 
                    if (ReorderedGraphVertId[1][meshEdgeId] >=
                        firstNonDirGraphVertId)
                    {
                        vwgts_perm[ReorderedGraphVertId[1][meshEdgeId] -
                                   firstNonDirGraphVertId] =
                            Dofs[1][meshEdgeId];
                    }
                }
            }
 
            bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
 
            m_nextLevelLocalToGlobalMap =
                MemoryManager<AssemblyMap>::AllocateSharedPtr(this,
                                                              bottomUpGraph);
        }
    }
}
 
void CoupledLocalToGlobalC0ContMap::FindEdgeIdToAddMeanPressure(
    std::vector<std::map<int, int>> &ReorderedGraphVertId, size_t &nel,
    const LocalRegions::ExpansionVector &locExpVector, size_t &edgeId,
    [[maybe_unused]] size_t &vertId, int &firstNonDirGraphVertId,
    std::map<int, int> &IsDirEdgeDof,
    MultiRegions::BottomUpSubStructuredGraphSharedPtr &bottomUpGraph,
    Array<OneD, int> &AddMeanPressureToEdgeId)
{
    size_t i, j, k;
 
    // Make list of homogeneous graph edges to elmt mappings
    Array<TwoD, int> EdgeIdToElmts(ReorderedGraphVertId[1].size(), 2, -1);
    std::map<int, int> HomGraphEdgeIdToEdgeId;
 
    for (i = 0; i < nel; ++i)
    {
        size_t nVert = locExpVector[i]->GetNverts();
        for (j = 0; j < nVert; ++j)
        {
            edgeId =
                (locExpVector[i]->as<LocalRegions::Expansion2D>()->GetGeom2D())
                    ->GetEid(j);
 
            // note second condition stops us using mixed boundary condition
            if ((ReorderedGraphVertId[1][edgeId] >= firstNonDirGraphVertId) &&
                (IsDirEdgeDof.count(edgeId) == 0))
            {
                HomGraphEdgeIdToEdgeId[ReorderedGraphVertId[1][edgeId] -
                                       firstNonDirGraphVertId] = edgeId;
 
                if (EdgeIdToElmts[edgeId][0] == -1)
                {
                    EdgeIdToElmts[edgeId][0] = i;
                }
                else
                {
                    EdgeIdToElmts[edgeId][1] = i;
                }
            }
        }
    }
 
    // Start at second to last level and find edge on boundary
    // to attach element
    size_t nlevels = bottomUpGraph->GetNlevels();
 
    // determine a default edge to attach pressure modes to
    // which is part of the inner solve;
    int defedge = -1;
 
    std::vector<MultiRegions::SubGraphSharedPtr> bndgraphs =
        bottomUpGraph->GetInteriorBlocks(nlevels);
    for (i = 0; i < bndgraphs.size(); ++i)
    {
        int GlobIdOffset = bndgraphs[i]->GetIdOffset();
        size_t nVert     = bndgraphs[i]->GetNverts();
 
        for (j = 0; j < nVert; ++j)
        {
            // find edge in graph vert list
            if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
            {
                edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
 
                if (defedge == -1)
                {
                    defedge = edgeId;
                    break;
                }
            }
        }
        if (defedge != -1)
        {
            break;
        }
    }
 
    for (size_t n = 1; n < nlevels; ++n)
    {
        // produce a map with a key that is the element id
        // that contains which next level patch it belongs to
        std::vector<MultiRegions::SubGraphSharedPtr> bndgraphs =
            bottomUpGraph->GetInteriorBlocks(n + 1);
 
        // Fill next level graph  of adjacent elements and their level
        std::map<int, int> ElmtInBndry;
 
        for (i = 0; i < bndgraphs.size(); ++i)
        {
            int GlobIdOffset = bndgraphs[i]->GetIdOffset();
            size_t nVert     = bndgraphs[i]->GetNverts();
 
            for (j = 0; j < nVert; ++j)
            {
                // find edge in graph vert list
                if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
                {
                    edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
 
                    if (EdgeIdToElmts[edgeId][0] != -1)
                    {
                        ElmtInBndry[EdgeIdToElmts[edgeId][0]] = i;
                    }
                    if (EdgeIdToElmts[edgeId][1] != -1)
                    {
                        ElmtInBndry[EdgeIdToElmts[edgeId][1]] = i;
                    }
                }
            }
        }
 
        // Now search interior patches in this level for edges
        // that share the same element as a boundary edge and
        // assign this elmt that boundary edge
        std::vector<MultiRegions::SubGraphSharedPtr> intgraphs =
            bottomUpGraph->GetInteriorBlocks(n);
        for (i = 0; i < intgraphs.size(); ++i)
        {
            int GlobIdOffset = intgraphs[i]->GetIdOffset();
            bool SetEdge     = false;
            int elmtid       = 0;
            size_t nVert     = intgraphs[i]->GetNverts();
            for (j = 0; j < nVert; ++j)
            {
                // Check to see if graph vert is an edge
                if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
                {
                    edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
 
                    for (k = 0; k < 2; ++k)
                    {
                        // relevant edge id
                        elmtid = EdgeIdToElmts[edgeId][k];
 
                        if (elmtid != -1)
                        {
                            auto mapIt = ElmtInBndry.find(elmtid);
 
                            if (mapIt != ElmtInBndry.end())
                            {
                                // now find a edge in the next level boundary
                                // graph
                                int GlobIdOffset1 =
                                    bndgraphs[mapIt->second]->GetIdOffset();
                                for (int l = 0;
                                     l < bndgraphs[mapIt->second]->GetNverts();
                                     ++l)
                                {
                                    // find edge in graph vert list
                                    if (HomGraphEdgeIdToEdgeId.count(
                                            GlobIdOffset1 + l) != 0)
                                    {
                                        AddMeanPressureToEdgeId[elmtid] =
                                            defedge;
 
                                        SetEdge = true;
                                        break;
                                    }
                                }
                            }
                        }
                    }
                }
            }
 
            // if we have failed to find matching edge in next
            // level patch boundary then set last found elmt
            // associated to this interior patch to the
            // default edget value
            if (SetEdge == false)
            {
                if (elmtid == -1) // find an elmtid in patch
                {
                    size_t nVert = intgraphs[i]->GetNverts();
                    for (j = 0; j < nVert; ++j)
                    {
                        if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
                        {
                            edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
                            for (k = 0; k < 2; ++k)
                            {
                                // relevant edge id
                                elmtid = EdgeIdToElmts[edgeId][k];
                                if (elmtid != -1)
                                {
                                    break;
                                }
                            }
                        }
                        if (elmtid != -1)
                        {
                            break;
                        }
                    }
                }
                if (AddMeanPressureToEdgeId[elmtid] == -1)
                {
                    AddMeanPressureToEdgeId[elmtid] = defedge;
                }
            }
        }
    }
}
 
} // namespace Nektar