AssemblyMap.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: AssemblyMap.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: Assembly (e.g. local to global) base mapping routines
//
///////////////////////////////////////////////////////////////////////////////
 
#include <boost/core/ignore_unused.hpp>
 
#include <MultiRegions/AssemblyMap/AssemblyMap.h>
 
using namespace std;
 
namespace Nektar
{
namespace MultiRegions
{
/**
 * @class AssemblyMap
 * This class acts as a base class for constructing mappings between
 * local, global and boundary degrees of freedom. It holds the storage
 * for the maps and provides the accessors needed to retrieve them.
 *
 * There are two derived classes: AssemblyMapCG and
 * AssemblyMapDG. These perform the actual construction of the
 * maps within their specific contexts.
 *
 */
 
/// Rounds a double precision number to an integer.
int RoundNekDoubleToInt(NekDouble x)
{
    return int(x > 0.0 ? x + 0.5 : x - 0.5);
}
 
/// Rounds an array of double precision numbers to integers.
void RoundNekDoubleToInt(const Array<OneD, const NekDouble> inarray,
                         Array<OneD, int> outarray)
{
    int size = inarray.size();
    ASSERTL1(outarray.size() >= size, "Array sizes not compatible");
 
    NekDouble x;
    for (int i = 0; i < size; i++)
    {
        x           = inarray[i];
        outarray[i] = int(x > 0.0 ? x + 0.5 : x - 0.5);
    }
}
 
/**
 * Initialises an empty mapping.
 */
AssemblyMap::AssemblyMap()
    : m_session(), m_comm(), m_hash(0), m_numLocalBndCoeffs(0),
      m_numGlobalBndCoeffs(0), m_numLocalDirBndCoeffs(0),
      m_numGlobalDirBndCoeffs(0), m_solnType(eNoSolnType),
      m_bndSystemBandWidth(0), m_successiveRHS(0),
      m_linSysIterSolver("ConjugateGradient"), m_gsh(0), m_bndGsh(0)
{
}
 
AssemblyMap::AssemblyMap(const LibUtilities::SessionReaderSharedPtr &pSession,
                         const LibUtilities::CommSharedPtr &comm,
                         const std::string variable)
    : m_session(pSession), m_comm(comm), m_hash(0), m_numLocalBndCoeffs(0),
      m_numGlobalBndCoeffs(0), m_numLocalDirBndCoeffs(0),
      m_numGlobalDirBndCoeffs(0), m_bndSystemBandWidth(0), m_successiveRHS(0),
      m_linSysIterSolver("ConjugateGradient"), m_gsh(0), m_bndGsh(0)
{
    // Default value from Solver Info
    m_solnType =
        pSession->GetSolverInfoAsEnum<GlobalSysSolnType>("GlobalSysSoln");
    m_preconType =
        pSession->GetSolverInfoAsEnum<PreconditionerType>("Preconditioner");
 
    // Override values with data from GlobalSysSolnInfo section
    if (pSession->DefinesGlobalSysSolnInfo(variable, "GlobalSysSoln"))
    {
        std::string sysSoln =
            pSession->GetGlobalSysSolnInfo(variable, "GlobalSysSoln");
        m_solnType = pSession->GetValueAsEnum<GlobalSysSolnType>(
            "GlobalSysSoln", sysSoln);
    }
 
    if (pSession->DefinesGlobalSysSolnInfo(variable, "Preconditioner"))
    {
        std::string precon =
            pSession->GetGlobalSysSolnInfo(variable, "Preconditioner");
        m_preconType = pSession->GetValueAsEnum<PreconditionerType>(
            "Preconditioner", precon);
    }
 
    if (pSession->DefinesGlobalSysSolnInfo(variable,
                                           "IterativeSolverTolerance"))
    {
        m_iterativeTolerance = boost::lexical_cast<NekDouble>(
            pSession->GetGlobalSysSolnInfo(variable, "IterativeSolverTolerance")
                .c_str());
    }
    else
    {
        pSession->LoadParameter("IterativeSolverTolerance",
                                m_iterativeTolerance,
                                NekConstants::kNekIterativeTol);
    }
 
    if (pSession->DefinesGlobalSysSolnInfo(variable, "MaxIterations"))
    {
        m_maxIterations = boost::lexical_cast<int>(
            pSession->GetGlobalSysSolnInfo(variable, "MaxIterations").c_str());
    }
    else
    {
        pSession->LoadParameter("MaxIterations", m_maxIterations, 5000);
    }
 
    if (pSession->DefinesGlobalSysSolnInfo(variable, "SuccessiveRHS"))
    {
        m_successiveRHS = boost::lexical_cast<int>(
            pSession->GetGlobalSysSolnInfo(variable, "SuccessiveRHS").c_str());
    }
    else
    {
        pSession->LoadParameter("SuccessiveRHS", m_successiveRHS, 0);
    }
 
    if (pSession->DefinesGlobalSysSolnInfo(variable, "LinSysIterSolver"))
    {
        m_linSysIterSolver =
            pSession->GetGlobalSysSolnInfo(variable, "LinSysIterSolver");
    }
    else if (pSession->DefinesSolverInfo("LinSysIterSolver"))
    {
        m_linSysIterSolver = pSession->GetSolverInfo("LinSysIterSolver");
    }
    else
    {
        m_linSysIterSolver = "ConjugateGradient";
    }
}
 
/**
 * Create a new level of mapping using the information in
 * multiLevelGraph and performing the following steps:
 */
AssemblyMap::AssemblyMap(
    AssemblyMap *oldLevelMap,
    const BottomUpSubStructuredGraphSharedPtr &multiLevelGraph)
    : m_session(oldLevelMap->m_session), m_comm(oldLevelMap->GetComm()),
      m_hash(0), m_solnType(oldLevelMap->m_solnType),
      m_preconType(oldLevelMap->m_preconType),
      m_maxIterations(oldLevelMap->m_maxIterations),
      m_iterativeTolerance(oldLevelMap->m_iterativeTolerance),
      m_successiveRHS(oldLevelMap->m_successiveRHS),
      m_linSysIterSolver(oldLevelMap->m_linSysIterSolver),
      m_gsh(oldLevelMap->m_gsh), m_bndGsh(oldLevelMap->m_bndGsh),
      m_lowestStaticCondLevel(oldLevelMap->m_lowestStaticCondLevel)
{
    int i;
    int j;
    int cnt;
 
    //--------------------------------------------------------------
    // -- Extract information from the input argument
    int numGlobalBndCoeffsOld    = oldLevelMap->GetNumGlobalBndCoeffs();
    int numGlobalDirBndCoeffsOld = oldLevelMap->GetNumGlobalDirBndCoeffs();
    int numLocalBndCoeffsOld     = oldLevelMap->GetNumLocalBndCoeffs();
    int numLocalDirBndCoeffsOld  = oldLevelMap->GetNumLocalDirBndCoeffs();
    bool signChangeOld           = oldLevelMap->GetSignChange();
 
    int staticCondLevelOld        = oldLevelMap->GetStaticCondLevel();
    int numPatchesOld             = oldLevelMap->GetNumPatches();
    GlobalSysSolnType solnTypeOld = oldLevelMap->GetGlobalSysSolnType();
    const Array<OneD, const unsigned int> &numLocalBndCoeffsPerPatchOld =
        oldLevelMap->GetNumLocalBndCoeffsPerPatch();
    //--------------------------------------------------------------
 
    //--------------------------------------------------------------
    int newLevel = staticCondLevelOld + 1;
    /** - STEP 1: setup a mask array to determine to which patch
     *          of the new level every patch of the current
     *          level belongs.  To do so we make four arrays,
     *          #gloPatchMask, #globHomPatchMask,
     *          #locPatchMask_NekDouble and #locPatchMask.
     *          These arrays are then used to check which local
     *          dofs of the old level belong to which patch of
     *          the new level
     */
    Array<OneD, NekDouble> globPatchMask(numGlobalBndCoeffsOld, -1.0);
    Array<OneD, NekDouble> globHomPatchMask(globPatchMask +
                                            numGlobalDirBndCoeffsOld);
    Array<OneD, NekDouble> locPatchMask_NekDouble(numLocalBndCoeffsOld, -3.0);
    Array<OneD, int> locPatchMask(numLocalBndCoeffsOld);
 
    // Fill the array globPatchMask as follows:
    // - The first part (i.e. the glob bnd dofs) is filled with the
    //   value -1
    // - The second part (i.e. the glob interior dofs) is numbered
    //   according to the patch it belongs to (i.e. dofs in first block
    //   all are numbered 0, the second block numbered are 1, etc...)
    multiLevelGraph->MaskPatches(newLevel, globHomPatchMask);
 
    // Map from Global Dofs to Local Dofs
    // As a result, we know for each local dof whether
    // it is mapped to the boundary of the next level, or to which
    // patch. Based upon this, we can than later associate every patch
    // of the current level with a patch in the next level.
    oldLevelMap->GlobalToLocalBndWithoutSign(globPatchMask,
                                             locPatchMask_NekDouble);
 
    // Convert the result to an array of integers rather than doubles
    RoundNekDoubleToInt(locPatchMask_NekDouble, locPatchMask);
 
    /** - STEP 2: We calculate how many local bnd dofs of the
     *  old level belong to the boundaries of each patch at
     *  the new level. We need this information to set up the
     *  mapping between different levels.
     */
 
    // Retrieve the number of patches at the next level
    int numPatchesWithIntNew =
        multiLevelGraph->GetNpatchesWithInterior(newLevel);
    int numPatchesNew = numPatchesWithIntNew;
 
    // Allocate memory to store the number of local dofs associated to
    // each of elemental boundaries of these patches
    std::map<int, int> numLocalBndCoeffsPerPatchNew;
    for (int i = 0; i < numPatchesNew; i++)
    {
        numLocalBndCoeffsPerPatchNew[i] = 0;
    }
 
    int minval;
    int maxval;
    int curPatch;
    for (i = cnt = 0; i < numPatchesOld; i++)
    {
        // For every patch at the current level, the mask array
        // locPatchMask should be filled with
        // - the same (positive) number for each entry (which will
        //   correspond to the patch at the next level it belongs to)
        // - the same (positive) number for each entry, except some
        //   entries that are -1 (the enties correspond to -1, will be
        //   mapped to the local boundary of the next level patch given
        //   by the positive number)
        // - -1 for all entries. In this case, we will make an
        //   additional patch only consisting of boundaries at the next
        //   level
        minval =
            *min_element(&locPatchMask[cnt],
                         &locPatchMask[cnt] + numLocalBndCoeffsPerPatchOld[i]);
        maxval =
            *max_element(&locPatchMask[cnt],
                         &locPatchMask[cnt] + numLocalBndCoeffsPerPatchOld[i]);
        ASSERTL0((minval == maxval) || (minval == -1),
                 "These values should never be the same");
 
        if (maxval == -1)
        {
            curPatch                               = numPatchesNew;
            numLocalBndCoeffsPerPatchNew[curPatch] = 0;
            numPatchesNew++;
        }
        else
        {
            curPatch = maxval;
        }
 
        for (j = 0; j < numLocalBndCoeffsPerPatchOld[i]; j++)
        {
            ASSERTL0((locPatchMask[cnt] == maxval) ||
                         (locPatchMask[cnt] == minval),
                     "These values should never be the same");
            if (locPatchMask[cnt] == -1)
            {
                ++numLocalBndCoeffsPerPatchNew[curPatch];
            }
            cnt++;
        }
    }
 
    // Count how many local dofs of the old level are mapped
    // to the local boundary dofs of the new level
    m_numLocalCoeffs            = 0;
    m_numLocalBndCoeffs         = 0;
    m_numPatches                = numLocalBndCoeffsPerPatchNew.size();
    m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
    m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches, 0u);
    multiLevelGraph->GetNintDofsPerPatch(newLevel, m_numLocalIntCoeffsPerPatch);
 
    for (int i = 0; i < m_numPatches; i++)
    {
        m_numLocalBndCoeffsPerPatch[i] =
            (unsigned int)numLocalBndCoeffsPerPatchNew[i];
        m_numLocalBndCoeffs += m_numLocalBndCoeffsPerPatch[i];
        m_numLocalCoeffs +=
            m_numLocalBndCoeffsPerPatch[i] + m_numLocalIntCoeffsPerPatch[i];
    }
 
    // Also initialise some more data members
    m_solnType = solnTypeOld;
    ASSERTL1(m_solnType == eDirectMultiLevelStaticCond ||
                 m_solnType == eIterativeMultiLevelStaticCond ||
                 m_solnType == eXxtMultiLevelStaticCond ||
                 m_solnType == ePETScMultiLevelStaticCond,
             "This method should only be called for in "
             "case of multi-level static condensation.");
    m_staticCondLevel       = newLevel;
    m_signChange            = signChangeOld;
    m_numLocalDirBndCoeffs  = numLocalDirBndCoeffsOld;
    m_numGlobalDirBndCoeffs = numGlobalDirBndCoeffsOld;
    m_numGlobalBndCoeffs =
        multiLevelGraph->GetInteriorOffset(newLevel) + m_numGlobalDirBndCoeffs;
    m_numGlobalCoeffs =
        multiLevelGraph->GetNumGlobalDofs(newLevel) + m_numGlobalDirBndCoeffs;
    m_localToGlobalBndMap = Array<OneD, int>(m_numLocalBndCoeffs);
 
    m_localToLocalBndMap = Array<OneD, int>(m_numLocalBndCoeffs);
    m_localToLocalIntMap =
        Array<OneD, int>(m_numLocalCoeffs - m_numLocalBndCoeffs);
 
    // local to bnd map is just a copy
    for (int i = 0; i < m_numLocalBndCoeffs; ++i)
    {
        m_localToLocalBndMap[i] = i;
    }
 
    // local to int map is just a copy plus offset
    for (int i = m_numLocalBndCoeffs; i < m_numLocalCoeffs; ++i)
    {
        m_localToLocalIntMap[i - m_numLocalBndCoeffs] = i;
    }
 
    if (m_signChange)
    {
        m_localToGlobalBndSign = Array<OneD, NekDouble>(m_numLocalBndCoeffs);
    }
 
    m_patchMapFromPrevLevel =
        MemoryManager<PatchMap>::AllocateSharedPtr(numLocalBndCoeffsOld);
 
    m_globalToUniversalBndMap = Array<OneD, int>(
        m_numGlobalBndCoeffs, oldLevelMap->GetGlobalToUniversalBndMap());
    m_globalToUniversalBndMapUnique = Array<OneD, int>(
        m_numGlobalBndCoeffs, oldLevelMap->GetGlobalToUniversalBndMapUnique());
 
    // Set up an offset array that denotes the offset of the local
    // boundary degrees of freedom of the next level
    Array<OneD, int> numLocalBndCoeffsPerPatchOffset(m_numPatches + 1, 0);
    for (int i = 1; i < m_numPatches + 1; i++)
    {
        numLocalBndCoeffsPerPatchOffset[i] +=
            numLocalBndCoeffsPerPatchOffset[i - 1] +
            numLocalBndCoeffsPerPatchNew[i - 1];
    }
 
    int additionalPatchCnt = numPatchesWithIntNew;
    int newid;
    int blockid;
    bool isBndDof;
    NekDouble sign;
    Array<OneD, int> bndDofPerPatchCnt(m_numPatches, 0);
 
    for (i = cnt = 0; i < numPatchesOld; i++)
    {
        minval =
            *min_element(&locPatchMask[cnt],
                         &locPatchMask[cnt] + numLocalBndCoeffsPerPatchOld[i]);
        maxval =
            *max_element(&locPatchMask[cnt],
                         &locPatchMask[cnt] + numLocalBndCoeffsPerPatchOld[i]);
        ASSERTL0((minval == maxval) || (minval == -1),
                 "These values should never be the same");
 
        if (maxval == -1)
        {
            curPatch = additionalPatchCnt;
            additionalPatchCnt++;
        }
        else
        {
            curPatch = maxval;
        }
 
        for (j = 0; j < numLocalBndCoeffsPerPatchOld[i]; j++)
        {
            ASSERTL0((locPatchMask[cnt] == maxval) ||
                         (locPatchMask[cnt] == minval),
                     "These values should never be the same");
 
            sign = oldLevelMap->GetLocalToGlobalBndSign(cnt);
 
            if (locPatchMask[cnt] == -1)
            {
                newid = numLocalBndCoeffsPerPatchOffset[curPatch];
 
                m_localToGlobalBndMap[newid] =
                    oldLevelMap->GetLocalToGlobalBndMap(cnt);
 
                if (m_signChange)
                {
                    m_localToGlobalBndSign[newid] = sign;
                }
 
                blockid  = bndDofPerPatchCnt[curPatch];
                isBndDof = true;
 
                numLocalBndCoeffsPerPatchOffset[curPatch]++;
                bndDofPerPatchCnt[curPatch]++;
            }
            else
            {
                newid = oldLevelMap->GetLocalToGlobalBndMap(cnt) -
                        m_numGlobalBndCoeffs + m_numLocalBndCoeffs;
 
                blockid =
                    oldLevelMap->GetLocalToGlobalBndMap(cnt) -
                    m_numGlobalDirBndCoeffs -
                    multiLevelGraph->GetInteriorOffset(newLevel, curPatch);
                isBndDof = false;
            }
 
            sign = isBndDof ? 1.0 : sign;
 
            m_patchMapFromPrevLevel->SetPatchMap(cnt, curPatch, blockid,
                                                 isBndDof, sign);
            cnt++;
        }
    }
 
    // set up local to local mapping from previous to new level
    m_patchMapFromPrevLevel->SetNewLevelMap(m_numLocalBndCoeffsPerPatch,
                                            m_numLocalIntCoeffsPerPatch);
 
    CalculateBndSystemBandWidth();
 
    // Postprocess the computed information - Update the old
    // level with the mapping to new level
    // oldLevelMap->SetLocalBndToNextLevelMap(oldLocalBndToNewLevelMap,oldLocalBndToNewLevelSign);
 
    // - Construct the next level mapping object
    if (m_staticCondLevel < (multiLevelGraph->GetNlevels() - 1))
    {
        m_nextLevelLocalToGlobalMap =
            MemoryManager<AssemblyMap>::AllocateSharedPtr(this,
                                                          multiLevelGraph);
    }
}
 
AssemblyMap::~AssemblyMap(void)
{
}
 
/**
 * The bandwidth calculated corresponds to what is referred to as
 * half-bandwidth.  If the elements of the matrix are designated as
 * a_ij, it corresponds to the maximum value of |i-j| for non-zero
 * a_ij.  As a result, the value also corresponds to the number of
 * sub- or super-diagonals.
 *
 * The bandwith can be calculated elementally as it corresponds to the
 * maximal elemental bandwith (i.e. the maximal difference in global
 * DOF index for every element).
 *
 * We here calculate the bandwith of the global boundary system (as
 * used for static condensation).
 */
void AssemblyMap::CalculateBndSystemBandWidth()
{
    int i, j;
    int cnt = 0;
    int locSize;
    int maxId;
    int minId;
    int bwidth = -1;
    for (i = 0; i < m_numPatches; ++i)
    {
        locSize = m_numLocalBndCoeffsPerPatch[i];
        maxId   = -1;
        minId   = m_numLocalBndCoeffs + 1;
        for (j = 0; j < locSize; j++)
        {
            if (m_localToGlobalBndMap[cnt + j] >= m_numGlobalDirBndCoeffs)
            {
                if (m_localToGlobalBndMap[cnt + j] > maxId)
                {
                    maxId = m_localToGlobalBndMap[cnt + j];
                }
 
                if (m_localToGlobalBndMap[cnt + j] < minId)
                {
                    minId = m_localToGlobalBndMap[cnt + j];
                }
            }
        }
        bwidth = (bwidth > (maxId - minId)) ? bwidth : (maxId - minId);
 
        cnt += locSize;
    }
 
    m_bndSystemBandWidth = bwidth;
}
 
int AssemblyMap::v_GetLocalToGlobalMap(const int i) const
{
    boost::ignore_unused(i);
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    return 0;
}
 
int AssemblyMap::v_GetGlobalToUniversalMap(const int i) const
{
    boost::ignore_unused(i);
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    return 0;
}
 
int AssemblyMap::v_GetGlobalToUniversalMapUnique(const int i) const
{
    boost::ignore_unused(i);
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    return 0;
}
 
const Array<OneD, const int> &AssemblyMap::v_GetLocalToGlobalMap()
{
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    static Array<OneD, const int> result;
    return result;
}
 
const Array<OneD, const int> &AssemblyMap::v_GetGlobalToUniversalMap()
{
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    static Array<OneD, const int> result;
    return result;
}
 
const Array<OneD, const int> &AssemblyMap::v_GetGlobalToUniversalMapUnique()
{
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    static Array<OneD, const int> result;
    return result;
}
 
NekDouble AssemblyMap::v_GetLocalToGlobalSign(const int i) const
{
    boost::ignore_unused(i);
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    return 0.0;
}
 
const Array<OneD, NekDouble> &AssemblyMap::v_GetLocalToGlobalSign() const
{
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    static Array<OneD, NekDouble> result;
    return result;
}
 
void AssemblyMap::v_LocalToGlobal(const Array<OneD, const NekDouble> &loc,
                                  Array<OneD, NekDouble> &global,
                                  bool useComm) const
{
    boost::ignore_unused(loc, global, useComm);
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
}
 
void AssemblyMap::v_LocalToGlobal(const NekVector<NekDouble> &loc,
                                  NekVector<NekDouble> &global,
                                  bool useComm) const
{
    boost::ignore_unused(loc, global, useComm);
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
}
 
void AssemblyMap::v_GlobalToLocal(const Array<OneD, const NekDouble> &global,
                                  Array<OneD, NekDouble> &loc) const
{
    boost::ignore_unused(loc, global);
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
}
 
void AssemblyMap::v_GlobalToLocal(const NekVector<NekDouble> &global,
                                  NekVector<NekDouble> &loc) const
{
    boost::ignore_unused(loc, global);
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
}
 
void AssemblyMap::v_Assemble(const Array<OneD, const NekDouble> &loc,
                             Array<OneD, NekDouble> &global) const
{
    boost::ignore_unused(loc, global);
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
}
 
void AssemblyMap::v_Assemble(const NekVector<NekDouble> &loc,
                             NekVector<NekDouble> &global) const
{
    boost::ignore_unused(loc, global);
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
}
 
void AssemblyMap::v_UniversalAssemble(Array<OneD, NekDouble> &pGlobal) const
{
    boost::ignore_unused(pGlobal);
    // Do nothing here since multi-level static condensation uses a
    // AssemblyMap and thus will call this routine in serial.
}
 
void AssemblyMap::v_UniversalAssemble(NekVector<NekDouble> &pGlobal) const
{
    boost::ignore_unused(pGlobal);
    // Do nothing here since multi-level static condensation uses a
    // AssemblyMap and thus will call this routine in serial.
}
 
void AssemblyMap::v_UniversalAssemble(Array<OneD, NekDouble> &pGlobal,
                                      int offset) const
{
    boost::ignore_unused(pGlobal, offset);
    // Do nothing here since multi-level static condensation uses a
    // AssemblyMap and thus will call this routine in serial.
}
 
int AssemblyMap::v_GetFullSystemBandWidth() const
{
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    return 0;
}
 
int AssemblyMap::v_GetNumNonDirVertexModes() const
{
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    return 0;
}
 
int AssemblyMap::v_GetNumNonDirEdgeModes() const
{
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    return 0;
}
 
int AssemblyMap::v_GetNumNonDirFaceModes() const
{
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    return 0;
}
 
int AssemblyMap::v_GetNumDirEdges() const
{
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    return 0;
}
 
int AssemblyMap::v_GetNumDirFaces() const
{
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    return 0;
}
 
int AssemblyMap::v_GetNumNonDirEdges() const
{
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    return 0;
}
 
int AssemblyMap::v_GetNumNonDirFaces() const
{
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    return 0;
}
 
const Array<OneD, const int> &AssemblyMap::v_GetExtraDirEdges()
{
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    static Array<OneD, const int> result;
    return result;
}
 
std::shared_ptr<AssemblyMap> AssemblyMap::v_LinearSpaceMap(
    const ExpList &locexp, GlobalSysSolnType solnType)
{
    boost::ignore_unused(locexp, solnType);
    NEKERROR(ErrorUtil::efatal, "Not defined for this type of mapping.");
    static std::shared_ptr<AssemblyMap> result;
    return result;
}
 
LibUtilities::CommSharedPtr AssemblyMap::GetComm()
{
    return m_comm;
}
 
size_t AssemblyMap::GetHash() const
{
    return m_hash;
}
 
int AssemblyMap::GetLocalToGlobalMap(const int i) const
{
    return v_GetLocalToGlobalMap(i);
}
 
int AssemblyMap::GetGlobalToUniversalMap(const int i) const
{
    return v_GetGlobalToUniversalMap(i);
}
 
int AssemblyMap::GetGlobalToUniversalMapUnique(const int i) const
{
    return v_GetGlobalToUniversalMapUnique(i);
}
 
const Array<OneD, const int> &AssemblyMap::GetLocalToGlobalMap()
{
    return v_GetLocalToGlobalMap();
}
 
const Array<OneD, const int> &AssemblyMap::GetGlobalToUniversalMap()
{
    return v_GetGlobalToUniversalMap();
}
 
const Array<OneD, const int> &AssemblyMap::GetGlobalToUniversalMapUnique()
{
    return v_GetGlobalToUniversalMapUnique();
}
 
NekDouble AssemblyMap::GetLocalToGlobalSign(const int i) const
{
    return v_GetLocalToGlobalSign(i);
}
 
const Array<OneD, NekDouble> &AssemblyMap::GetLocalToGlobalSign() const
{
    return v_GetLocalToGlobalSign();
}
 
void AssemblyMap::LocalToGlobal(const Array<OneD, const NekDouble> &loc,
                                Array<OneD, NekDouble> &global,
                                bool useComm) const
{
    v_LocalToGlobal(loc, global, useComm);
}
 
void AssemblyMap::LocalToGlobal(const NekVector<NekDouble> &loc,
                                NekVector<NekDouble> &global,
                                bool useComm) const
{
    v_LocalToGlobal(loc, global, useComm);
}
 
void AssemblyMap::GlobalToLocal(const Array<OneD, const NekDouble> &global,
                                Array<OneD, NekDouble> &loc) const
{
    v_GlobalToLocal(global, loc);
}
 
void AssemblyMap::GlobalToLocal(const NekVector<NekDouble> &global,
                                NekVector<NekDouble> &loc) const
{
    v_GlobalToLocal(global, loc);
}
 
void AssemblyMap::Assemble(const Array<OneD, const NekDouble> &loc,
                           Array<OneD, NekDouble> &global) const
{
    v_Assemble(loc, global);
}
 
void AssemblyMap::Assemble(const NekVector<NekDouble> &loc,
                           NekVector<NekDouble> &global) const
{
    v_Assemble(loc, global);
}
 
void AssemblyMap::UniversalAssemble(Array<OneD, NekDouble> &pGlobal) const
{
    v_UniversalAssemble(pGlobal);
}
 
void AssemblyMap::UniversalAssemble(NekVector<NekDouble> &pGlobal) const
{
    v_UniversalAssemble(pGlobal);
}
 
void AssemblyMap::UniversalAssemble(Array<OneD, NekDouble> &pGlobal,
                                    int offset) const
{
    v_UniversalAssemble(pGlobal, offset);
}
 
void AssemblyMap::PatchLocalToGlobal(const Array<OneD, const NekDouble> &loc,
                                     Array<OneD, NekDouble> &global) const
{
    Array<OneD, const NekDouble> local;
 
    Array<OneD, const int> map = m_patchMapFromPrevLevel->GetNewLevelMap();
    Array<OneD, const NekDouble> sign = m_patchMapFromPrevLevel->GetSign();
 
    if (global.data() == loc.data())
    {
        local = Array<OneD, NekDouble>(map.size(), loc.data());
    }
    else
    {
        local = loc; // create reference
    }
 
    Vmath::Scatr(map.size(), sign.get(), local.get(), map.get(), global.get());
}
 
void AssemblyMap::PatchGlobalToLocal(const Array<OneD, const NekDouble> &global,
                                     Array<OneD, NekDouble> &loc) const
{
    Array<OneD, const NekDouble> glo;
 
    Array<OneD, const int> map = m_patchMapFromPrevLevel->GetNewLevelMap();
    Array<OneD, const NekDouble> sign = m_patchMapFromPrevLevel->GetSign();
 
    if (global.data() == loc.data())
    {
        glo = Array<OneD, NekDouble>(global.size(), global.data());
    }
    else
    {
        glo = global; // create reference
    }
 
    Vmath::Gathr(map.size(), sign.get(), glo.get(), map.get(), loc.get());
}
 
void AssemblyMap::PatchAssemble(const Array<OneD, const NekDouble> &loc,
                                Array<OneD, NekDouble> &global) const
{
    Array<OneD, const NekDouble> local;
    Array<OneD, const int> map = m_patchMapFromPrevLevel->GetNewLevelMap();
    Array<OneD, const NekDouble> sign = m_patchMapFromPrevLevel->GetSign();
 
    if (global.data() == loc.data())
    {
        local = Array<OneD, NekDouble>(map.size(), loc.data());
    }
    else
    {
        local = loc; // create reference
    }
 
    // since we are calling mapping from level down from array
    // the m_numLocaBndCoeffs represents the size of the
    // boundary elements we need to assemble into
    Vmath::Zero(m_numLocalCoeffs, global.get(), 1);
 
    Vmath::Assmb(map.size(), sign.get(), local.get(), map.get(), global.get());
}
 
int AssemblyMap::GetFullSystemBandWidth() const
{
    return v_GetFullSystemBandWidth();
}
 
int AssemblyMap::GetNumNonDirVertexModes() const
{
    return v_GetNumNonDirVertexModes();
}
 
int AssemblyMap::GetNumNonDirEdgeModes() const
{
    return v_GetNumNonDirEdgeModes();
}
 
int AssemblyMap::GetNumNonDirFaceModes() const
{
    return v_GetNumNonDirFaceModes();
}
 
int AssemblyMap::GetNumDirEdges() const
{
    return v_GetNumDirEdges();
}
 
int AssemblyMap::GetNumDirFaces() const
{
    return v_GetNumDirFaces();
}
 
int AssemblyMap::GetNumNonDirEdges() const
{
    return v_GetNumNonDirEdges();
}
 
int AssemblyMap::GetNumNonDirFaces() const
{
    return v_GetNumNonDirFaces();
}
 
const Array<OneD, const int> &AssemblyMap::GetExtraDirEdges()
{
    return v_GetExtraDirEdges();
}
 
std::shared_ptr<AssemblyMap> AssemblyMap::LinearSpaceMap(
    const ExpList &locexp, GlobalSysSolnType solnType)
{
    return v_LinearSpaceMap(locexp, solnType);
}
 
int AssemblyMap::GetLocalToGlobalBndMap(const int i) const
{
    return m_localToGlobalBndMap[i];
}
 
const Array<OneD, const int> &AssemblyMap::GetLocalToGlobalBndMap(void)
{
    return m_localToGlobalBndMap;
}
 
bool AssemblyMap::GetSignChange()
{
    return m_signChange;
}
 
Array<OneD, const NekDouble> AssemblyMap::GetLocalToGlobalBndSign(void) const
{
    return m_localToGlobalBndSign;
}
 
const Array<OneD, const int> &AssemblyMap::GetGlobalToUniversalBndMap()
{
    return m_globalToUniversalBndMap;
}
 
const Array<OneD, const int> &AssemblyMap::GetGlobalToUniversalBndMapUnique()
{
    return m_globalToUniversalBndMapUnique;
}
 
NekDouble AssemblyMap::GetLocalToGlobalBndSign(const int i) const
{
    if (m_signChange)
    {
        return m_localToGlobalBndSign[i];
    }
    else
    {
        return 1.0;
    }
}
 
const Array<OneD, const int> &AssemblyMap::GetBndCondCoeffsToLocalCoeffsMap()
{
    return m_bndCondCoeffsToLocalCoeffsMap;
}
 
const Array<OneD, NekDouble> &AssemblyMap::GetBndCondCoeffsToLocalCoeffsSign()
{
    return m_bndCondCoeffsToLocalCoeffsSign;
}
 
const Array<OneD, const int> &AssemblyMap::GetBndCondCoeffsToLocalTraceMap()
{
    return m_bndCondCoeffsToLocalTraceMap;
}
 
int AssemblyMap::GetBndCondIDToGlobalTraceID(const int i)
{
    ASSERTL1(i < m_bndCondIDToGlobalTraceID.size(), "Index out of range.");
    return m_bndCondIDToGlobalTraceID[i];
}
 
const Array<OneD, const int> &AssemblyMap::GetBndCondIDToGlobalTraceID()
{
    return m_bndCondIDToGlobalTraceID;
}
 
int AssemblyMap::GetNumGlobalDirBndCoeffs() const
{
    return m_numGlobalDirBndCoeffs;
}
 
int AssemblyMap::GetNumLocalDirBndCoeffs() const
{
    return m_numLocalDirBndCoeffs;
}
 
int AssemblyMap::GetNumLocalBndCoeffs() const
{
    return m_numLocalBndCoeffs;
}
 
int AssemblyMap::GetNumGlobalBndCoeffs() const
{
    return m_numGlobalBndCoeffs;
}
 
int AssemblyMap::GetNumLocalCoeffs() const
{
    return m_numLocalCoeffs;
}
 
int AssemblyMap::GetNumGlobalCoeffs() const
{
    return m_numGlobalCoeffs;
}
 
bool AssemblyMap::GetSingularSystem() const
{
    return m_systemSingular;
}
 
void AssemblyMap::GlobalToLocalBnd(const NekVector<NekDouble> &global,
                                   NekVector<NekDouble> &loc, int offset) const
{
    GlobalToLocalBnd(global.GetPtr(), loc.GetPtr(), offset);
}
 
void AssemblyMap::GlobalToLocalBnd(const NekVector<NekDouble> &global,
                                   NekVector<NekDouble> &loc) const
{
    GlobalToLocalBnd(global.GetPtr(), loc.GetPtr());
}
 
void AssemblyMap::GlobalToLocalBnd(const Array<OneD, const NekDouble> &global,
                                   Array<OneD, NekDouble> &loc,
                                   int offset) const
{
    ASSERTL1(loc.size() >= m_numLocalBndCoeffs,
             "Local vector is not of correct dimension");
    ASSERTL1(global.size() >= m_numGlobalBndCoeffs - offset,
             "Global vector is not of correct dimension");
 
    // offset input data by length "offset" for Dirichlet boundary conditions.
    Array<OneD, NekDouble> tmp(m_numGlobalBndCoeffs, 0.0);
    Vmath::Vcopy(m_numGlobalBndCoeffs - offset, global.get(), 1,
                 tmp.get() + offset, 1);
 
    if (m_signChange)
    {
        Vmath::Gathr(m_numLocalBndCoeffs, m_localToGlobalBndSign.get(),
                     tmp.get(), m_localToGlobalBndMap.get(), loc.get());
    }
    else
    {
        Vmath::Gathr(m_numLocalBndCoeffs, tmp.get(),
                     m_localToGlobalBndMap.get(), loc.get());
    }
}
 
void AssemblyMap::GlobalToLocalBnd(const Array<OneD, const NekDouble> &global,
                                   Array<OneD, NekDouble> &loc) const
{
    ASSERTL1(loc.size() >= m_numLocalBndCoeffs,
             "Local vector is not of correct dimension");
    ASSERTL1(global.size() >= m_numGlobalBndCoeffs,
             "Global vector is not of correct dimension");
 
    if (m_signChange)
    {
        Vmath::Gathr(m_numLocalBndCoeffs, m_localToGlobalBndSign.get(),
                     global.get(), m_localToGlobalBndMap.get(), loc.get());
    }
    else
    {
        Vmath::Gathr(m_numLocalBndCoeffs, global.get(),
                     m_localToGlobalBndMap.get(), loc.get());
    }
}
 
void AssemblyMap::LocalBndToGlobal(const Array<OneD, const NekDouble> &loc,
                                   Array<OneD, NekDouble> &global, int offset,
                                   bool UseComm) const
{
    ASSERTL1(loc.size() >= m_numLocalBndCoeffs,
             "Local vector is not of correct dimension");
    ASSERTL1(global.size() >= m_numGlobalBndCoeffs - offset,
             "Global vector is not of correct dimension");
 
    // offset input data by length "offset" for Dirichlet boundary conditions.
    Array<OneD, NekDouble> tmp(m_numGlobalBndCoeffs, 0.0);
 
    if (m_signChange)
    {
        Vmath::Scatr(m_numLocalBndCoeffs, m_localToGlobalBndSign.get(),
                     loc.get(), m_localToGlobalBndMap.get(), tmp.get());
    }
    else
    {
        Vmath::Scatr(m_numLocalBndCoeffs, loc.get(),
                     m_localToGlobalBndMap.get(), tmp.get());
    }
 
    // Ensure each processor has unique value with a max gather.
    if (UseComm)
    {
        Gs::Gather(tmp, Gs::gs_max, m_bndGsh);
    }
    Vmath::Vcopy(m_numGlobalBndCoeffs - offset, tmp.get() + offset, 1,
                 global.get(), 1);
}
 
void AssemblyMap::LocalBndToGlobal(const Array<OneD, const NekDouble> &loc,
                                   Array<OneD, NekDouble> &global,
                                   bool UseComm) const
{
    ASSERTL1(loc.size() >= m_numLocalBndCoeffs,
             "Local vector is not of correct dimension");
    ASSERTL1(global.size() >= m_numGlobalBndCoeffs,
             "Global vector is not of correct dimension");
 
    if (m_signChange)
    {
        Vmath::Scatr(m_numLocalBndCoeffs, m_localToGlobalBndSign.get(),
                     loc.get(), m_localToGlobalBndMap.get(), global.get());
    }
    else
    {
        Vmath::Scatr(m_numLocalBndCoeffs, loc.get(),
                     m_localToGlobalBndMap.get(), global.get());
    }
    if (UseComm)
    {
        Gs::Gather(global, Gs::gs_max, m_bndGsh);
    }
}
 
void AssemblyMap::LocalToLocalBnd(const Array<OneD, const NekDouble> &local,
                                  Array<OneD, NekDouble> &locbnd) const
{
    ASSERTL1(locbnd.size() >= m_numLocalBndCoeffs,
             "LocBnd vector is not of correct dimension");
    ASSERTL1(local.size() >= m_numLocalCoeffs,
             "Local vector is not of correct dimension");
 
    Vmath::Gathr(m_numLocalBndCoeffs, local.get(), m_localToLocalBndMap.get(),
                 locbnd.get());
}
 
void AssemblyMap::LocalToLocalInt(const Array<OneD, const NekDouble> &local,
                                  Array<OneD, NekDouble> &locint) const
{
    ASSERTL1(locint.size() >= m_numLocalCoeffs - m_numLocalBndCoeffs,
             "Locint vector is not of correct dimension");
    ASSERTL1(local.size() >= m_numLocalCoeffs,
             "Local vector is not of correct dimension");
 
    Vmath::Gathr(m_numLocalCoeffs - m_numLocalBndCoeffs, local.get(),
                 m_localToLocalIntMap.get(), locint.get());
}
 
void AssemblyMap::LocalBndToLocal(const Array<OneD, const NekDouble> &locbnd,
                                  Array<OneD, NekDouble> &local) const
{
    ASSERTL1(locbnd.size() >= m_numLocalBndCoeffs,
             "LocBnd vector is not of correct dimension");
    ASSERTL1(local.size() >= m_numLocalCoeffs,
             "Local vector is not of correct dimension");
 
    Vmath::Scatr(m_numLocalBndCoeffs, locbnd.get(), m_localToLocalBndMap.get(),
                 local.get());
}
 
void AssemblyMap::LocalIntToLocal(const Array<OneD, const NekDouble> &locint,
                                  Array<OneD, NekDouble> &local) const
{
    ASSERTL1(locint.size() >= m_numLocalCoeffs - m_numLocalBndCoeffs,
             "LocBnd vector is not of correct dimension");
    ASSERTL1(local.size() >= m_numLocalCoeffs,
             "Local vector is not of correct dimension");
 
    Vmath::Scatr(m_numLocalCoeffs - m_numLocalBndCoeffs, locint.get(),
                 m_localToLocalIntMap.get(), local.get());
}
 
void AssemblyMap::AssembleBnd(const NekVector<NekDouble> &loc,
                              NekVector<NekDouble> &global, int offset) const
{
    AssembleBnd(loc.GetPtr(), global.GetPtr(), offset);
}
 
void AssemblyMap::AssembleBnd(const NekVector<NekDouble> &loc,
                              NekVector<NekDouble> &global) const
{
    AssembleBnd(loc.GetPtr(), global.GetPtr());
}
 
void AssemblyMap::AssembleBnd(const Array<OneD, const NekDouble> &loc,
                              Array<OneD, NekDouble> &global, int offset) const
{
    ASSERTL1(loc.size() >= m_numLocalBndCoeffs,
             "Local array is not of correct dimension");
    ASSERTL1(global.size() >= m_numGlobalBndCoeffs - offset,
             "Global array is not of correct dimension");
    Array<OneD, NekDouble> tmp(m_numGlobalBndCoeffs, 0.0);
 
    if (m_signChange)
    {
        Vmath::Assmb(m_numLocalBndCoeffs, m_localToGlobalBndSign.get(),
                     loc.get(), m_localToGlobalBndMap.get(), tmp.get());
    }
    else
    {
        Vmath::Assmb(m_numLocalBndCoeffs, loc.get(),
                     m_localToGlobalBndMap.get(), tmp.get());
    }
    UniversalAssembleBnd(tmp);
    Vmath::Vcopy(m_numGlobalBndCoeffs - offset, tmp.get() + offset, 1,
                 global.get(), 1);
}
 
void AssemblyMap::AssembleBnd(const Array<OneD, const NekDouble> &loc,
                              Array<OneD, NekDouble> &global) const
{
    ASSERTL1(loc.size() >= m_numLocalBndCoeffs,
             "Local vector is not of correct dimension");
    ASSERTL1(global.size() >= m_numGlobalBndCoeffs,
             "Global vector is not of correct dimension");
 
    Vmath::Zero(m_numGlobalBndCoeffs, global.get(), 1);
 
    if (m_signChange)
    {
        Vmath::Assmb(m_numLocalBndCoeffs, m_localToGlobalBndSign.get(),
                     loc.get(), m_localToGlobalBndMap.get(), global.get());
    }
    else
    {
        Vmath::Assmb(m_numLocalBndCoeffs, loc.get(),
                     m_localToGlobalBndMap.get(), global.get());
    }
    UniversalAssembleBnd(global);
}
 
void AssemblyMap::UniversalAssembleBnd(Array<OneD, NekDouble> &pGlobal) const
{
    ASSERTL1(pGlobal.size() >= m_numGlobalBndCoeffs, "Wrong size.");
    Gs::Gather(pGlobal, Gs::gs_add, m_bndGsh);
}
 
void AssemblyMap::UniversalAssembleBnd(NekVector<NekDouble> &pGlobal) const
{
    UniversalAssembleBnd(pGlobal.GetPtr());
}
 
void AssemblyMap::UniversalAssembleBnd(Array<OneD, NekDouble> &pGlobal,
                                       int offset) const
{
    Array<OneD, NekDouble> tmp(offset);
    if (offset > 0)
        Vmath::Vcopy(offset, pGlobal, 1, tmp, 1);
    UniversalAssembleBnd(pGlobal);
    if (offset > 0)
        Vmath::Vcopy(offset, tmp, 1, pGlobal, 1);
}
 
void AssemblyMap::UniversalAbsMaxBnd(Array<OneD, NekDouble> &bndvals)
{
    Gs::Gather(bndvals, Gs::gs_amax, m_dirBndGsh);
}
 
int AssemblyMap::GetBndSystemBandWidth() const
{
    return m_bndSystemBandWidth;
}
 
int AssemblyMap::GetStaticCondLevel() const
{
    return m_staticCondLevel;
}
 
int AssemblyMap::GetNumPatches() const
{
    return m_numPatches;
}
 
const Array<OneD, const unsigned int>
    &AssemblyMap::GetNumLocalBndCoeffsPerPatch()
{
    return m_numLocalBndCoeffsPerPatch;
}
 
const Array<OneD, const unsigned int>
    &AssemblyMap::GetNumLocalIntCoeffsPerPatch()
{
    return m_numLocalIntCoeffsPerPatch;
}
 
const AssemblyMapSharedPtr AssemblyMap::GetNextLevelLocalToGlobalMap() const
{
    return m_nextLevelLocalToGlobalMap;
}
 
const PatchMapSharedPtr &AssemblyMap::GetPatchMapFromPrevLevel(void) const
{
    return m_patchMapFromPrevLevel;
}
 
bool AssemblyMap::AtLastLevel() const
{
    return !(m_nextLevelLocalToGlobalMap.get());
}
 
GlobalSysSolnType AssemblyMap::GetGlobalSysSolnType() const
{
    return m_solnType;
}
 
PreconditionerType AssemblyMap::GetPreconType() const
{
    return m_preconType;
}
 
NekDouble AssemblyMap::GetIterativeTolerance() const
{
    return m_iterativeTolerance;
}
 
int AssemblyMap::GetMaxIterations() const
{
    return m_maxIterations;
}
 
int AssemblyMap::GetSuccessiveRHS() const
{
    return m_successiveRHS;
}
 
std::string AssemblyMap::GetLinSysIterSolver() const
{
    return m_linSysIterSolver;
}
 
void AssemblyMap::GlobalToLocalBndWithoutSign(
    const Array<OneD, const NekDouble> &global, Array<OneD, NekDouble> &loc)
{
    ASSERTL1(loc.size() >= m_numLocalBndCoeffs,
             "Local vector is not of correct dimension");
    ASSERTL1(global.size() >= m_numGlobalBndCoeffs,
             "Global vector is not of correct dimension");
 
    Vmath::Gathr(m_numLocalBndCoeffs, global.get(), m_localToGlobalBndMap.get(),
                 loc.get());
}
 
void AssemblyMap::PrintStats(std::ostream &out, std::string variable,
                             bool printHeader) const
{
    LibUtilities::CommSharedPtr vRowComm = m_session->GetComm()->GetRowComm();
    bool isRoot                          = vRowComm->GetRank() == 0;
    int n                                = vRowComm->GetSize();
    int i;
 
    // Determine number of global degrees of freedom.
    int globBndCnt = 0, globDirCnt = 0;
 
    for (i = 0; i < m_numGlobalBndCoeffs; ++i)
    {
        if (m_globalToUniversalBndMapUnique[i] > 0)
        {
            globBndCnt++;
 
            if (i < m_numGlobalDirBndCoeffs)
            {
                globDirCnt++;
            }
        }
    }
 
    int globCnt = m_numGlobalCoeffs - m_numGlobalBndCoeffs + globBndCnt;
 
    // Calculate maximum valency
    Array<OneD, NekDouble> tmpLoc(m_numLocalBndCoeffs, 1.0);
    Array<OneD, NekDouble> tmpGlob(m_numGlobalBndCoeffs, 0.0);
 
    Vmath::Assmb(m_numLocalBndCoeffs, tmpLoc.get(), m_localToGlobalBndMap.get(),
                 tmpGlob.get());
    UniversalAssembleBnd(tmpGlob);
 
    int totGlobDof     = globCnt;
    int totGlobBndDof  = globBndCnt;
    int totGlobDirDof  = globDirCnt;
    int totLocalDof    = m_numLocalCoeffs;
    int totLocalBndDof = m_numLocalBndCoeffs;
    int totLocalDirDof = m_numLocalDirBndCoeffs;
 
    int meanValence = 0;
    int maxValence  = 0;
    int minValence  = 10000000;
    for (int i = 0; i < m_numGlobalBndCoeffs; ++i)
    {
        if (!m_globalToUniversalBndMapUnique[i])
        {
            continue;
        }
 
        if (tmpGlob[i] > maxValence)
        {
            maxValence = tmpGlob[i];
        }
        if (tmpGlob[i] < minValence)
        {
            minValence = tmpGlob[i];
        }
        meanValence += tmpGlob[i];
    }
 
    vRowComm->AllReduce(maxValence, LibUtilities::ReduceMax);
    vRowComm->AllReduce(minValence, LibUtilities::ReduceMin);
    vRowComm->AllReduce(meanValence, LibUtilities::ReduceSum);
    vRowComm->AllReduce(totGlobDof, LibUtilities::ReduceSum);
    vRowComm->AllReduce(totGlobBndDof, LibUtilities::ReduceSum);
    vRowComm->AllReduce(totGlobDirDof, LibUtilities::ReduceSum);
    vRowComm->AllReduce(totLocalDof, LibUtilities::ReduceSum);
    vRowComm->AllReduce(totLocalBndDof, LibUtilities::ReduceSum);
    vRowComm->AllReduce(totLocalDirDof, LibUtilities::ReduceSum);
 
    meanValence /= totGlobBndDof;
 
    if (isRoot)
    {
        if (printHeader)
        {
            out << "Assembly map statistics for field " << variable << ":"
                << endl;
        }
 
        out << "  - Number of local/global dof             : " << totLocalDof
            << " " << totGlobDof << endl;
        out << "  - Number of local/global boundary dof    : " << totLocalBndDof
            << " " << totGlobBndDof << endl;
        out << "  - Number of local/global Dirichlet dof   : " << totLocalDirDof
            << " " << totGlobDirDof << endl;
        out << "  - dof valency (min/max/mean)             : " << minValence
            << " " << maxValence << " " << meanValence << endl;
 
        if (n > 1)
        {
            NekDouble mean = m_numLocalCoeffs, mean2 = mean * mean;
            NekDouble minval = mean, maxval = mean;
            Array<OneD, NekDouble> tmp(1);
 
            for (i = 1; i < n; ++i)
            {
                vRowComm->Recv(i, tmp);
                mean += tmp[0];
                mean2 += tmp[0] * tmp[0];
 
                if (tmp[0] > maxval)
                {
                    maxval = tmp[0];
                }
                if (tmp[0] < minval)
                {
                    minval = tmp[0];
                }
            }
 
            if (maxval > 0.1)
            {
                out << "  - Local dof dist. (min/max/mean/dev)     : " << minval
                    << " " << maxval << " " << (mean / n) << " "
                    << sqrt(mean2 / n - mean * mean / n / n) << endl;
            }
 
            vRowComm->Block();
 
            mean = minval = maxval = m_numLocalBndCoeffs;
            mean2                  = mean * mean;
 
            for (i = 1; i < n; ++i)
            {
                vRowComm->Recv(i, tmp);
                mean += tmp[0];
                mean2 += tmp[0] * tmp[0];
 
                if (tmp[0] > maxval)
                {
                    maxval = tmp[0];
                }
                if (tmp[0] < minval)
                {
                    minval = tmp[0];
                }
            }
 
            out << "  - Local bnd dof dist. (min/max/mean/dev) : " << minval
                << " " << maxval << " " << (mean / n) << " "
                << sqrt(mean2 / n - mean * mean / n / n) << endl;
        }
    }
    else
    {
        Array<OneD, NekDouble> tmp(1);
        tmp[0] = m_numLocalCoeffs;
        vRowComm->Send(0, tmp);
        vRowComm->Block();
        tmp[0] = m_numLocalBndCoeffs;
        vRowComm->Send(0, tmp);
    }
 
    // Either we have no more levels in the static condensation, or we
    // are not multi-level.
    if (!m_nextLevelLocalToGlobalMap)
    {
        return;
    }
 
    int level                = 2;
    AssemblyMapSharedPtr tmp = m_nextLevelLocalToGlobalMap;
    while (tmp->m_nextLevelLocalToGlobalMap)
    {
        tmp = tmp->m_nextLevelLocalToGlobalMap;
        ++level;
    }
 
    // Print out multi-level static condensation information.
    if (n > 1)
    {
        if (isRoot)
        {
            NekDouble mean = level, mean2 = mean * mean;
            int minval = level, maxval = level;
 
            Array<OneD, NekDouble> tmpRecv(1);
            for (i = 1; i < n; ++i)
            {
                vRowComm->Recv(i, tmpRecv);
                mean += tmpRecv[0];
                mean2 += tmpRecv[0] * tmpRecv[0];
 
                if (tmpRecv[0] > maxval)
                {
                    maxval = (int)(tmpRecv[0] + 0.5);
                }
                if (tmpRecv[0] < minval)
                {
                    minval = (int)(tmpRecv[0] + 0.5);
                }
            }
 
            out << "  - M-level sc. dist. (min/max/mean/dev)   : " << minval
                << " " << maxval << " " << (mean / n) << " "
                << sqrt(mean2 / n - mean * mean / n / n) << endl;
        }
        else
        {
            Array<OneD, NekDouble> tmpSend(1);
            tmpSend[0] = level;
            vRowComm->Send(0, tmpSend);
        }
    }
    else
    {
        out << "  - Number of static cond. levels          : " << level << endl;
    }
 
    if (isRoot)
    {
        out << "Stats at lowest static cond. level:" << endl;
    }
    tmp->PrintStats(out, variable, false);
}
} // namespace MultiRegions
} // namespace Nektar