DisContField.cpp

Go to the documentation of this file.

//////////////////////////////////////////////////////////////////////////////
//
// File: DisContField.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).
//
// License for the specific language governing rights and limitations under
// 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: Field definition for 1D domain with boundary
// conditions using LDG-H
//
///////////////////////////////////////////////////////////////////////////////
 
#include <LibUtilities/BasicUtils/ParseUtils.h>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <LocalRegions/Expansion0D.h>
#include <LocalRegions/Expansion1D.h>
#include <LocalRegions/HexExp.h>
#include <LocalRegions/QuadExp.h>
#include <LocalRegions/TetExp.h>
#include <LocalRegions/TriExp.h>
#include <MultiRegions/DisContField.h>
#include <SpatialDomains/MeshGraph.h>
#include <StdRegions/StdSegExp.h>
 
using namespace std;
 
namespace Nektar::MultiRegions
{
/**
 * @class DisContField
 * This class augments the list of local expansions inherited from
 * ExpList with boundary conditions. Inter-element boundaries are
 * handled using an discontinuous Galerkin scheme.
 */
 
/**
 * Constructs an empty expansion list with no boundary conditions.
 */
DisContField::DisContField()
    : ExpList(), m_bndConditions(), m_bndCondExpansions(),
      m_trace(NullExpListSharedPtr)
{
}
 
/**
 * An expansion list for the boundary expansions is generated first for
 * the field. These are subsequently evaluated for time zero. The trace
 * map is then constructed.
 * @param   graph1D     A mesh, containing information about the domain
 *                      and the spectral/hp element expansions.
 * @param   bcs         Information about the enforced boundary
 *                      conditions.
 * @param   variable    The session variable associated with the
 *                      boundary conditions to enforce.
 * @param   solnType    Type of global system to use.
 */
DisContField::DisContField(const LibUtilities::SessionReaderSharedPtr &pSession,
                           const SpatialDomains::MeshGraphSharedPtr &graph,
                           const std::string &variable, const bool SetUpJustDG,
                           const bool DeclareCoeffPhysArrays,
                           const Collections::ImplementationType ImpType,
                           const std::string bcvariable)
    : ExpList(pSession, graph, DeclareCoeffPhysArrays, variable, ImpType),
      m_trace(NullExpListSharedPtr)
{
    std::string bcvar;
    if (bcvariable == "NotSet")
    {
        bcvar = variable;
    }
    else
    {
        bcvar = bcvariable;
    }
 
    if (bcvar.compare("DefaultVar") != 0)
    {
        SpatialDomains::BoundaryConditions bcs(m_session, graph);
 
        GenerateBoundaryConditionExpansion(graph, bcs, bcvar,
                                           DeclareCoeffPhysArrays, ImpType);
        if (DeclareCoeffPhysArrays)
        {
            EvaluateBoundaryConditions(0.0, bcvar);
        }
        ApplyGeomInfo();
        FindPeriodicTraces(bcs, bcvar);
    }
 
    int i, cnt;
    Array<OneD, int> ElmtID, TraceID;
    GetBoundaryToElmtMap(ElmtID, TraceID);
 
    for (cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
    {
        MultiRegions::ExpListSharedPtr locExpList;
        locExpList = m_bndCondExpansions[i];
 
        for (int e = 0; e < locExpList->GetExpSize(); ++e)
        {
            (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(TraceID[cnt + e],
                                                   locExpList->GetExp(e));
            locExpList->GetExp(e)->SetAdjacentElementExp(
                TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
        }
        cnt += m_bndCondExpansions[i]->GetExpSize();
    }
 
    if (SetUpJustDG)
    {
        SetUpDG(variable, ImpType);
    }
    else
    {
        if (m_session->DefinesSolverInfo("PROJECTION"))
        {
            std::string ProjectStr = m_session->GetSolverInfo("PROJECTION");
            if ((ProjectStr == "MixedCGDG") ||
                (ProjectStr == "Mixed_CG_Discontinuous"))
            {
                SetUpDG(variable);
            }
            else
            {
                SetUpPhysNormals();
            }
        }
        else
        {
            SetUpPhysNormals();
        }
    }
}
 
/**
 * @brief Set up all DG member variables and maps.
 */
void DisContField::SetUpDG(const std::string variable,
                           const Collections::ImplementationType ImpType)
{
    // Check for multiple calls
    if (m_trace != NullExpListSharedPtr)
    {
        return;
    }
 
    // Set up matrix map
    m_globalBndMat = MemoryManager<GlobalLinSysMap>::AllocateSharedPtr();
 
    // Set up trace space
    m_trace = MemoryManager<ExpList>::AllocateSharedPtr(
        m_session, m_bndCondExpansions, m_bndConditions, *m_exp, m_graph,
        m_comm, true, "DefaultVar", ImpType);
 
    PeriodicMap periodicTraces = (m_expType == e1D)   ? m_periodicVerts
                                 : (m_expType == e2D) ? m_periodicEdges
                                                      : m_periodicFaces;
 
    m_traceMap = MemoryManager<AssemblyMapDG>::AllocateSharedPtr(
        m_session, m_graph, m_trace, *this, m_bndCondExpansions,
        m_bndConditions, periodicTraces, variable);
 
    if (m_session->DefinesCmdLineArgument("verbose"))
    {
        m_traceMap->PrintStats(std::cout, variable);
    }
 
    Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
        m_traceMap->GetElmtToTrace();
 
    // Scatter trace to elements. For each element, we find
    // the trace  associated to each elemental trace. The element
    // then retains a pointer to the elemental trace, to
    // ensure uniqueness of normals when retrieving from two
    // adjoining elements which do not lie in a plane.
    for (int i = 0; i < m_exp->size(); ++i)
    {
        for (int j = 0; j < (*m_exp)[i]->GetNtraces(); ++j)
        {
            LocalRegions::ExpansionSharedPtr exp = elmtToTrace[i][j];
            ;
 
            exp->SetAdjacentElementExp(j, (*m_exp)[i]);
            (*m_exp)[i]->SetTraceExp(j, exp);
        }
    }
 
    // Set up physical normals
    SetUpPhysNormals();
 
    // Create interface exchange object
    m_interfaceMap =
        MemoryManager<InterfaceMapDG>::AllocateSharedPtr(m_graph, m_trace);
 
    int cnt, n;
 
    // Identify boundary trace
    for (cnt = 0, n = 0; n < m_bndCondExpansions.size(); ++n)
    {
        if (m_bndConditions[n]->GetBoundaryConditionType() !=
            SpatialDomains::ePeriodic)
        {
            for (int v = 0; v < m_bndCondExpansions[n]->GetExpSize(); ++v)
            {
                m_boundaryTraces.insert(
                    m_traceMap->GetBndCondIDToGlobalTraceID(cnt + v));
            }
            cnt += m_bndCondExpansions[n]->GetExpSize();
        }
    }
 
    // Set up information for periodic boundary conditions.
    std::unordered_map<int, pair<int, int>> perTraceToExpMap;
    for (cnt = n = 0; n < m_exp->size(); ++n)
    {
        for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
        {
            auto it = periodicTraces.find((*m_exp)[n]->GetGeom()->GetTid(v));
 
            if (it != periodicTraces.end())
            {
                perTraceToExpMap[it->first] = make_pair(n, v);
            }
        }
    }
 
    // Set up left-adjacent tracelist.
    m_leftAdjacentTraces.resize(cnt);
 
    // count size of trace
    for (cnt = n = 0; n < m_exp->size(); ++n)
    {
        for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
        {
            m_leftAdjacentTraces[cnt] = IsLeftAdjacentTrace(n, v);
        }
    }
 
    // Set up mapping to copy Fwd of periodic bcs to Bwd of other edge.
    for (cnt = n = 0; n < m_exp->size(); ++n)
    {
        for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
        {
            int GeomId = (*m_exp)[n]->GetGeom()->GetTid(v);
 
            // Check to see if this trace is periodic.
            auto it = periodicTraces.find(GeomId);
 
            if (it != periodicTraces.end())
            {
                const PeriodicEntity &ent = it->second[0];
                auto it2                  = perTraceToExpMap.find(ent.id);
 
                if (it2 == perTraceToExpMap.end())
                {
                    if (m_session->GetComm()->GetRowComm()->GetSize() > 1 &&
                        !ent.isLocal)
                    {
                        continue;
                    }
                    else
                    {
                        ASSERTL1(false, "Periodic trace not found!");
                    }
                }
 
                ASSERTL1(m_leftAdjacentTraces[cnt],
                         "Periodic trace in non-forward space?");
 
                int offset = m_trace->GetPhys_Offset(
                    (m_traceMap->GetElmtToTrace())[n][v]->GetElmtId());
 
                int offset2 = m_trace->GetPhys_Offset(
                    (m_traceMap->GetElmtToTrace())[it2->second.first]
                                                  [it2->second.second]
                                                      ->GetElmtId());
 
                switch (m_expType)
                {
                    case e1D:
                    {
                        m_periodicFwdCopy.push_back(offset);
                        m_periodicBwdCopy.push_back(offset2);
                    }
                    break;
                    case e2D:
                    {
                        // Calculate relative orientations between trace to
                        // calculate copying map.
                        int nquad = elmtToTrace[n][v]->GetNumPoints(0);
 
                        vector<int> tmpBwd(nquad);
                        vector<int> tmpFwd(nquad);
 
                        if (ent.orient == StdRegions::eForwards)
                        {
                            for (int i = 0; i < nquad; ++i)
                            {
                                tmpBwd[i] = offset2 + i;
                                tmpFwd[i] = offset + i;
                            }
                        }
                        else
                        {
                            for (int i = 0; i < nquad; ++i)
                            {
                                tmpBwd[i] = offset2 + i;
                                tmpFwd[i] = offset + nquad - i - 1;
                            }
                        }
 
                        for (int i = 0; i < nquad; ++i)
                        {
                            m_periodicFwdCopy.push_back(tmpFwd[i]);
                            m_periodicBwdCopy.push_back(tmpBwd[i]);
                        }
                    }
                    break;
                    case e3D:
                    {
                        // Calculate relative orientations between faces to
                        // calculate copying map.
                        int nquad1 = elmtToTrace[n][v]->GetNumPoints(0);
                        int nquad2 = elmtToTrace[n][v]->GetNumPoints(1);
 
                        vector<int> tmpBwd(nquad1 * nquad2);
                        vector<int> tmpFwd(nquad1 * nquad2);
 
                        if (ent.orient ==
                                StdRegions::eDir1FwdDir2_Dir2FwdDir1 ||
                            ent.orient ==
                                StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
                            ent.orient ==
                                StdRegions::eDir1FwdDir2_Dir2BwdDir1 ||
                            ent.orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
                        {
                            for (int i = 0; i < nquad2; ++i)
                            {
                                for (int j = 0; j < nquad1; ++j)
                                {
                                    tmpBwd[i * nquad1 + j] =
                                        offset2 + i * nquad1 + j;
                                    tmpFwd[i * nquad1 + j] =
                                        offset + j * nquad2 + i;
                                }
                            }
                        }
                        else
                        {
                            for (int i = 0; i < nquad2; ++i)
                            {
                                for (int j = 0; j < nquad1; ++j)
                                {
                                    tmpBwd[i * nquad1 + j] =
                                        offset2 + i * nquad1 + j;
                                    tmpFwd[i * nquad1 + j] =
                                        offset + i * nquad1 + j;
                                }
                            }
                        }
 
                        if (ent.orient ==
                                StdRegions::eDir1BwdDir1_Dir2FwdDir2 ||
                            ent.orient ==
                                StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
                            ent.orient ==
                                StdRegions::eDir1FwdDir2_Dir2BwdDir1 ||
                            ent.orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
                        {
                            // Reverse x direction
                            for (int i = 0; i < nquad2; ++i)
                            {
                                for (int j = 0; j < nquad1 / 2; ++j)
                                {
                                    swap(tmpFwd[i * nquad1 + j],
                                         tmpFwd[i * nquad1 + nquad1 - j - 1]);
                                }
                            }
                        }
 
                        if (ent.orient ==
                                StdRegions::eDir1FwdDir1_Dir2BwdDir2 ||
                            ent.orient ==
                                StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
                            ent.orient ==
                                StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
                            ent.orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
                        {
                            // Reverse y direction
                            for (int j = 0; j < nquad1; ++j)
                            {
                                for (int i = 0; i < nquad2 / 2; ++i)
                                {
                                    swap(tmpFwd[i * nquad1 + j],
                                         tmpFwd[(nquad2 - i - 1) * nquad1 + j]);
                                }
                            }
                        }
 
                        for (int i = 0; i < nquad1 * nquad2; ++i)
                        {
                            m_periodicFwdCopy.push_back(tmpFwd[i]);
                            m_periodicBwdCopy.push_back(tmpBwd[i]);
                        }
                    }
                    break;
                    default:
                        ASSERTL1(false, "not set up");
                }
            }
        }
    }
 
    m_locTraceToTraceMap = MemoryManager<LocTraceToTraceMap>::AllocateSharedPtr(
        *this, m_trace, elmtToTrace, m_leftAdjacentTraces);
}
 
/**
 * @brief This routine determines if an element is to the "left"
 * of the adjacent trace, which arises from the idea there is a local
 * normal direction between two elements (i.e. on the trace)
 * and one elements would then be the left.
 *
 * This is typically required since we only define one normal
 * direction along the trace which goes from the left adjacent
 * element to the right adjacent element.  It is also useful
 * in DG discretisations to set up a local Riemann problem
 * where we need to have the concept of a local normal
 * direction which is unique between two elements.
 *
 * There are two cases to be checked:
 *
 * 1) First is the trace on a boundary condition (perioidic or
 * otherwise) or on a partition boundary. If a partition
 * boundary then trace is always considered to be the left
 * adjacent trace and the normal is pointing outward of the
 * soltuion domain. We have to perform an additional case for
 * a periodic boundary where wer chose the element with the
 * lowest global id. If the trace is on a parallel partition
 * we use a member of the traceMap where one side is chosen to
 * contribute to a unique map and have a value which is not
 * -1 so this is the identifier for the left adjacent side
 *
 * 2) If the element is a elemental boundary on one element
 * the left adjacent element is defined by a link to the left
 * element from the trace expansion and this is consistent
 * with the defitiion of the normal which is determined by the
 * largest id (in contrast to the periodic boundary definition
 * !). This reversal of convention does not really matter
 * providing the normals are defined consistently.
 */
 
bool DisContField::IsLeftAdjacentTrace(const int n, const int e)
{
    LocalRegions::ExpansionSharedPtr traceEl =
        m_traceMap->GetElmtToTrace()[n][e];
 
    PeriodicMap periodicTraces = (m_expType == e1D)   ? m_periodicVerts
                                 : (m_expType == e2D) ? m_periodicEdges
                                                      : m_periodicFaces;
 
    bool fwd = true;
    if (traceEl->GetLeftAdjacentElementTrace() == -1 ||
        traceEl->GetRightAdjacentElementTrace() == -1)
    {
        // Boundary edge (1 connected element). Do nothing in
        // serial.
        auto it = m_boundaryTraces.find(traceEl->GetElmtId());
 
        // If the edge does not have a boundary condition set on
        // it, then assume it is a partition edge or periodic.
        if (it == m_boundaryTraces.end())
        {
            fwd = true;
        }
    }
    else if (traceEl->GetLeftAdjacentElementTrace() != -1 &&
             traceEl->GetRightAdjacentElementTrace() != -1)
    {
        // Non-boundary edge (2 connected elements).
        fwd = (traceEl->GetLeftAdjacentElementExp().get()) == (*m_exp)[n].get();
    }
    else
    {
        ASSERTL2(false, "Unconnected trace element!");
    }
 
    return fwd;
}
 
// Given all boundary regions for the whole solution determine
// which ones (if any) are part of domain and proivde all other
// conditions are given as UserDefined Dirichlet.
SpatialDomains::BoundaryConditionsSharedPtr DisContField::GetDomainBCs(
    const SpatialDomains::CompositeMap &domain,
    const SpatialDomains::BoundaryConditions &Allbcs,
    [[maybe_unused]] const std::string &variable)
{
    SpatialDomains::BoundaryConditionsSharedPtr returnval;
 
    returnval =
        MemoryManager<SpatialDomains::BoundaryConditions>::AllocateSharedPtr();
 
    map<int, int> GeometryToRegionsMap;
 
    const SpatialDomains::BoundaryRegionCollection &bregions =
        Allbcs.GetBoundaryRegions();
    const SpatialDomains::BoundaryConditionCollection &bconditions =
        Allbcs.GetBoundaryConditions();
 
    // Set up a map of all boundary regions
    for (auto &it : bregions)
    {
        for (auto &bregionIt : *it.second)
        {
            // can assume that all regions only contain one point in 1D
            // Really do not need loop above
            int id = bregionIt.second->m_geomVec[0]->GetGlobalID();
            GeometryToRegionsMap[id] = it.first;
        }
    }
 
    map<int, SpatialDomains::Geometry *> EndOfDomain;
 
    // Now find out which points in domain have only one vertex
    for (auto &domIt : domain)
    {
        SpatialDomains::CompositeSharedPtr geomvector = domIt.second;
        for (int i = 0; i < geomvector->m_geomVec.size(); ++i)
        {
            for (int j = 0; j < 2; ++j)
            {
                int vid = geomvector->m_geomVec[i]->GetVid(j);
                if (EndOfDomain.count(vid) == 0)
                {
                    EndOfDomain[vid] = geomvector->m_geomVec[i]->GetVertex(j);
                }
                else
                {
                    EndOfDomain.erase(vid);
                }
            }
        }
    }
    ASSERTL1(EndOfDomain.size() == 2, "Did not find two ends of domain");
 
    for (auto &regIt : EndOfDomain)
    {
        if (GeometryToRegionsMap.count(regIt.first) != 0)
        {
            // Set up boundary condition up
            auto iter = GeometryToRegionsMap.find(regIt.first);
            ASSERTL1(iter != GeometryToRegionsMap.end(),
                     "Failied to find GeometryToRegionMap");
 
            int regionId      = iter->second;
            auto bregionsIter = bregions.find(regionId);
            ASSERTL1(bregionsIter != bregions.end(),
                     "Failed to find boundary region");
 
            SpatialDomains::BoundaryRegionShPtr breg = bregionsIter->second;
            returnval->AddBoundaryRegions(regionId, breg);
 
            auto bconditionsIter = bconditions.find(regionId);
            ASSERTL1(bconditionsIter != bconditions.end(),
                     "Failed to find boundary collection");
 
            SpatialDomains::BoundaryConditionMapShPtr bcond =
                bconditionsIter->second;
            returnval->AddBoundaryConditions(regionId, bcond);
        }
    }
 
    return returnval;
}
 
/**
 * Constructor for use in multidomain computations where a
 * domain list can be passed instead of graph1D
 *
 * @param   domain  Subdomain specified in the inputfile from
 *                  which the DisContField is set up
 */
DisContField::DisContField(const LibUtilities::SessionReaderSharedPtr &pSession,
                           const SpatialDomains::MeshGraphSharedPtr &graph1D,
                           const SpatialDomains::CompositeMap &domain,
                           const SpatialDomains::BoundaryConditions &Allbcs,
                           const std::string &variable,
                           const LibUtilities::CommSharedPtr &comm,
                           bool SetToOneSpaceDimension,
                           const Collections::ImplementationType ImpType)
    : ExpList(pSession, domain, graph1D, true, variable, SetToOneSpaceDimension,
              comm, ImpType)
{
    if (variable.compare("DefaultVar") != 0)
    {
        SpatialDomains::BoundaryConditionsSharedPtr DomBCs =
            GetDomainBCs(domain, Allbcs, variable);
 
        GenerateBoundaryConditionExpansion(m_graph, *DomBCs, variable, true,
                                           ImpType);
        EvaluateBoundaryConditions(0.0, variable);
        ApplyGeomInfo();
        FindPeriodicTraces(*DomBCs, variable);
    }
 
    SetUpDG(variable);
}
 
/**
 * Constructs a field as a copy of an existing field.
 * @param   In          Existing DisContField object to copy.
 */
DisContField::DisContField(const DisContField &In,
                           const bool DeclareCoeffPhysArrays)
    : ExpList(In, DeclareCoeffPhysArrays), m_bndConditions(In.m_bndConditions),
      m_bndCondExpansions(In.m_bndCondExpansions),
      m_bndCondBndWeight(In.m_bndCondBndWeight),
      m_interfaceMap(In.m_interfaceMap), m_globalBndMat(In.m_globalBndMat),
      m_traceMap(In.m_traceMap), m_boundaryTraces(In.m_boundaryTraces),
      m_periodicVerts(In.m_periodicVerts),
      m_periodicFwdCopy(In.m_periodicFwdCopy),
      m_periodicBwdCopy(In.m_periodicBwdCopy),
      m_leftAdjacentTraces(In.m_leftAdjacentTraces),
      m_locTraceToTraceMap(In.m_locTraceToTraceMap)
{
    if (In.m_trace) // independent trace space
    {
        m_trace = MemoryManager<ExpList>::AllocateSharedPtr(
            *In.m_trace, DeclareCoeffPhysArrays);
    }
 
    if (In.m_locElmtTrace)
    {
        m_locElmtTrace = MemoryManager<ExpList>::AllocateSharedPtr(
            *In.m_locElmtTrace, DeclareCoeffPhysArrays);
    }
}
 
/*
 * @brief Copy type constructor which declares new boundary conditions
 * and re-uses mapping info and trace space if possible
 */
DisContField::DisContField(const DisContField &In,
                           const SpatialDomains::MeshGraphSharedPtr &graph,
                           const std::string &variable, const bool SetUpJustDG,
                           const bool DeclareCoeffPhysArrays)
    : ExpList(In, DeclareCoeffPhysArrays)
{
    m_trace = NullExpListSharedPtr;
 
    // Set up boundary conditions for this variable.
    // Do not set up BCs if default variable
    if (variable.compare("DefaultVar") != 0)
    {
        SpatialDomains::BoundaryConditions bcs(m_session, graph);
        GenerateBoundaryConditionExpansion(In.m_bndCondExpansions, bcs,
                                           variable);
 
        if (DeclareCoeffPhysArrays)
        {
            EvaluateBoundaryConditions(0.0, variable);
        }
 
        if (!SameTypeOfBoundaryConditions(In))
        {
            // Find periodic edges for this variable.
            FindPeriodicTraces(bcs, variable);
 
            if (SetUpJustDG)
            {
                SetUpDG(variable);
            }
            else
            {
                // set elmt edges to point to robin bc edges if required
                int i, cnt = 0;
                Array<OneD, int> ElmtID, TraceID;
                GetBoundaryToElmtMap(ElmtID, TraceID);
 
                for (i = 0; i < m_bndCondExpansions.size(); ++i)
                {
                    MultiRegions::ExpListSharedPtr locExpList;
 
                    int e;
                    locExpList = m_bndCondExpansions[i];
 
                    for (e = 0; e < locExpList->GetExpSize(); ++e)
                    {
                        (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(
                            TraceID[cnt + e], locExpList->GetExp(e));
                        locExpList->GetExp(e)->SetAdjacentElementExp(
                            TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
                    }
 
                    cnt += m_bndCondExpansions[i]->GetExpSize();
                }
 
                if (m_session->DefinesSolverInfo("PROJECTION"))
                {
                    std::string ProjectStr =
                        m_session->GetSolverInfo("PROJECTION");
 
                    if ((ProjectStr == "MixedCGDG") ||
                        (ProjectStr == "Mixed_CG_Discontinuous"))
                    {
                        SetUpDG(variable);
                    }
                    else
                    {
                        SetUpPhysNormals();
                    }
                }
                else
                {
                    SetUpPhysNormals();
                }
            }
        }
        else
        {
            m_globalBndMat       = In.m_globalBndMat;
            m_trace              = In.m_trace;
            m_locElmtTrace       = In.m_locElmtTrace;
            m_traceMap           = In.m_traceMap;
            m_interfaceMap       = In.m_interfaceMap;
            m_locTraceToTraceMap = In.m_locTraceToTraceMap;
            m_periodicVerts      = In.m_periodicVerts;
            m_periodicEdges      = In.m_periodicEdges;
            m_periodicFaces      = In.m_periodicFaces;
            m_periodicFwdCopy    = In.m_periodicFwdCopy;
            m_periodicBwdCopy    = In.m_periodicBwdCopy;
            m_boundaryTraces     = In.m_boundaryTraces;
            m_leftAdjacentTraces = In.m_leftAdjacentTraces;
 
            if (!SetUpJustDG)
            {
                // set elmt edges to point to robin bc edges if required
                int i, cnt = 0;
                Array<OneD, int> ElmtID, TraceID;
                GetBoundaryToElmtMap(ElmtID, TraceID);
 
                for (i = 0; i < m_bndCondExpansions.size(); ++i)
                {
                    MultiRegions::ExpListSharedPtr locExpList;
 
                    int e;
                    locExpList = m_bndCondExpansions[i];
 
                    for (e = 0; e < locExpList->GetExpSize(); ++e)
                    {
                        (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(
                            TraceID[cnt + e], locExpList->GetExp(e));
                        locExpList->GetExp(e)->SetAdjacentElementExp(
                            TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
                    }
                    cnt += m_bndCondExpansions[i]->GetExpSize();
                }
 
                SetUpPhysNormals();
            }
        }
    }
}
 
/**
 * Constructs a field as a copy of an existing explist field.
 * @param   In          Existing ExpList object to copy.
 */
DisContField::DisContField(const ExpList &In) : ExpList(In)
{
}
 
/**
 *
 */
DisContField::~DisContField()
{
}
 
/**
 * \brief This function discretises the boundary conditions by setting
 * up a list of one-dimensions lower boundary expansions.
 *
 * According to their boundary region, the separate  boundary
 * expansions are bundled together in an object of the class
 *
 * @param   graph       A mesh, containing information about the domain
 *                      and the spectral/hp element expansions.
 * @param   bcs         Information about the enforced boundary
 *                      conditions.
 * @param   variable    The session variable associated with the
 *                      boundary conditions to enforce.
 * @param DeclareCoeffPhysArrays bool to identify if array
 *                               space should be setup.
 *                               Default is true.
 */
void DisContField::GenerateBoundaryConditionExpansion(
    const SpatialDomains::MeshGraphSharedPtr &graph,
    const SpatialDomains::BoundaryConditions &bcs, const std::string variable,
    const bool DeclareCoeffPhysArrays,
    const Collections::ImplementationType ImpType)
{
    int cnt = 0;
    SpatialDomains::BoundaryConditionShPtr bc;
    MultiRegions::ExpListSharedPtr locExpList;
    const SpatialDomains::BoundaryRegionCollection &bregions =
        bcs.GetBoundaryRegions();
    const SpatialDomains::BoundaryConditionCollection &bconditions =
        bcs.GetBoundaryConditions();
 
    std::set<int> ProcessBnd;
    if (m_session->DefinesTag("CreateBndRegions"))
    {
        // evaluate bnd regions to be generated
        vector<unsigned int> bndRegions;
        ASSERTL0(ParseUtils::GenerateVector(
                     m_session->GetTag("CreateBndRegions"), bndRegions),
                 "Failed to interpret bnd values string");
 
        for (auto &bnd : bndRegions)
        {
            if (bregions.count(bnd) ==
                1) // this boundary may not exist on processor
            {
                ProcessBnd.insert(bnd);
            }
        }
    }
    else
    {
        for (auto &it : bregions)
        {
            ProcessBnd.insert(it.first);
        }
    }
 
    m_bndCondExpansions =
        Array<OneD, MultiRegions::ExpListSharedPtr>(ProcessBnd.size());
    m_bndConditions =
        Array<OneD, SpatialDomains::BoundaryConditionShPtr>(ProcessBnd.size());
 
    m_bndCondBndWeight = Array<OneD, NekDouble>{bregions.size(), 0.0};
 
    // count the number of non-periodic boundary points
    for (auto &it : bregions)
    {
        // check to see if reduced bnd regions have been set in FieldConvert.
        if (ProcessBnd.count(it.first) == 0)
        {
            continue;
        }
 
        bc = GetBoundaryCondition(bconditions, it.first, variable);
 
        locExpList = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
            m_session, *(it.second), graph, DeclareCoeffPhysArrays, variable,
            false, bc->GetComm(), ImpType);
 
        m_bndCondExpansions[cnt] = locExpList;
        m_bndConditions[cnt]     = bc;
 
        std::string type = m_bndConditions[cnt]->GetUserDefined();
 
        // Set up normals on non-Dirichlet boundary conditions. Second
        // two conditions ideally should be in local solver setup (when
        // made into factory)
        if (bc->GetBoundaryConditionType() != SpatialDomains::eDirichlet ||
            boost::iequals(type, "I") || boost::iequals(type, "CalcBC"))
        {
            SetUpPhysNormals();
        }
        cnt++;
    }
}
 
/**
 * \brief This function discretises the boundary conditions using an already
 * setup expansion list of one-dimensions lower boundary expansions.
 *
 * According to their boundary region, the separate  boundary
 * expansions are bundled together in an object of the class
 *
 * @param   In          An array of boundary conditions from  an alread setup
 *                      expansion and the spectral/hp element expansions.
 * @param   bcs         Information about the enforced boundary
 *                      conditions.
 * @param   variable    The session variable associated with the
 *                      boundary conditions to enforce.
 * @param DeclareCoeffPhysArrays bool to identify if array
 *                               space should be setup.
 *                               Default is true.
 */
void DisContField::GenerateBoundaryConditionExpansion(
    const Array<OneD, const MultiRegions::ExpListSharedPtr> &In,
    const SpatialDomains::BoundaryConditions &bcs, const std::string variable,
    const bool DeclareCoeffPhysArrays,
    [[maybe_unused]] const Collections::ImplementationType ImpType)
{
    int cnt = 0;
    SpatialDomains::BoundaryConditionShPtr bc;
    MultiRegions::ExpListSharedPtr locExpList;
    const SpatialDomains::BoundaryRegionCollection &bregions =
        bcs.GetBoundaryRegions();
    const SpatialDomains::BoundaryConditionCollection &bconditions =
        bcs.GetBoundaryConditions();
 
    std::set<int> ProcessBnd;
    if (m_session->DefinesTag("CreateBndRegions"))
    {
        // evaluate bnd regions to be generated
        vector<unsigned int> bndRegions;
        ASSERTL0(ParseUtils::GenerateVector(
                     m_session->GetTag("CreateBndRegions"), bndRegions),
                 "Failed to interpret bnd values string");
 
        for (auto &bnd : bndRegions)
        {
            if (bregions.count(bnd) == 1)
            {
                ProcessBnd.insert(bnd);
            }
        }
    }
    else
    {
        for (auto &it : bregions)
        {
            ProcessBnd.insert(it.first);
        }
    }
 
    m_bndCondExpansions =
        Array<OneD, MultiRegions::ExpListSharedPtr>(ProcessBnd.size());
    m_bndConditions =
        Array<OneD, SpatialDomains::BoundaryConditionShPtr>(ProcessBnd.size());
 
    m_bndCondBndWeight = Array<OneD, NekDouble>{bregions.size(), 0.0};
 
    // count the number of non-periodic boundary points
    for (auto &it : bregions)
    {
        // check to see if reduced bnd regions have been set in FieldConvert.
        if (ProcessBnd.count(it.first) == 0)
        {
            continue;
        }
 
        bc = GetBoundaryCondition(bconditions, it.first, variable);
 
        // make a copy of existing Bc passed to generate function
        locExpList = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
            *(In[cnt]), DeclareCoeffPhysArrays);
 
        m_bndCondExpansions[cnt] = locExpList;
        m_bndConditions[cnt]     = bc;
 
        std::string type = m_bndConditions[cnt]->GetUserDefined();
 
        // Set up normals on non-Dirichlet boundary conditions. Second
        // two conditions ideally should be in local solver setup (when
        // made into factory)
        if (bc->GetBoundaryConditionType() != SpatialDomains::eDirichlet ||
            boost::iequals(type, "I") || boost::iequals(type, "CalcBC"))
        {
            SetUpPhysNormals();
        }
        cnt++;
    }
}
 
/**
 * @brief Determine the periodic faces, edges and vertices for the given
 * graph.
 *
 * @param   bcs         Information about the boundary conditions.
 * @param   variable    Specifies the field.
 */
void DisContField::FindPeriodicTraces(
    const SpatialDomains::BoundaryConditions &bcs, const std::string variable)
{
    const SpatialDomains::BoundaryRegionCollection &bregions =
        bcs.GetBoundaryRegions();
    const SpatialDomains::BoundaryConditionCollection &bconditions =
        bcs.GetBoundaryConditions();
 
    LibUtilities::CommSharedPtr vComm = m_comm->GetRowComm();
 
    switch (m_expType)
    {
        case e1D:
        {
            int i, region1ID, region2ID;
 
            SpatialDomains::BoundaryConditionShPtr locBCond;
            map<int, int> BregionToVertMap;
 
            // Construct list of all periodic Region and their
            // global vertex on this process.
            for (auto &it : bregions)
            {
                locBCond =
                    GetBoundaryCondition(bconditions, it.first, variable);
 
                if (locBCond->GetBoundaryConditionType() !=
                    SpatialDomains::ePeriodic)
                {
                    continue;
                }
                int id =
                    it.second->begin()->second->m_geomVec[0]->GetGlobalID();
 
                BregionToVertMap[it.first] = id;
            }
 
            set<int> islocal;
 
            int n = vComm->GetSize();
            int p = vComm->GetRank();
 
            Array<OneD, int> nregions(n, 0);
            nregions[p] = BregionToVertMap.size();
            vComm->AllReduce(nregions, LibUtilities::ReduceSum);
 
            int totRegions = Vmath::Vsum(n, nregions, 1);
 
            Array<OneD, int> regOffset(n, 0);
 
            for (i = 1; i < n; ++i)
            {
                regOffset[i] = regOffset[i - 1] + nregions[i - 1];
            }
 
            Array<OneD, int> bregmap(totRegions, 0);
            Array<OneD, int> bregid(totRegions, 0);
 
            i = regOffset[p];
            for (auto &iit : BregionToVertMap)
            {
                bregid[i]    = iit.first;
                bregmap[i++] = iit.second;
                islocal.insert(iit.first);
            }
 
            vComm->AllReduce(bregmap, LibUtilities::ReduceSum);
            vComm->AllReduce(bregid, LibUtilities::ReduceSum);
 
            for (int i = 0; i < totRegions; ++i)
            {
                BregionToVertMap[bregid[i]] = bregmap[i];
            }
 
            // Construct list of all periodic pairs local to this process.
            for (auto &it : bregions)
            {
                locBCond =
                    GetBoundaryCondition(bconditions, it.first, variable);
 
                if (locBCond->GetBoundaryConditionType() !=
                    SpatialDomains::ePeriodic)
                {
                    continue;
                }
 
                // Identify periodic boundary region IDs.
                region1ID = it.first;
                region2ID =
                    std::static_pointer_cast<
                        SpatialDomains::PeriodicBoundaryCondition>(locBCond)
                        ->m_connectedBoundaryRegion;
 
                ASSERTL0(BregionToVertMap.count(region1ID) != 0,
                         "Cannot determine vertex of region1ID");
 
                ASSERTL0(BregionToVertMap.count(region2ID) != 0,
                         "Cannot determine vertex of region2ID");
 
                PeriodicEntity ent(BregionToVertMap[region2ID],
                                   StdRegions::eNoOrientation,
                                   islocal.count(region2ID) != 0);
 
                m_periodicVerts[BregionToVertMap[region1ID]].push_back(ent);
            }
        }
        break;
        case e2D:
        {
            int region1ID, region2ID, i, j, k, cnt;
            SpatialDomains::BoundaryConditionShPtr locBCond;
 
            SpatialDomains::CompositeOrdering compOrder =
                m_graph->GetCompositeOrdering();
            SpatialDomains::BndRegionOrdering bndRegOrder =
                m_graph->GetBndRegionOrdering();
            SpatialDomains::CompositeMap compMap = m_graph->GetComposites();
 
            // Unique collection of pairs of periodic composites
            // (i.e. if composites 1 and 2 are periodic then this
            // map will contain either the pair (1,2) or (2,1) but
            // not both).
            map<int, int> perComps;
            map<int, vector<int>> allVerts;
            set<int> locVerts;
            map<int, pair<int, StdRegions::Orientation>> allEdges;
 
            // Set up a set of all local verts and edges.
            for (i = 0; i < (*m_exp).size(); ++i)
            {
                for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
                {
                    int id = (*m_exp)[i]->GetGeom()->GetVid(j);
                    locVerts.insert(id);
                }
            }
 
            // Construct list of all periodic pairs local to this process.
            for (auto &it : bregions)
            {
                locBCond =
                    GetBoundaryCondition(bconditions, it.first, variable);
 
                if (locBCond->GetBoundaryConditionType() !=
                    SpatialDomains::ePeriodic)
                {
                    continue;
                }
 
                // Identify periodic boundary region IDs.
                region1ID = it.first;
                region2ID =
                    std::static_pointer_cast<
                        SpatialDomains::PeriodicBoundaryCondition>(locBCond)
                        ->m_connectedBoundaryRegion;
 
                // From this identify composites. Note that in
                // serial this will be an empty map.
                int cId1, cId2;
                if (vComm->GetSize() == 1)
                {
                    cId1 = it.second->begin()->first;
                    cId2 = bregions.find(region2ID)->second->begin()->first;
                }
                else
                {
                    cId1 = bndRegOrder.find(region1ID)->second[0];
                    cId2 = bndRegOrder.find(region2ID)->second[0];
                }
 
                ASSERTL0(it.second->size() == 1,
                         "Boundary region " + std::to_string(region1ID) +
                             " should only contain 1 composite.");
 
                vector<unsigned int> tmpOrder;
 
                // Construct set containing all periodic edgesd on
                // this process
                SpatialDomains::CompositeSharedPtr c =
                    it.second->begin()->second;
 
                for (i = 0; i < c->m_geomVec.size(); ++i)
                {
                    SpatialDomains::SegGeom *segGeom =
                        dynamic_cast<SpatialDomains::SegGeom *>(
                            c->m_geomVec[i]);
                    ASSERTL0(segGeom, "Unable to cast to SegGeom*");
 
                    SpatialDomains::GeometryLinkSharedPtr elmt =
                        m_graph->GetElementsFromEdge(segGeom);
                    ASSERTL0(elmt->size() == 1,
                             "The periodic boundaries belong to "
                             "more than one element of the mesh");
 
                    SpatialDomains::Geometry2D *geom =
                        dynamic_cast<SpatialDomains::Geometry2D *>(
                            elmt->at(0).first);
 
                    allEdges[c->m_geomVec[i]->GetGlobalID()] =
                        make_pair(elmt->at(0).second,
                                  geom->GetEorient(elmt->at(0).second));
 
                    // In serial mesh partitioning will not have occurred so
                    // need to fill composite ordering map manually.
                    if (vComm->GetSize() == 1)
                    {
                        tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
                    }
 
                    vector<int> vertList(2);
                    vertList[0] = segGeom->GetVid(0);
                    vertList[1] = segGeom->GetVid(1);
                    allVerts[c->m_geomVec[i]->GetGlobalID()] = vertList;
                }
 
                if (vComm->GetSize() == 1)
                {
                    compOrder[it.second->begin()->first] = tmpOrder;
                }
 
                // See if we already have either region1 or
                // region2 stored in perComps map.
                if (perComps.count(cId1) == 0)
                {
                    if (perComps.count(cId2) == 0)
                    {
                        perComps[cId1] = cId2;
                    }
                    else
                    {
                        std::stringstream ss;
                        ss << "Boundary region " << cId2 << " should be "
                           << "periodic with " << perComps[cId2] << " but "
                           << "found " << cId1 << " instead!";
                        ASSERTL0(perComps[cId2] == cId1, ss.str());
                    }
                }
                else
                {
                    std::stringstream ss;
                    ss << "Boundary region " << cId1 << " should be "
                       << "periodic with " << perComps[cId1] << " but "
                       << "found " << cId2 << " instead!";
                    ASSERTL0(perComps[cId1] == cId1, ss.str());
                }
            }
 
            // In case of periodic partition being split over many composites
            // may not have all periodic matches so check this
            int idmax = -1;
            for (auto &cIt : perComps)
            {
                idmax = max(idmax, cIt.first);
                idmax = max(idmax, cIt.second);
            }
            vComm->AllReduce(idmax, LibUtilities::ReduceMax);
            idmax++;
            Array<OneD, int> perid(idmax, -1);
            for (auto &cIt : perComps)
            {
                perid[cIt.first] = cIt.second;
            }
            vComm->AllReduce(perid, LibUtilities::ReduceMax);
            // update all partitions
            for (int i = 0; i < idmax; ++i)
            {
                if (perid[i] > -1)
                {
                    // skip if complemetary relationship has
                    // already been speficied in list
                    if (perComps.count(perid[i]))
                    {
                        continue;
                    }
                    perComps[i] = perid[i];
                }
            }
 
            // Process local edge list to obtain relative edge
            // orientations.
            int n = vComm->GetSize();
            int p = vComm->GetRank();
            int totEdges;
            Array<OneD, int> edgecounts(n, 0);
            Array<OneD, int> edgeoffset(n, 0);
            Array<OneD, int> vertoffset(n, 0);
 
            edgecounts[p] = allEdges.size();
            vComm->AllReduce(edgecounts, LibUtilities::ReduceSum);
 
            edgeoffset[0] = 0;
            for (i = 1; i < n; ++i)
            {
                edgeoffset[i] = edgeoffset[i - 1] + edgecounts[i - 1];
            }
 
            totEdges = Vmath::Vsum(n, edgecounts, 1);
            Array<OneD, int> edgeIds(totEdges, 0);
            Array<OneD, int> edgeIdx(totEdges, 0);
            Array<OneD, int> edgeOrient(totEdges, 0);
            Array<OneD, int> edgeVerts(totEdges, 0);
 
            auto sIt = allEdges.begin();
 
            for (i = 0; sIt != allEdges.end(); ++sIt)
            {
                edgeIds[edgeoffset[p] + i]     = sIt->first;
                edgeIdx[edgeoffset[p] + i]     = sIt->second.first;
                edgeOrient[edgeoffset[p] + i]  = sIt->second.second;
                edgeVerts[edgeoffset[p] + i++] = allVerts[sIt->first].size();
            }
 
            vComm->AllReduce(edgeIds, LibUtilities::ReduceSum);
            vComm->AllReduce(edgeIdx, LibUtilities::ReduceSum);
            vComm->AllReduce(edgeOrient, LibUtilities::ReduceSum);
            vComm->AllReduce(edgeVerts, LibUtilities::ReduceSum);
 
            // Calculate number of vertices on each processor.
            Array<OneD, int> procVerts(n, 0);
            int nTotVerts;
 
            // Note if there are no periodic edges at all calling Vsum will
            // cause a segfault.
            if (totEdges > 0)
            {
                nTotVerts = Vmath::Vsum(totEdges, edgeVerts, 1);
            }
            else
            {
                nTotVerts = 0;
            }
 
            for (i = 0; i < n; ++i)
            {
                if (edgecounts[i] > 0)
                {
                    procVerts[i] = Vmath::Vsum(edgecounts[i],
                                               edgeVerts + edgeoffset[i], 1);
                }
                else
                {
                    procVerts[i] = 0;
                }
            }
            vertoffset[0] = 0;
 
            for (i = 1; i < n; ++i)
            {
                vertoffset[i] = vertoffset[i - 1] + procVerts[i - 1];
            }
 
            Array<OneD, int> traceIds(nTotVerts, 0);
            for (i = 0, sIt = allEdges.begin(); sIt != allEdges.end(); ++sIt)
            {
                for (j = 0; j < allVerts[sIt->first].size(); ++j)
                {
                    traceIds[vertoffset[p] + i++] = allVerts[sIt->first][j];
                }
            }
 
            vComm->AllReduce(traceIds, LibUtilities::ReduceSum);
 
            // For simplicity's sake create a map of edge id -> orientation.
            map<int, StdRegions::Orientation> orientMap;
            map<int, vector<int>> vertMap;
 
            for (cnt = i = 0; i < totEdges; ++i)
            {
                ASSERTL0(orientMap.count(edgeIds[i]) == 0,
                         "Already found edge in orientation map!");
 
                // Work out relative orientations. To avoid having
                // to exchange vertex locations, we figure out if
                // the edges are backwards or forwards orientated
                // with respect to a forwards orientation that is
                // CCW. Since our local geometries are
                // forwards-orientated with respect to the
                // Cartesian axes, we need to invert the
                // orientation for the top and left edges of a
                // quad and the left edge of a triangle.
                StdRegions::Orientation o =
                    (StdRegions::Orientation)edgeOrient[i];
 
                if (edgeIdx[i] > 1)
                {
                    o = o == StdRegions::eForwards ? StdRegions::eBackwards
                                                   : StdRegions::eForwards;
                }
 
                orientMap[edgeIds[i]] = o;
 
                vector<int> verts(edgeVerts[i]);
 
                for (j = 0; j < edgeVerts[i]; ++j)
                {
                    verts[j] = traceIds[cnt++];
                }
                vertMap[edgeIds[i]] = verts;
            }
 
            // Go through list of composites and figure out which
            // edges are parallel from original ordering in
            // session file. This includes composites which are
            // not necessarily on this process.
            map<int, int> allCompPairs;
 
            // Store temporary map of periodic vertices which will hold all
            // periodic vertices on the entire mesh so that doubly periodic
            // vertices can be counted properly across partitions. Local
            // vertices are copied into m_periodicVerts at the end of the
            // function.
            PeriodicMap periodicVerts;
 
            for (auto &cIt : perComps)
            {
                SpatialDomains::CompositeSharedPtr c[2];
                const int id1    = cIt.first;
                const int id2    = cIt.second;
                std::string id1s = std::to_string(id1);
                std::string id2s = std::to_string(id2);
 
                if (compMap.count(id1) > 0)
                {
                    c[0] = compMap[id1];
                }
 
                if (compMap.count(id2) > 0)
                {
                    c[1] = compMap[id2];
                }
 
                // Loop over composite ordering to construct list of all
                // periodic edges regardless of whether they are on this
                // process.
                map<int, int> compPairs;
 
                ASSERTL0(compOrder.count(id1) > 0,
                         "Unable to find composite " + id1s + " in order map.");
                ASSERTL0(compOrder.count(id2) > 0,
                         "Unable to find composite " + id2s + " in order map.");
                WARNINGL0(
                    compOrder[id1].size() == compOrder[id2].size(),
                    "Periodic composites " + id1s + " and " + id2s +
                        " should have the same number of elements. Have " +
                        std::to_string(compOrder[id1].size()) + " vs " +
                        std::to_string(compOrder[id2].size()));
                ASSERTL0(compOrder[id1].size() > 0, "Periodic composites " +
                                                        id1s + " and " + id2s +
                                                        " are empty!");
 
                // TODO: Add more checks.
                for (i = 0; i < compOrder[id1].size(); ++i)
                {
                    int eId1 = compOrder[id1][i];
                    int eId2 = compOrder[id2][i];
 
                    ASSERTL0(compPairs.count(eId1) == 0, "Already paired.");
 
                    if (compPairs.count(eId2) != 0)
                    {
                        ASSERTL0(compPairs[eId2] == eId1, "Pairing incorrect");
                    }
                    compPairs[eId1] = eId2;
                }
 
                // Construct set of all edges that we have
                // locally on this processor.
                set<int> locEdges;
                for (i = 0; i < 2; ++i)
                {
                    if (!c[i])
                    {
                        continue;
                    }
 
                    if (c[i]->m_geomVec.size() > 0)
                    {
                        for (j = 0; j < c[i]->m_geomVec.size(); ++j)
                        {
                            locEdges.insert(c[i]->m_geomVec[j]->GetGlobalID());
                        }
                    }
                }
 
                // Loop over all edges in the geometry composite.
                for (auto &pIt : compPairs)
                {
                    int ids[2]    = {pIt.first, pIt.second};
                    bool local[2] = {locEdges.count(pIt.first) > 0,
                                     locEdges.count(pIt.second) > 0};
 
                    ASSERTL0(orientMap.count(ids[0]) > 0 &&
                                 orientMap.count(ids[1]) > 0,
                             "Unable to find edge in orientation map.");
 
                    allCompPairs[pIt.first]  = pIt.second;
                    allCompPairs[pIt.second] = pIt.first;
 
                    for (i = 0; i < 2; ++i)
                    {
                        if (!local[i])
                        {
                            continue;
                        }
 
                        int other = (i + 1) % 2;
 
                        StdRegions::Orientation o =
                            orientMap[ids[i]] == orientMap[ids[other]]
                                ? StdRegions::eBackwards
                                : StdRegions::eForwards;
 
                        PeriodicEntity ent(ids[other], o, local[other]);
                        m_periodicEdges[ids[i]].push_back(ent);
                    }
 
                    for (i = 0; i < 2; ++i)
                    {
                        int other = (i + 1) % 2;
 
                        StdRegions::Orientation o =
                            orientMap[ids[i]] == orientMap[ids[other]]
                                ? StdRegions::eBackwards
                                : StdRegions::eForwards;
 
                        // Determine periodic vertices.
                        vector<int> perVerts1 = vertMap[ids[i]];
                        vector<int> perVerts2 = vertMap[ids[other]];
 
                        map<int, pair<int, bool>> tmpMap;
                        if (o == StdRegions::eForwards)
                        {
                            tmpMap[perVerts1[0]] = make_pair(
                                perVerts2[0], locVerts.count(perVerts2[0]) > 0);
                            tmpMap[perVerts1[1]] = make_pair(
                                perVerts2[1], locVerts.count(perVerts2[1]) > 0);
                        }
                        else
                        {
                            tmpMap[perVerts1[0]] = make_pair(
                                perVerts2[1], locVerts.count(perVerts2[1]) > 0);
                            tmpMap[perVerts1[1]] = make_pair(
                                perVerts2[0], locVerts.count(perVerts2[0]) > 0);
                        }
 
                        for (auto &mIt : tmpMap)
                        {
                            // See if this vertex has been recorded already.
                            PeriodicEntity ent2(mIt.second.first,
                                                StdRegions::eNoOrientation,
                                                mIt.second.second);
                            auto perIt = periodicVerts.find(mIt.first);
 
                            if (perIt == periodicVerts.end())
                            {
                                periodicVerts[mIt.first].push_back(ent2);
                                perIt = periodicVerts.find(mIt.first);
                            }
                            else
                            {
                                bool doAdd = true;
                                for (j = 0; j < perIt->second.size(); ++j)
                                {
                                    if (perIt->second[j].id == mIt.second.first)
                                    {
                                        doAdd = false;
                                        break;
                                    }
                                }
 
                                if (doAdd)
                                {
                                    perIt->second.push_back(ent2);
                                }
                            }
                        }
                    }
                }
            }
 
            // Search for periodic vertices and edges which are not in
            // a periodic composite but lie in this process. First, loop
            // over all information we have from other processors.
            for (cnt = i = 0; i < totEdges; ++i)
            {
                int edgeId = edgeIds[i];
 
                ASSERTL0(allCompPairs.count(edgeId) > 0,
                         "Unable to find matching periodic edge.");
 
                int perEdgeId = allCompPairs[edgeId];
 
                for (j = 0; j < edgeVerts[i]; ++j, ++cnt)
                {
                    int vId = traceIds[cnt];
 
                    auto perId = periodicVerts.find(vId);
 
                    if (perId == periodicVerts.end())
                    {
                        // This vertex is not included in the
                        // map. Figure out which vertex it is
                        // supposed to be periodic with. perEdgeId
                        // is the edge ID which is periodic with
                        // edgeId. The logic is much the same as
                        // the loop above.
                        int perVertexId =
                            orientMap[edgeId] == orientMap[perEdgeId]
                                ? vertMap[perEdgeId][(j + 1) % 2]
                                : vertMap[perEdgeId][j];
 
                        PeriodicEntity ent(perVertexId,
                                           StdRegions::eNoOrientation,
                                           locVerts.count(perVertexId) > 0);
 
                        periodicVerts[vId].push_back(ent);
                    }
                }
            }
 
            // Loop over all periodic vertices to complete connectivity
            // information.
            for (auto &perIt : periodicVerts)
            {
                // Loop over associated vertices.
                for (i = 0; i < perIt.second.size(); ++i)
                {
                    auto perIt2 = periodicVerts.find(perIt.second[i].id);
                    ASSERTL0(perIt2 != periodicVerts.end(),
                             "Couldn't find periodic vertex.");
 
                    for (j = 0; j < perIt2->second.size(); ++j)
                    {
                        if (perIt2->second[j].id == perIt.first)
                        {
                            continue;
                        }
 
                        bool doAdd = true;
                        for (k = 0; k < perIt.second.size(); ++k)
                        {
                            if (perIt2->second[j].id == perIt.second[k].id)
                            {
                                doAdd = false;
                                break;
                            }
                        }
 
                        if (doAdd)
                        {
                            perIt.second.push_back(perIt2->second[j]);
                        }
                    }
                }
            }
 
            // Do one final loop over periodic vertices to remove non-local
            // vertices from map.
            for (auto &perIt : periodicVerts)
            {
                if (locVerts.count(perIt.first) > 0)
                {
                    m_periodicVerts.insert(perIt);
                }
            }
        }
        break;
        case e3D:
        {
            SpatialDomains::CompositeOrdering compOrder =
                m_graph->GetCompositeOrdering();
            SpatialDomains::BndRegionOrdering bndRegOrder =
                m_graph->GetBndRegionOrdering();
            SpatialDomains::CompositeMap compMap = m_graph->GetComposites();
 
            // perComps: Stores a unique collection of pairs of periodic
            // composites (i.e. if composites 1 and 2 are periodic then this map
            // will contain either the pair (1,2) or (2,1) but not both).
            //
            // The four maps allVerts, allCoord, allEdges and allOrient map a
            // periodic face to a vector containing the vertex ids of the face;
            // their coordinates; the edge ids of the face; and their
            // orientation within that face respectively.
            //
            // Finally the three sets locVerts, locEdges and locFaces store any
            // vertices, edges and faces that belong to a periodic composite and
            // lie on this process.
            map<int, RotPeriodicInfo> rotComp;
            map<int, int> perComps;
            map<int, vector<int>> allVerts;
            map<int, std::vector<SpatialDomains::PointGeom *>> allCoord;
            map<int, vector<int>> allEdges;
            map<int, vector<StdRegions::Orientation>> allOrient;
            set<int> locVerts;
            set<int> locEdges;
            set<int> locFaces;
 
            int region1ID, region2ID, i, j, k, cnt;
            SpatialDomains::BoundaryConditionShPtr locBCond;
 
            // Set up a set of all local verts and edges.
            for (i = 0; i < (*m_exp).size(); ++i)
            {
                for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
                {
                    int id = (*m_exp)[i]->GetGeom()->GetVid(j);
                    locVerts.insert(id);
                }
 
                for (j = 0; j < (*m_exp)[i]->GetGeom()->GetNumEdges(); ++j)
                {
                    int id = (*m_exp)[i]->GetGeom()->GetEid(j);
                    locEdges.insert(id);
                }
            }
 
            // Begin by populating the perComps map. We loop over all periodic
            // boundary conditions and determine the composite associated with
            // it, then fill out the all* maps.
            for (auto &it : bregions)
            {
 
                locBCond =
                    GetBoundaryCondition(bconditions, it.first, variable);
 
                if (locBCond->GetBoundaryConditionType() !=
                    SpatialDomains::ePeriodic)
                {
                    continue;
                }
 
                // Identify periodic boundary region IDs.
                region1ID = it.first;
                region2ID =
                    std::static_pointer_cast<
                        SpatialDomains::PeriodicBoundaryCondition>(locBCond)
                        ->m_connectedBoundaryRegion;
 
                // Check the region only contains a single composite.
                ASSERTL0(it.second->size() == 1,
                         "Boundary region " + std::to_string(region1ID) +
                             " should only contain 1 composite.");
 
                // From this identify composites by looking at the original
                // boundary region ordering. Note that in serial the mesh
                // partitioner is not run, so this map will be empty and
                // therefore needs to be populated by using the corresponding
                // boundary region.
                int cId1, cId2;
                if (vComm->GetSize() == 1)
                {
                    cId1 = it.second->begin()->first;
                    cId2 = bregions.find(region2ID)->second->begin()->first;
                }
                else
                {
                    cId1 = bndRegOrder.find(region1ID)->second[0];
                    cId2 = bndRegOrder.find(region2ID)->second[0];
                }
 
                // check to see if boundary is rotationally aligned
                if (boost::icontains(locBCond->GetUserDefined(), "Rotated"))
                {
                    vector<string> tmpstr;
 
                    boost::split(tmpstr, locBCond->GetUserDefined(),
                                 boost::is_any_of(":"));
 
                    if (boost::iequals(tmpstr[0], "Rotated"))
                    {
                        ASSERTL1(tmpstr.size() > 2,
                                 "Expected Rotated user defined string to "
                                 "contain direction and rotation angle "
                                 "and optionally a tolerance, "
                                 "i.e. Rotated:dir:PI/2:1e-6");
 
                        ASSERTL1((tmpstr[1] == "x") || (tmpstr[1] == "y") ||
                                     (tmpstr[1] == "z"),
                                 "Rotated Dir is "
                                 "not specified as x,y or z");
 
                        RotPeriodicInfo RotInfo;
                        RotInfo.m_dir = (tmpstr[1] == "x")   ? 0
                                        : (tmpstr[1] == "y") ? 1
                                                             : 2;
 
                        LibUtilities::Interpreter strEval;
                        int ExprId      = strEval.DefineFunction("", tmpstr[2]);
                        RotInfo.m_angle = strEval.Evaluate(ExprId);
 
                        if (tmpstr.size() == 4)
                        {
                            try
                            {
                                RotInfo.m_tol = std::stod(tmpstr[3]);
                            }
                            catch (...)
                            {
                                NEKERROR(ErrorUtil::efatal,
                                         "failed to cast tolerance input "
                                         "to a double value in Rotated"
                                         "boundary information");
                            }
                        }
                        else
                        {
                            RotInfo.m_tol = 1e-8;
                        }
                        rotComp[cId1] = RotInfo;
                    }
                }
 
                SpatialDomains::CompositeSharedPtr c =
                    it.second->begin()->second;
 
                vector<unsigned int> tmpOrder;
 
                // store the rotation info of this
 
                // From the composite, we now construct the allVerts, allEdges
                // and allCoord map so that they can be transferred across
                // processors. We also populate the locFaces set to store a
                // record of all faces local to this process.
                for (i = 0; i < c->m_geomVec.size(); ++i)
                {
                    SpatialDomains::Geometry2D *faceGeom =
                        dynamic_cast<SpatialDomains::Geometry2D *>(
                            c->m_geomVec[i]);
                    ASSERTL1(faceGeom, "Unable to cast to Geometry2D*");
 
                    // Get geometry ID of this face and store in locFaces.
                    int faceId = c->m_geomVec[i]->GetGlobalID();
                    locFaces.insert(faceId);
 
                    // In serial, mesh partitioning will not have occurred so
                    // need to fill composite ordering map manually.
                    if (vComm->GetSize() == 1)
                    {
                        tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
                    }
 
                    // Loop over vertices and edges of the face to populate
                    // allVerts, allEdges and allCoord maps.
                    vector<int> vertList, edgeList;
                    std::vector<SpatialDomains::PointGeom *> coordVec;
                    vector<StdRegions::Orientation> orientVec;
                    for (j = 0; j < faceGeom->GetNumVerts(); ++j)
                    {
                        vertList.push_back(faceGeom->GetVid(j));
                        edgeList.push_back(faceGeom->GetEid(j));
                        coordVec.push_back(faceGeom->GetVertex(j));
                        orientVec.push_back(faceGeom->GetEorient(j));
                    }
 
                    allVerts[faceId]  = vertList;
                    allEdges[faceId]  = edgeList;
                    allCoord[faceId]  = coordVec;
                    allOrient[faceId] = orientVec;
                }
 
                // In serial, record the composite ordering in compOrder for
                // later in the routine.
                if (vComm->GetSize() == 1)
                {
                    compOrder[it.second->begin()->first] = tmpOrder;
                }
 
                // See if we already have either region1 or region2 stored in
                // perComps map already and do a sanity check to ensure regions
                // are mutually periodic.
                if (perComps.count(cId1) == 0)
                {
                    if (perComps.count(cId2) == 0)
                    {
                        perComps[cId1] = cId2;
                    }
                    else
                    {
                        std::stringstream ss;
                        ss << "Boundary region " << cId2 << " should be "
                           << "periodic with " << perComps[cId2] << " but "
                           << "found " << cId1 << " instead!";
                        ASSERTL0(perComps[cId2] == cId1, ss.str());
                    }
                }
                else
                {
                    std::stringstream ss;
                    ss << "Boundary region " << cId1 << " should be "
                       << "periodic with " << perComps[cId1] << " but "
                       << "found " << cId2 << " instead!";
                    ASSERTL0(perComps[cId1] == cId1, ss.str());
                }
            }
 
            // In case of periodic partition being split over many composites
            // may not have all periodic matches so check this
            int idmax = -1;
            for (auto &cIt : perComps)
            {
                idmax = max(idmax, cIt.first);
                idmax = max(idmax, cIt.second);
            }
            vComm->AllReduce(idmax, LibUtilities::ReduceMax);
            idmax++;
            Array<OneD, int> perid(idmax, -1);
            for (auto &cIt : perComps)
            {
                perid[cIt.first] = cIt.second;
            }
            vComm->AllReduce(perid, LibUtilities::ReduceMax);
            // update all partitions
            for (int i = 0; i < idmax; ++i)
            {
                if (perid[i] > -1)
                {
                    // skip if equivlaent relationship has
                    // already been speficied in list
                    if (perComps.count(perid[i]))
                    {
                        continue;
                    }
                    perComps[i] = perid[i];
                }
            }
 
            // The next routines process local face lists to
            // exchange vertices,
            // edges and faces.
            int n = vComm->GetSize();
            int p = vComm->GetRank();
            int totFaces;
            Array<OneD, int> facecounts(n, 0);
            Array<OneD, int> vertcounts(n, 0);
            Array<OneD, int> faceoffset(n, 0);
            Array<OneD, int> vertoffset(n, 0);
 
            Array<OneD, int> rotcounts(n, 0);
            Array<OneD, int> rotoffset(n, 0);
 
            rotcounts[p] = rotComp.size();
            vComm->AllReduce(rotcounts, LibUtilities::ReduceSum);
            int totrot = Vmath::Vsum(n, rotcounts, 1);
 
            if (totrot)
            {
                for (i = 1; i < n; ++i)
                {
                    rotoffset[i] = rotoffset[i - 1] + rotcounts[i - 1];
                }
 
                Array<OneD, int> compid(totrot, 0);
                Array<OneD, int> rotdir(totrot, 0);
                Array<OneD, NekDouble> rotangle(totrot, 0.0);
                Array<OneD, NekDouble> rottol(totrot, 0.0);
 
                // fill in rotational informaiton
                auto rIt = rotComp.begin();
 
                for (i = 0; rIt != rotComp.end(); ++rIt)
                {
                    compid[rotoffset[p] + i]   = rIt->first;
                    rotdir[rotoffset[p] + i]   = rIt->second.m_dir;
                    rotangle[rotoffset[p] + i] = rIt->second.m_angle;
                    rottol[rotoffset[p] + i++] = rIt->second.m_tol;
                }
 
                vComm->AllReduce(compid, LibUtilities::ReduceSum);
                vComm->AllReduce(rotdir, LibUtilities::ReduceSum);
                vComm->AllReduce(rotangle, LibUtilities::ReduceSum);
                vComm->AllReduce(rottol, LibUtilities::ReduceSum);
 
                // Fill in full rotational composite list
                for (i = 0; i < totrot; ++i)
                {
                    RotPeriodicInfo rinfo(rotdir[i], rotangle[i], rottol[i]);
 
                    rotComp[compid[i]] = rinfo;
                }
            }
 
            // First exchange the number of faces on each process.
            facecounts[p] = locFaces.size();
            vComm->AllReduce(facecounts, LibUtilities::ReduceSum);
 
            // Set up an offset map to allow us to distribute face IDs to all
            // processors.
            faceoffset[0] = 0;
            for (i = 1; i < n; ++i)
            {
                faceoffset[i] = faceoffset[i - 1] + facecounts[i - 1];
            }
 
            // Calculate total number of faces.
            totFaces = Vmath::Vsum(n, facecounts, 1);
 
            // faceIds holds face IDs for each periodic face. faceVerts holds
            // the number of vertices in this face.
            Array<OneD, int> faceIds(totFaces, 0);
            Array<OneD, int> faceVerts(totFaces, 0);
 
            // Process p writes IDs of its faces into position faceoffset[p] of
            // faceIds which allows us to perform an AllReduce to distribute
            // information amongst processors.
            auto sIt = locFaces.begin();
            for (i = 0; sIt != locFaces.end(); ++sIt)
            {
                faceIds[faceoffset[p] + i]     = *sIt;
                faceVerts[faceoffset[p] + i++] = allVerts[*sIt].size();
            }
 
            vComm->AllReduce(faceIds, LibUtilities::ReduceSum);
            vComm->AllReduce(faceVerts, LibUtilities::ReduceSum);
 
            // procVerts holds number of vertices (and also edges since each
            // face is 2D) on each process.
            Array<OneD, int> procVerts(n, 0);
            int nTotVerts;
 
            // Note if there are no periodic faces at all calling Vsum will
            // cause a segfault.
            if (totFaces > 0)
            {
                // Calculate number of vertices on each processor.
                nTotVerts = Vmath::Vsum(totFaces, faceVerts, 1);
            }
            else
            {
                nTotVerts = 0;
            }
 
            for (i = 0; i < n; ++i)
            {
                if (facecounts[i] > 0)
                {
                    procVerts[i] = Vmath::Vsum(facecounts[i],
                                               faceVerts + faceoffset[i], 1);
                }
                else
                {
                    procVerts[i] = 0;
                }
            }
 
            // vertoffset is defined in the same manner as edgeoffset
            // beforehand.
            vertoffset[0] = 0;
            for (i = 1; i < n; ++i)
            {
                vertoffset[i] = vertoffset[i - 1] + procVerts[i - 1];
            }
 
            // At this point we exchange all vertex IDs, edge IDs and vertex
            // coordinates for each face. The coordinates are necessary because
            // we need to calculate relative face orientations between periodic
            // faces to determined edge and vertex connectivity.
            Array<OneD, int> vertIds(nTotVerts, 0);
            Array<OneD, int> edgeIds(nTotVerts, 0);
            Array<OneD, int> edgeOrt(nTotVerts, 0);
            Array<OneD, NekDouble> vertX(nTotVerts, 0.0);
            Array<OneD, NekDouble> vertY(nTotVerts, 0.0);
            Array<OneD, NekDouble> vertZ(nTotVerts, 0.0);
 
            for (cnt = 0, sIt = locFaces.begin(); sIt != locFaces.end(); ++sIt)
            {
                for (j = 0; j < allVerts[*sIt].size(); ++j)
                {
                    int vertId                     = allVerts[*sIt][j];
                    vertIds[vertoffset[p] + cnt]   = vertId;
                    vertX[vertoffset[p] + cnt]     = (*allCoord[*sIt][j])(0);
                    vertY[vertoffset[p] + cnt]     = (*allCoord[*sIt][j])(1);
                    vertZ[vertoffset[p] + cnt]     = (*allCoord[*sIt][j])(2);
                    edgeIds[vertoffset[p] + cnt]   = allEdges[*sIt][j];
                    edgeOrt[vertoffset[p] + cnt++] = allOrient[*sIt][j];
                }
            }
 
            vComm->AllReduce(vertIds, LibUtilities::ReduceSum);
            vComm->AllReduce(vertX, LibUtilities::ReduceSum);
            vComm->AllReduce(vertY, LibUtilities::ReduceSum);
            vComm->AllReduce(vertZ, LibUtilities::ReduceSum);
            vComm->AllReduce(edgeIds, LibUtilities::ReduceSum);
            vComm->AllReduce(edgeOrt, LibUtilities::ReduceSum);
 
            // Finally now we have all of this information, we construct maps
            // which make accessing the information easier. These are
            // conceptually the same as all* maps at the beginning of the
            // routine, but now hold information for all periodic vertices.
            map<int, vector<int>> vertMap;
            map<int, vector<int>> edgeMap;
            map<int, std::vector<SpatialDomains::PointGeom *>> coordMap;
 
            // These final two maps are required for determining the relative
            // orientation of periodic edges. vCoMap associates vertex IDs with
            // their coordinates, and eIdMap maps an edge ID to the two vertices
            // which construct it.
            map<int, SpatialDomains::PointGeomUniquePtr> vCoMap;
            map<int, pair<int, int>> eIdMap;
 
            for (cnt = i = 0; i < totFaces; ++i)
            {
                vector<int> edges(faceVerts[i]);
                vector<int> verts(faceVerts[i]);
                std::vector<SpatialDomains::PointGeom *> coord(faceVerts[i]);
 
                // Keep track of cnt to enable correct edge vertices to be
                // inserted into eIdMap.
                int tmp = cnt;
                for (j = 0; j < faceVerts[i]; ++j, ++cnt)
                {
                    edges[j] = edgeIds[cnt];
                    verts[j] = vertIds[cnt];
 
                    auto vIt = vCoMap.find(vertIds[cnt]);
 
                    if (vIt == vCoMap.end())
                    {
                        auto pt = ObjPoolManager<SpatialDomains::PointGeom>::
                            AllocateUniquePtr(3, verts[j], vertX[cnt],
                                              vertY[cnt], vertZ[cnt]);
                        coord[j]             = pt.get();
                        vCoMap[vertIds[cnt]] = std::move(pt);
                    }
                    else
                    {
                        coord[j] = vIt->second.get();
                    }
 
                    // Try to insert edge into the eIdMap to avoid re-inserting.
                    auto testIns = eIdMap.insert(make_pair(
                        edgeIds[cnt],
                        make_pair(vertIds[tmp + j],
                                  vertIds[tmp + ((j + 1) % faceVerts[i])])));
 
                    if (testIns.second == false)
                    {
                        continue;
                    }
 
                    // If the edge is reversed with respect to the face, then
                    // swap the edges so that we have the original ordering of
                    // the edge in the 3D element. This is necessary to properly
                    // determine edge orientation. Note that the logic relies on
                    // the fact that all edge forward directions are CCW
                    // orientated: we use a tensor product ordering for 2D
                    // elements so need to reverse this for edge IDs 2 and 3.
                    StdRegions::Orientation edgeOrient =
                        static_cast<StdRegions::Orientation>(edgeOrt[cnt]);
                    if (j > 1)
                    {
                        edgeOrient = edgeOrient == StdRegions::eForwards
                                         ? StdRegions::eBackwards
                                         : StdRegions::eForwards;
                    }
 
                    if (edgeOrient == StdRegions::eBackwards)
                    {
                        swap(testIns.first->second.first,
                             testIns.first->second.second);
                    }
                }
 
                vertMap[faceIds[i]]  = verts;
                edgeMap[faceIds[i]]  = edges;
                coordMap[faceIds[i]] = coord;
            }
 
            // Go through list of composites and figure out which edges are
            // parallel from original ordering in session file. This includes
            // composites which are not necessarily on this process.
 
            // Store temporary map of periodic vertices which will hold all
            // periodic vertices on the entire mesh so that doubly periodic
            // vertices/edges can be counted properly across partitions. Local
            // vertices/edges are copied into m_periodicVerts and
            // m_periodicEdges at the end of the function.
            PeriodicMap periodicVerts, periodicEdges;
 
            // Construct two maps which determine how vertices and edges of
            // faces connect given a specific face orientation. The key of the
            // map is the number of vertices in the face, used to determine
            // difference between tris and quads.
            map<int, map<StdRegions::Orientation, vector<int>>> vmap;
            map<int, map<StdRegions::Orientation, vector<int>>> emap;
 
            map<StdRegions::Orientation, vector<int>> quadVertMap;
            quadVertMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2, 3};
            quadVertMap[StdRegions::eDir1FwdDir1_Dir2BwdDir2] = {3, 2, 1, 0};
            quadVertMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {1, 0, 3, 2};
            quadVertMap[StdRegions::eDir1BwdDir1_Dir2BwdDir2] = {2, 3, 0, 1};
            quadVertMap[StdRegions::eDir1FwdDir2_Dir2FwdDir1] = {0, 3, 2, 1};
            quadVertMap[StdRegions::eDir1FwdDir2_Dir2BwdDir1] = {1, 2, 3, 0};
            quadVertMap[StdRegions::eDir1BwdDir2_Dir2FwdDir1] = {3, 0, 1, 2};
            quadVertMap[StdRegions::eDir1BwdDir2_Dir2BwdDir1] = {2, 1, 0, 3};
 
            map<StdRegions::Orientation, vector<int>> quadEdgeMap;
            quadEdgeMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2, 3};
            quadEdgeMap[StdRegions::eDir1FwdDir1_Dir2BwdDir2] = {2, 1, 0, 3};
            quadEdgeMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {0, 3, 2, 1};
            quadEdgeMap[StdRegions::eDir1BwdDir1_Dir2BwdDir2] = {2, 3, 0, 1};
            quadEdgeMap[StdRegions::eDir1FwdDir2_Dir2FwdDir1] = {3, 2, 1, 0};
            quadEdgeMap[StdRegions::eDir1FwdDir2_Dir2BwdDir1] = {1, 2, 3, 0};
            quadEdgeMap[StdRegions::eDir1BwdDir2_Dir2FwdDir1] = {3, 0, 1, 2};
            quadEdgeMap[StdRegions::eDir1BwdDir2_Dir2BwdDir1] = {1, 0, 3, 2};
 
            map<StdRegions::Orientation, vector<int>> triVertMap;
            triVertMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2};
            triVertMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {1, 0, 2};
 
            map<StdRegions::Orientation, vector<int>> triEdgeMap;
            triEdgeMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2};
            triEdgeMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {0, 2, 1};
 
            vmap[3] = triVertMap;
            vmap[4] = quadVertMap;
            emap[3] = triEdgeMap;
            emap[4] = quadEdgeMap;
 
            map<int, int> allCompPairs;
 
            // Collect composite id's of each periodic face for use if rotation
            // is required
            map<int, int> fIdToCompId;
 
            // Finally we have enough information to populate the periodic
            // vertex, edge and face maps. Begin by looping over all pairs of
            // periodic composites to determine pairs of periodic faces.
            for (auto &cIt : perComps)
            {
                SpatialDomains::CompositeSharedPtr c[2];
                const int id1    = cIt.first;
                const int id2    = cIt.second;
                std::string id1s = std::to_string(id1);
                std::string id2s = std::to_string(id2);
 
                if (compMap.count(id1) > 0)
                {
                    c[0] = compMap[id1];
                }
 
                if (compMap.count(id2) > 0)
                {
                    c[1] = compMap[id2];
                }
 
                // Loop over composite ordering to construct list of all
                // periodic faces, regardless of whether they are on this
                // process.
                map<int, int> compPairs;
 
                ASSERTL0(compOrder.count(id1) > 0,
                         "Unable to find composite " + id1s + " in order map.");
                ASSERTL0(compOrder.count(id2) > 0,
                         "Unable to find composite " + id2s + " in order map.");
                WARNINGL0(
                    compOrder[id1].size() == compOrder[id2].size(),
                    "Periodic composites " + id1s + " and " + id2s +
                        " should have the same number of elements. Have " +
                        std::to_string(compOrder[id1].size()) + " vs " +
                        std::to_string(compOrder[id2].size()));
                ASSERTL0(compOrder[id1].size() > 0, "Periodic composites " +
                                                        id1s + " and " + id2s +
                                                        " are empty!");
 
                // Look up composite ordering to determine pairs.
                for (i = 0; i < compOrder[id1].size(); ++i)
                {
                    int eId1 = compOrder[id1][i];
                    int eId2 = compOrder[id2][i];
 
                    ASSERTL0(compPairs.count(eId1) == 0, "Already paired.");
 
                    // Sanity check that the faces are mutually periodic.
                    if (compPairs.count(eId2) != 0)
                    {
                        ASSERTL0(compPairs[eId2] == eId1, "Pairing incorrect");
                    }
                    compPairs[eId1] = eId2;
 
                    // store  a map of face ids to composite ids
                    fIdToCompId[eId1] = id1;
                    fIdToCompId[eId2] = id2;
                }
 
                // Now that we have all pairs of periodic faces, loop over the
                // ones local on this process and populate face/edge/vertex
                // maps.
                for (auto &pIt : compPairs)
                {
                    int ids[2]    = {pIt.first, pIt.second};
                    bool local[2] = {locFaces.count(pIt.first) > 0,
                                     locFaces.count(pIt.second) > 0};
 
                    ASSERTL0(coordMap.count(ids[0]) > 0 &&
                                 coordMap.count(ids[1]) > 0,
                             "Unable to find face in coordinate map");
 
                    allCompPairs[pIt.first]  = pIt.second;
                    allCompPairs[pIt.second] = pIt.first;
 
                    // Loop up coordinates of the faces, check they have the
                    // same number of vertices.
                    std::vector<SpatialDomains::PointGeom *> tmpVec[2] = {
                        coordMap[ids[0]], coordMap[ids[1]]};
 
                    ASSERTL0(tmpVec[0].size() == tmpVec[1].size(),
                             "Two periodic faces have different number "
                             "of vertices!");
 
                    // o will store relative orientation of faces. Note that in
                    // some transpose cases (Dir1FwdDir2_Dir2BwdDir1 and
                    // Dir1BwdDir1_Dir2FwdDir1) it seems orientation will be
                    // different going from face1->face2 instead of face2->face1
                    // (check this).
                    StdRegions::Orientation o;
                    bool rotbnd     = false;
                    int dir         = 0;
                    NekDouble angle = 0.0;
                    NekDouble sign  = 0.0;
                    NekDouble tol   = 1e-8;
 
                    // check to see if perioid boundary is rotated
                    if (rotComp.count(fIdToCompId[pIt.first]))
                    {
                        rotbnd = true;
                        dir    = rotComp[fIdToCompId[pIt.first]].m_dir;
                        angle  = rotComp[fIdToCompId[pIt.first]].m_angle;
                        tol    = rotComp[fIdToCompId[pIt.first]].m_tol;
                    }
 
                    // Record periodic faces.
                    for (i = 0; i < 2; ++i)
                    {
                        if (!local[i])
                        {
                            continue;
                        }
 
                        // Reference to the other face.
                        int other = (i + 1) % 2;
 
                        // angle is set up for i = 0 to i = 1
                        sign = (i == 0) ? 1.0 : -1.0;
 
                        // Calculate relative face orientation.
                        if (tmpVec[0].size() == 3)
                        {
                            o = SpatialDomains::TriGeom::GetFaceOrientation(
                                {tmpVec[i][0], tmpVec[i][1], tmpVec[i][2]},
                                {tmpVec[other][0], tmpVec[other][1],
                                 tmpVec[other][2]},
                                rotbnd, dir, sign * angle, tol);
                        }
                        else
                        {
                            o = SpatialDomains::QuadGeom::GetFaceOrientation(
                                {tmpVec[i][0], tmpVec[i][1], tmpVec[i][2],
                                 tmpVec[i][3]},
                                {tmpVec[other][0], tmpVec[other][1],
                                 tmpVec[other][2], tmpVec[other][3]},
                                rotbnd, dir, sign * angle, tol);
                        }
 
                        // Record face ID, orientation and whether other face is
                        // local.
                        PeriodicEntity ent(ids[other], o, local[other]);
                        m_periodicFaces[ids[i]].push_back(ent);
                    }
 
                    int nFaceVerts = vertMap[ids[0]].size();
 
                    // Determine periodic vertices.
                    for (i = 0; i < 2; ++i)
                    {
                        int other = (i + 1) % 2;
 
                        // angle is set up for i = 0 to i = 1
                        sign = (i == 0) ? 1.0 : -1.0;
 
                        // Calculate relative face orientation.
                        if (tmpVec[0].size() == 3)
                        {
                            o = SpatialDomains::TriGeom::GetFaceOrientation(
                                {tmpVec[i][0], tmpVec[i][1], tmpVec[i][2]},
                                {tmpVec[other][0], tmpVec[other][1],
                                 tmpVec[other][2]},
                                rotbnd, dir, sign * angle, tol);
                        }
                        else
                        {
                            o = SpatialDomains::QuadGeom::GetFaceOrientation(
                                {tmpVec[i][0], tmpVec[i][1], tmpVec[i][2],
                                 tmpVec[i][3]},
                                {tmpVec[other][0], tmpVec[other][1],
                                 tmpVec[other][2], tmpVec[other][3]},
                                rotbnd, dir, sign * angle, tol);
                        }
 
                        if (nFaceVerts == 3)
                        {
                            ASSERTL0(
                                o == StdRegions::eDir1FwdDir1_Dir2FwdDir2 ||
                                    o == StdRegions::eDir1BwdDir1_Dir2FwdDir2,
                                "Unsupported face orientation for face " +
                                    std::to_string(ids[i]));
                        }
 
                        // Look up vertices for this face.
                        vector<int> per1 = vertMap[ids[i]];
                        vector<int> per2 = vertMap[ids[other]];
 
                        // tmpMap will hold the pairs of vertices which are
                        // periodic.
                        map<int, pair<int, bool>> tmpMap;
 
                        // Use vmap to determine which vertices connect given
                        // the orientation o.
                        for (j = 0; j < nFaceVerts; ++j)
                        {
                            int v = vmap[nFaceVerts][o][j];
                            tmpMap[per1[j]] =
                                make_pair(per2[v], locVerts.count(per2[v]) > 0);
                        }
 
                        // Now loop over tmpMap to associate periodic vertices.
                        for (auto &mIt : tmpMap)
                        {
                            PeriodicEntity ent2(mIt.second.first,
                                                StdRegions::eNoOrientation,
                                                mIt.second.second);
 
                            // See if this vertex has been recorded already.
                            auto perIt = periodicVerts.find(mIt.first);
 
                            if (perIt == periodicVerts.end())
                            {
                                // Vertex is new - just record this entity as
                                // usual.
                                periodicVerts[mIt.first].push_back(ent2);
                                perIt = periodicVerts.find(mIt.first);
                            }
                            else
                            {
                                // Vertex is known - loop over the vertices
                                // inside the record and potentially add vertex
                                // mIt.second to the list.
                                for (k = 0; k < perIt->second.size(); ++k)
                                {
                                    if (perIt->second[k].id == mIt.second.first)
                                    {
                                        break;
                                    }
                                }
 
                                if (k == perIt->second.size())
                                {
                                    perIt->second.push_back(ent2);
                                }
                            }
                        }
                    }
 
                    // Determine periodic edges. Logic is the same as above,
                    // and perhaps should be condensed to avoid replication.
                    for (i = 0; i < 2; ++i)
                    {
                        int other = (i + 1) % 2;
 
                        // angle is set up for i = 0 to i = 1
                        sign = (i == 0) ? 1.0 : -1.0;
 
                        if (tmpVec[0].size() == 3)
                        {
                            o = SpatialDomains::TriGeom::GetFaceOrientation(
                                {tmpVec[i][0], tmpVec[i][1], tmpVec[i][2]},
                                {tmpVec[other][0], tmpVec[other][1],
                                 tmpVec[other][2]},
                                rotbnd, dir, sign * angle, tol);
                        }
                        else
                        {
                            o = SpatialDomains::QuadGeom::GetFaceOrientation(
                                {tmpVec[i][0], tmpVec[i][1], tmpVec[i][2],
                                 tmpVec[i][3]},
                                {tmpVec[other][0], tmpVec[other][1],
                                 tmpVec[other][2], tmpVec[other][3]},
                                rotbnd, dir, sign * angle, tol);
                        }
 
                        vector<int> per1 = edgeMap[ids[i]];
                        vector<int> per2 = edgeMap[ids[other]];
 
                        map<int, pair<int, bool>> tmpMap;
 
                        for (j = 0; j < nFaceVerts; ++j)
                        {
                            int e = emap[nFaceVerts][o][j];
                            tmpMap[per1[j]] =
                                make_pair(per2[e], locEdges.count(per2[e]) > 0);
                        }
 
                        for (auto &mIt : tmpMap)
                        {
                            // Note we assume orientation of edges is forwards -
                            // this may be reversed later.
                            PeriodicEntity ent2(mIt.second.first,
                                                StdRegions::eForwards,
                                                mIt.second.second);
                            auto perIt = periodicEdges.find(mIt.first);
 
                            if (perIt == periodicEdges.end())
                            {
                                periodicEdges[mIt.first].push_back(ent2);
                                perIt = periodicEdges.find(mIt.first);
                            }
                            else
                            {
                                for (k = 0; k < perIt->second.size(); ++k)
                                {
                                    if (perIt->second[k].id == mIt.second.first)
                                    {
                                        break;
                                    }
                                }
 
                                if (k == perIt->second.size())
                                {
                                    perIt->second.push_back(ent2);
                                }
                            }
                        }
                    }
                }
            }
 
            // Make sure that the nubmer of face pairs and the
            // face Id to composite Id map match in size
            ASSERTL1(allCompPairs.size() == fIdToCompId.size(),
                     "At this point the size of allCompPairs "
                     "should have been the same as fIdToCompId");
 
            // also will need an edge id to composite id at
            // end of routine
            map<int, int> eIdToCompId;
 
            // Search for periodic vertices and edges which are not
            // in a periodic composite but lie in this process. First,
            // loop over all information we have from other
            // processors.
            for (cnt = i = 0; i < totFaces; ++i)
            {
                bool rotbnd     = false;
                int dir         = 0;
                NekDouble angle = 0.0;
                NekDouble tol   = 1e-8;
 
                int faceId = faceIds[i];
 
                ASSERTL0(allCompPairs.count(faceId) > 0,
                         "Unable to find matching periodic face. faceId = " +
                             std::to_string(faceId));
 
                int perFaceId = allCompPairs[faceId];
 
                // check to see if periodic boundary is rotated
                ASSERTL1((rotComp.size() == 0) || fIdToCompId.count(faceId) > 0,
                         "Face " + std::to_string(faceId) +
                             " not found in fIdtoCompId map");
                if (rotComp.count(fIdToCompId[faceId]))
                {
                    rotbnd = true;
                    dir    = rotComp[fIdToCompId[faceId]].m_dir;
                    angle  = rotComp[fIdToCompId[faceId]].m_angle;
                    tol    = rotComp[fIdToCompId[faceId]].m_tol;
                }
 
                for (j = 0; j < faceVerts[i]; ++j, ++cnt)
                {
                    int vId = vertIds[cnt];
 
                    auto perId = periodicVerts.find(vId);
 
                    if (perId == periodicVerts.end())
                    {
 
                        // This vertex is not included in the
                        // map. Figure out which vertex it is supposed
                        // to be periodic with. perFaceId is the face
                        // ID which is periodic with faceId. The logic
                        // is much the same as the loop above.
                        std::vector<SpatialDomains::PointGeom *> tmpVec[2] = {
                            coordMap[faceId], coordMap[perFaceId]};
 
                        int nFaceVerts = tmpVec[0].size();
                        StdRegions::Orientation o =
                            nFaceVerts == 3
                                ? SpatialDomains::TriGeom::GetFaceOrientation(
                                      {tmpVec[0][0], tmpVec[0][1],
                                       tmpVec[0][2]},
                                      {tmpVec[1][0], tmpVec[1][1],
                                       tmpVec[1][2]},
                                      rotbnd, dir, angle, tol)
                                : SpatialDomains::QuadGeom::GetFaceOrientation(
                                      {tmpVec[0][0], tmpVec[0][1], tmpVec[0][2],
                                       tmpVec[0][3]},
                                      {tmpVec[1][0], tmpVec[1][1], tmpVec[1][2],
                                       tmpVec[1][3]},
                                      rotbnd, dir, angle, tol);
 
                        // Use vmap to determine which vertex of the other face
                        // should be periodic with this one.
                        int perVertexId =
                            vertMap[perFaceId][vmap[nFaceVerts][o][j]];
 
                        PeriodicEntity ent(perVertexId,
                                           StdRegions::eNoOrientation,
                                           locVerts.count(perVertexId) > 0);
 
                        periodicVerts[vId].push_back(ent);
                    }
 
                    int eId = edgeIds[cnt];
 
                    perId = periodicEdges.find(eId);
 
                    // this map is required at very end to determine rotation of
                    // edges.
                    if (rotbnd)
                    {
                        eIdToCompId[eId] = fIdToCompId[faceId];
                    }
 
                    if (perId == periodicEdges.end())
                    {
                        // This edge is not included in the map. Figure
                        // out which edge it is supposed to be periodic
                        // with. perFaceId is the face ID which is
                        // periodic with faceId. The logic is much the
                        // same as the loop above.
                        std::vector<SpatialDomains::PointGeom *> tmpVec[2] = {
                            coordMap[faceId], coordMap[perFaceId]};
 
                        int nFaceEdges = tmpVec[0].size();
                        StdRegions::Orientation o =
                            nFaceEdges == 3
                                ? SpatialDomains::TriGeom::GetFaceOrientation(
                                      {tmpVec[0][0], tmpVec[0][1],
                                       tmpVec[0][2]},
                                      {tmpVec[1][0], tmpVec[1][1],
                                       tmpVec[1][2]},
                                      rotbnd, dir, angle, tol)
                                : SpatialDomains::QuadGeom::GetFaceOrientation(
                                      {tmpVec[0][0], tmpVec[0][1], tmpVec[0][2],
                                       tmpVec[0][3]},
                                      {tmpVec[1][0], tmpVec[1][1], tmpVec[1][2],
                                       tmpVec[1][3]},
                                      rotbnd, dir, angle, tol);
 
                        // Use emap to determine which edge of the other
                        // face should be periodic with this one.
                        int perEdgeId =
                            edgeMap[perFaceId][emap[nFaceEdges][o][j]];
 
                        PeriodicEntity ent(perEdgeId, StdRegions::eForwards,
                                           locEdges.count(perEdgeId) > 0);
 
                        periodicEdges[eId].push_back(ent);
 
                        // this map is required at very end to
                        // determine rotation of edges.
                        if (rotbnd)
                        {
                            eIdToCompId[perEdgeId] = fIdToCompId[perFaceId];
                        }
                    }
                }
            }
 
            // Finally, we must loop over the periodicVerts and periodicEdges
            // map to complete connectivity information.
            for (auto &perIt : periodicVerts)
            {
                // For each vertex that is periodic with this one...
                for (i = 0; i < perIt.second.size(); ++i)
                {
                    // Find the vertex in the periodicVerts map...
                    auto perIt2 = periodicVerts.find(perIt.second[i].id);
                    ASSERTL0(perIt2 != periodicVerts.end(),
                             "Couldn't find periodic vertex.");
 
                    // Now search through this vertex's list and make sure that
                    // we have a record of any vertices which aren't in the
                    // original list.
                    for (j = 0; j < perIt2->second.size(); ++j)
                    {
                        if (perIt2->second[j].id == perIt.first)
                        {
                            continue;
                        }
 
                        for (k = 0; k < perIt.second.size(); ++k)
                        {
                            if (perIt2->second[j].id == perIt.second[k].id)
                            {
                                break;
                            }
                        }
 
                        if (k == perIt.second.size())
                        {
                            perIt.second.push_back(perIt2->second[j]);
                        }
                    }
                }
            }
 
            for (auto &perIt : periodicEdges)
            {
                for (i = 0; i < perIt.second.size(); ++i)
                {
                    auto perIt2 = periodicEdges.find(perIt.second[i].id);
                    ASSERTL0(perIt2 != periodicEdges.end(),
                             "Couldn't find periodic edge.");
 
                    for (j = 0; j < perIt2->second.size(); ++j)
                    {
                        if (perIt2->second[j].id == perIt.first)
                        {
                            continue;
                        }
 
                        for (k = 0; k < perIt.second.size(); ++k)
                        {
                            if (perIt2->second[j].id == perIt.second[k].id)
                            {
                                break;
                            }
                        }
 
                        if (k == perIt.second.size())
                        {
                            perIt.second.push_back(perIt2->second[j]);
                        }
                    }
                }
            }
 
            // Loop over periodic edges to determine relative edge orientations.
            for (auto &perIt : periodicEdges)
            {
                bool rotbnd     = false;
                int dir         = 0;
                NekDouble angle = 0.0;
                NekDouble tol   = 1e-8;
 
                // Find edge coordinates
                auto eIt                       = eIdMap.find(perIt.first);
                SpatialDomains::PointGeom v[2] = {*vCoMap[eIt->second.first],
                                                  *vCoMap[eIt->second.second]};
 
                // check to see if perioid boundary is rotated
                if (rotComp.count(eIdToCompId[perIt.first]))
                {
                    rotbnd = true;
                    dir    = rotComp[eIdToCompId[perIt.first]].m_dir;
                    angle  = rotComp[eIdToCompId[perIt.first]].m_angle;
                    tol    = rotComp[eIdToCompId[perIt.first]].m_tol;
                }
 
                // Loop over each edge, and construct a vector that takes us
                // from one vertex to another. Use this to figure out which
                // vertex maps to which.
                for (i = 0; i < perIt.second.size(); ++i)
                {
                    eIt = eIdMap.find(perIt.second[i].id);
 
                    SpatialDomains::PointGeom w[2] = {
                        *vCoMap[eIt->second.first],
                        *vCoMap[eIt->second.second]};
 
                    int vMap[2] = {-1, -1};
                    if (rotbnd)
                    {
 
                        SpatialDomains::PointGeom r;
 
                        r.Rotate(v[0], dir, angle);
 
                        if (r.dist(w[0]) < tol)
                        {
                            vMap[0] = 0;
                        }
                        else
                        {
                            r.Rotate(v[1], dir, angle);
                            if (r.dist(w[0]) < tol)
                            {
                                vMap[0] = 1;
                            }
                            else
                            {
                                NEKERROR(ErrorUtil::efatal,
                                         "Unable to align rotationally "
                                         "periodic edge vertex");
                            }
                        }
                    }
                    else // translation test
                    {
                        NekDouble cx =
                            0.5 * (w[0](0) - v[0](0) + w[1](0) - v[1](0));
                        NekDouble cy =
                            0.5 * (w[0](1) - v[0](1) + w[1](1) - v[1](1));
                        NekDouble cz =
                            0.5 * (w[0](2) - v[0](2) + w[1](2) - v[1](2));
 
                        for (j = 0; j < 2; ++j)
                        {
                            NekDouble x = v[j](0);
                            NekDouble y = v[j](1);
                            NekDouble z = v[j](2);
                            for (k = 0; k < 2; ++k)
                            {
                                NekDouble x1 = w[k](0) - cx;
                                NekDouble y1 = w[k](1) - cy;
                                NekDouble z1 = w[k](2) - cz;
 
                                if (sqrt((x1 - x) * (x1 - x) +
                                         (y1 - y) * (y1 - y) +
                                         (z1 - z) * (z1 - z)) < 1e-8)
                                {
                                    vMap[k] = j;
                                    break;
                                }
                            }
                        }
 
                        // Sanity check the map.
                        ASSERTL0(vMap[0] >= 0 && vMap[1] >= 0,
                                 "Unable to align periodic edge vertex.");
                        ASSERTL0((vMap[0] == 0 || vMap[0] == 1) &&
                                     (vMap[1] == 0 || vMap[1] == 1) &&
                                     (vMap[0] != vMap[1]),
                                 "Unable to align periodic edge vertex.");
                    }
 
                    // If 0 -> 0 then edges are aligned already; otherwise
                    // reverse the orientation.
                    if (vMap[0] != 0)
                    {
                        perIt.second[i].orient = StdRegions::eBackwards;
                    }
                }
            }
 
            // Do one final loop over periodic vertices/edges to remove
            // non-local vertices/edges from map.
            for (auto &perIt : periodicVerts)
            {
                if (locVerts.count(perIt.first) > 0)
                {
                    m_periodicVerts.insert(perIt);
                }
            }
 
            for (auto &perIt : periodicEdges)
            {
                if (locEdges.count(perIt.first) > 0)
                {
                    m_periodicEdges.insert(perIt);
                }
            }
        }
        break;
        default:
            ASSERTL1(false, "Not setup for this expansion");
            break;
    }
}
 
/**
 *
 */
GlobalLinSysSharedPtr DisContField::GetGlobalBndLinSys(
    const GlobalLinSysKey &mkey)
{
    ASSERTL0(mkey.GetMatrixType() == StdRegions::eHybridDGHelmBndLam,
             "Routine currently only tested for HybridDGHelmholtz");
 
    ASSERTL1(mkey.GetGlobalSysSolnType() != eDirectFullMatrix,
             "Full matrix global systems are not supported for HDG "
             "expansions");
 
    ASSERTL1(mkey.GetGlobalSysSolnType() == m_traceMap->GetGlobalSysSolnType(),
             "The local to global map is not set up for the requested "
             "solution type");
 
    GlobalLinSysSharedPtr glo_matrix;
    auto matrixIter = m_globalBndMat->find(mkey);
 
    if (matrixIter == m_globalBndMat->end())
    {
        glo_matrix              = GenGlobalBndLinSys(mkey, m_traceMap);
        (*m_globalBndMat)[mkey] = glo_matrix;
    }
    else
    {
        glo_matrix = matrixIter->second;
    }
 
    return glo_matrix;
}
 
/**
 * For each boundary region, checks that the types and number of
 * boundary expansions in that region match.
 * @param   In          Field to compare with.
 * @return True if boundary conditions match.
 */
bool DisContField::SameTypeOfBoundaryConditions(const DisContField &In)
{
    int i;
    bool returnval = true;
 
    for (i = 0; i < m_bndConditions.size(); ++i)
    {
        // check to see if boundary condition type is the same
        // and there are the same number of boundary
        // conditions in the boundary definition.
        if ((m_bndConditions[i]->GetBoundaryConditionType() !=
             In.m_bndConditions[i]->GetBoundaryConditionType()) ||
            (m_bndCondExpansions[i]->GetExpSize() !=
             In.m_bndCondExpansions[i]->GetExpSize()))
        {
            returnval = false;
            break;
        }
    }
 
    // Compare with all other processes. Return true only if all
    // processes report having the same boundary conditions.
    int vSame = (returnval ? 1 : 0);
    m_comm->GetRowComm()->AllReduce(vSame, LibUtilities::ReduceMin);
 
    return (vSame == 1);
}
 
vector<bool> &DisContField::GetNegatedFluxNormal(void)
{
    if (m_negatedFluxNormal.size() == 0)
    {
        Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
            &elmtToTrace = m_traceMap->GetElmtToTrace();
 
        m_negatedFluxNormal.resize(2 * GetExpSize());
 
        for (int i = 0; i < GetExpSize(); ++i)
        {
 
            for (int v = 0; v < 2; ++v)
            {
                LocalRegions::Expansion0DSharedPtr vertExp =
                    elmtToTrace[i][v]->as<LocalRegions::Expansion0D>();
 
                if (vertExp->GetLeftAdjacentElementExp()
                        ->GetGeom()
                        ->GetGlobalID() !=
                    (*m_exp)[i]->GetGeom()->GetGlobalID())
                {
                    m_negatedFluxNormal[2 * i + v] = true;
                }
                else
                {
                    m_negatedFluxNormal[2 * i + v] = false;
                }
            }
        }
    }
 
    return m_negatedFluxNormal;
}
 
ExpListSharedPtr &DisContField::v_GetTrace()
{
    if (m_trace == NullExpListSharedPtr)
    {
        SetUpDG();
    }
 
    return m_trace;
}
 
AssemblyMapDGSharedPtr &DisContField::v_GetTraceMap(void)
{
    return m_traceMap;
}
 
InterfaceMapDGSharedPtr &DisContField::v_GetInterfaceMap(void)
{
    return m_interfaceMap;
}
 
const LocTraceToTraceMapSharedPtr &DisContField::v_GetLocTraceToTraceMap(
    void) const
{
    return m_locTraceToTraceMap;
}
 
// Return the internal vector which identifieds if trace
// is left adjacent definiing which trace the normal
// points otwards from
std::vector<bool> &DisContField::v_GetLeftAdjacentTraces(void)
{
    return m_leftAdjacentTraces;
}
 
const Array<OneD, const MultiRegions::ExpListSharedPtr> &DisContField::
    v_GetBndCondExpansions()
{
    return m_bndCondExpansions;
}
 
const Array<OneD, const SpatialDomains::BoundaryConditionShPtr> &DisContField::
    v_GetBndConditions()
{
    return m_bndConditions;
}
 
MultiRegions::ExpListSharedPtr &DisContField::v_UpdateBndCondExpansion(int i)
{
    return m_bndCondExpansions[i];
}
 
Array<OneD, SpatialDomains::BoundaryConditionShPtr> &DisContField::
    v_UpdateBndConditions()
{
    return m_bndConditions;
}
 
void DisContField::v_PeriodicBwdCopy(const Array<OneD, const NekDouble> &Fwd,
                                     Array<OneD, NekDouble> &Bwd)
{
    for (int n = 0; n < m_periodicFwdCopy.size(); ++n)
    {
        Bwd[m_periodicBwdCopy[n]] = Fwd[m_periodicFwdCopy[n]];
    }
}
 
void DisContField::v_PeriodicBwdRot(Array<OneD, Array<OneD, NekDouble>> &Bwd)
{
    int cnt = 0;
    for (int n = 0; n < m_bndConditions.size(); ++n)
    {
        // check to see if boundary is rotationally aligned
        if (boost::icontains(m_bndConditions[n]->GetUserDefined(), "Rotated"))
        {
            vector<string> tmpstr;
 
            boost::split(tmpstr, m_bndConditions[n]->GetUserDefined(),
                         boost::is_any_of(":"));
 
            ASSERTL1(tmpstr.size() > 2,
                     "Expected Rotated user defined string to "
                     "contain direction and rotation angle "
                     "and optionally a tolerance, "
                     "i.e. Rotated:dir:PI/2:1e-6");
 
            ASSERTL1((tmpstr[1] == "x") || (tmpstr[1] == "y") ||
                         (tmpstr[1] == "z"),
                     "Rotated Dir is "
                     "not specified as x,y or z");
 
            RotPeriodicInfo RotInfo;
            RotInfo.m_dir = (tmpstr[1] == "x") ? 0 : (tmpstr[1] == "y") ? 1 : 2;
 
            LibUtilities::Interpreter strEval;
            int ExprId      = strEval.DefineFunction("", tmpstr[2]);
            RotInfo.m_angle = strEval.Evaluate(ExprId);
 
            auto ne = m_bndCondExpansions[n]->GetExpSize();
 
            // Loop over each element in the boundary condition and rotate
            // velocity
            for (int e = 0; e < ne; ++e)
            {
                auto id2 = m_trace->GetPhys_Offset(
                    m_traceMap->GetPerBndCondIDToGlobalTraceID(e + cnt));
                int npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                Rotate(Bwd, RotInfo.m_dir, -RotInfo.m_angle, id2, npts);
            }
            cnt += ne;
        }
    }
}
 
void DisContField::v_PeriodicDeriveBwdRot(TensorOfArray3D<NekDouble> &Bwd)
{
    int cnt = 0;
    for (int n = 0; n < m_bndConditions.size(); ++n)
    {
        // check to see if boundary is rotationally aligned
        if (boost::icontains(m_bndConditions[n]->GetUserDefined(), "Rotated"))
        {
            vector<string> tmpstr;
 
            boost::split(tmpstr, m_bndConditions[n]->GetUserDefined(),
                         boost::is_any_of(":"));
 
            ASSERTL1(tmpstr.size() > 2,
                     "Expected Rotated user defined string to "
                     "contain direction and rotation angle "
                     "and optionally a tolerance, "
                     "i.e. Rotated:dir:PI/2:1e-6");
 
            ASSERTL1((tmpstr[1] == "x") || (tmpstr[1] == "y") ||
                         (tmpstr[1] == "z"),
                     "Rotated Dir is "
                     "not specified as x,y or z");
 
            RotPeriodicInfo RotInfo;
            RotInfo.m_dir = (tmpstr[1] == "x") ? 0 : (tmpstr[1] == "y") ? 1 : 2;
 
            LibUtilities::Interpreter strEval;
            int ExprId      = strEval.DefineFunction("", tmpstr[2]);
            RotInfo.m_angle = strEval.Evaluate(ExprId);
 
            auto ne = m_bndCondExpansions[n]->GetExpSize();
            // Loop over each element in the boundary condition and rotate
            // velocity
            for (int e = 0; e < ne; ++e)
            {
                auto id2 = m_trace->GetPhys_Offset(
                    m_traceMap->GetPerBndCondIDToGlobalTraceID(e + cnt));
                int npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                DeriveRotate(Bwd, RotInfo.m_dir, -RotInfo.m_angle, id2, npts);
            }
            cnt += ne;
        }
    }
}
 
void DisContField::v_RotLocalBwdTrace(Array<OneD, Array<OneD, NekDouble>> &Bwd)
{
    int nDim = GetCoordim(0);
    for (auto &interfaceTrace : m_interfaceMap->GetLocalInterface())
    {
 
        interfaceTrace->RotLocalBwdTrace(Bwd, nDim);
    }
}
 
void DisContField::v_RotLocalBwdDeriveTrace(TensorOfArray3D<NekDouble> &Bwd)
{
    int nDim = GetCoordim(0);
    for (auto &interfaceTrace : m_interfaceMap->GetLocalInterface())
    {
 
        interfaceTrace->RotLocalBwdDeriveTrace(Bwd, nDim);
    }
}
 
void DisContField::v_GetFwdBwdTracePhys(Array<OneD, NekDouble> &Fwd,
                                        Array<OneD, NekDouble> &Bwd)
{
    v_GetFwdBwdTracePhys(m_phys, Fwd, Bwd);
}
 
/**
 * @brief Obtain a copy of the periodic edges and vertices for this
 * field.
 */
void DisContField::v_GetPeriodicEntities(PeriodicMap &periodicVerts,
                                         PeriodicMap &periodicEdges,
                                         PeriodicMap &periodicFaces)
{
    periodicVerts = m_periodicVerts;
    periodicEdges = m_periodicEdges;
    periodicFaces = m_periodicFaces;
}
/**
 * \brief This method extracts the "forward" and "backward" trace
 * data from the array \a field and puts the data into output
 * vectors \a Fwd and \a Bwd.
 *
 * We first define the convention which defines "forwards" and
 * "backwards". First an association is made between the vertex/edge/face of
 * each element and its corresponding vertex/edge/face in the trace space
 * using the mapping #m_traceMap. The element can either be
 * left-adjacent or right-adjacent to this trace face (see
 * Expansion2D::GetLeftAdjacentElementExp). Boundary faces are
 * always left-adjacent since left-adjacency is populated first.
 *
 * If the element is left-adjacent we extract the trace data
 * from \a field into the forward trace space \a Fwd; otherwise,
 * we place it in the backwards trace space \a Bwd. In this way,
 * we form a unique set of trace normals since these are always
 * extracted from left-adjacent elements.
 *
 * \param field is a NekDouble array which contains the fielddata
 * from which we wish to extract the backward and forward
 * orientated trace/face arrays.
 *
 * \return Updates a NekDouble array \a Fwd and \a Bwd
 */
void DisContField::GetFwdBwdTracePhys(
    const Array<OneD, const NekDouble> &field, Array<OneD, NekDouble> &Fwd,
    Array<OneD, NekDouble> &Bwd,
    const Array<OneD, const SpatialDomains::BoundaryConditionShPtr> &BndCond,
    const Array<OneD, const ExpListSharedPtr> &BndCondExp)
{
    ASSERTL1((*m_exp)[0]->GetBasis(0)->GetBasisType() !=
                 LibUtilities::eGauss_Lagrange,
             "Routine not set up for Gauss Lagrange points distribution");
 
    // blocked routine
    Array<OneD, NekDouble> tracevals(m_locTraceToTraceMap->GetNLocTracePts());
 
    m_locTraceToTraceMap->LocTracesFromField(field, tracevals);
    m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, Fwd);
 
    Array<OneD, NekDouble> invals =
        tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
    m_locTraceToTraceMap->InterpLocTracesToTrace(1, invals, Bwd);
 
    DisContField::v_PeriodicBwdCopy(Fwd, Bwd);
 
    FillBwdWithBoundCond(Fwd, Bwd, BndCond, BndCondExp, false);
 
    // Do parallel exchange for forwards/backwards spaces.
    m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
 
    // Do exchange of interface traces (local and parallel)
    // We may have to split this out into separate local and
    // parallel for IP method???
    m_interfaceMap->ExchangeTrace(Fwd, Bwd);
}
 
void DisContField::v_GetFwdBwdTracePhys(
    const Array<OneD, const NekDouble> &field, Array<OneD, NekDouble> &Fwd,
    Array<OneD, NekDouble> &Bwd,
    bool FillBnd, // these should be template params
    bool PutFwdInBwdOnBCs, bool DoExchange)
{
    // Is this zeroing necessary?
    // Zero forward/backward vectors.
    Vmath::Zero(Fwd.size(), Fwd, 1);
    Vmath::Zero(Bwd.size(), Bwd, 1);
 
    // Basis definition on each element
    LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
    if ((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
    {
        // blocked routine
        Array<OneD, NekDouble> tracevals(
            m_locTraceToTraceMap->GetNLocTracePts());
 
        m_locTraceToTraceMap->LocTracesFromField(field, tracevals);
        m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, Fwd);
 
        Array<OneD, NekDouble> invals =
            tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
        m_locTraceToTraceMap->InterpLocTracesToTrace(1, invals, Bwd);
    }
    else
    {
        // Loop over elements and collect forward expansion
        auto nexp = GetExpSize();
        Array<OneD, NekDouble> e_tmp;
        LocalRegions::ExpansionSharedPtr exp;
 
        Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
            &elmtToTrace = m_traceMap->GetElmtToTrace();
 
        for (int n = 0, cnt = 0; n < nexp; ++n)
        {
            exp              = (*m_exp)[n];
            auto phys_offset = GetPhys_Offset(n);
 
            for (int e = 0; e < exp->GetNtraces(); ++e, ++cnt)
            {
                auto offset =
                    m_trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
 
                e_tmp =
                    (m_leftAdjacentTraces[cnt]) ? Fwd + offset : Bwd + offset;
 
                exp->GetTracePhysVals(e, elmtToTrace[n][e], field + phys_offset,
                                      e_tmp);
            }
        }
    }
 
    DisContField::v_PeriodicBwdCopy(Fwd, Bwd);
 
    if (FillBnd)
    {
        v_FillBwdWithBoundCond(Fwd, Bwd, PutFwdInBwdOnBCs);
    }
 
    if (DoExchange)
    {
        // Do parallel exchange for forwards/backwards spaces.
        m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
 
        // Do exchange of interface traces (local and parallel)
        // We may have to split this out into separate local and
        // parallel for IP method???
        m_interfaceMap->ExchangeTrace(Fwd, Bwd);
    }
}
 
void DisContField::v_FillBwdWithBoundCond(const Array<OneD, NekDouble> &Fwd,
                                          Array<OneD, NekDouble> &Bwd,
                                          bool PutFwdInBwdOnBCs)
{
    FillBwdWithBoundCond(Fwd, Bwd, m_bndConditions, m_bndCondExpansions,
                         PutFwdInBwdOnBCs);
}
 
/** function allowing different BCs to be passed to routine
 */
void DisContField::FillBwdWithBoundCond(
    const Array<OneD, NekDouble> &Fwd, Array<OneD, NekDouble> &Bwd,
    const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>
        &bndConditions,
    const Array<OneD, const ExpListSharedPtr> &bndCondExpansions,
    bool PutFwdInBwdOnBCs)
{
 
    // Fill boundary conditions into missing elements
    if (PutFwdInBwdOnBCs) // just set Bwd value to be Fwd value on BCs
    {
        // Fill boundary conditions into missing elements
        for (int n = 0, cnt = 0; n < bndCondExpansions.size(); ++n)
        {
            if (bndConditions[n]->GetBoundaryConditionType() ==
                SpatialDomains::eDirichlet)
            {
                auto ne = bndCondExpansions[n]->GetExpSize();
                for (int e = 0; e < ne; ++e)
                {
                    auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                    auto id2  = m_trace->GetPhys_Offset(
                        m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
                    Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
                }
 
                cnt += ne;
            }
            else if (bndConditions[n]->GetBoundaryConditionType() ==
                         SpatialDomains::eNeumann ||
                     bndConditions[n]->GetBoundaryConditionType() ==
                         SpatialDomains::eRobin)
            {
                auto ne = bndCondExpansions[n]->GetExpSize();
                for (int e = 0; e < ne; ++e)
                {
                    auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                    auto id2  = m_trace->GetPhys_Offset(
                        m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
                    Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
                }
                cnt += ne;
            }
            else if (bndConditions[n]->GetBoundaryConditionType() !=
                     SpatialDomains::ePeriodic)
            {
                ASSERTL0(false, "Method not set up for this "
                                "boundary condition.");
            }
        }
    }
    else
    {
        for (int n = 0, cnt = 0; n < bndCondExpansions.size(); ++n)
        {
            if (bndConditions[n]->GetBoundaryConditionType() ==
                SpatialDomains::eDirichlet)
            {
                auto ne = bndCondExpansions[n]->GetExpSize();
                for (int e = 0; e < ne; ++e)
                {
                    auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                    auto id1  = bndCondExpansions[n]->GetPhys_Offset(e);
                    auto id2  = m_trace->GetPhys_Offset(
                        m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
 
                    Vmath::Vcopy(npts, &(bndCondExpansions[n]->GetPhys())[id1],
                                 1, &Bwd[id2], 1);
                }
                cnt += ne;
            }
            else if (bndConditions[n]->GetBoundaryConditionType() ==
                         SpatialDomains::eNeumann ||
                     bndConditions[n]->GetBoundaryConditionType() ==
                         SpatialDomains::eRobin)
            {
                auto ne = m_bndCondExpansions[n]->GetExpSize();
                for (int e = 0; e < ne; ++e)
                {
                    auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                    auto id1  = bndCondExpansions[n]->GetPhys_Offset(e);
 
                    ASSERTL0((bndCondExpansions[n]->GetPhys())[id1] == 0.0,
                             "Method not set up for non-zero "
                             "Neumann boundary condition");
                    auto id2 = m_trace->GetPhys_Offset(
                        m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
 
                    Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
                }
 
                cnt += ne;
            }
            else if (bndConditions[n]->GetBoundaryConditionType() ==
                     SpatialDomains::eNotDefined)
            {
                // Do nothing
            }
            else if (bndConditions[n]->GetBoundaryConditionType() !=
                     SpatialDomains::ePeriodic)
            {
                NEKERROR(ErrorUtil::efatal,
                         "Method not set up for this boundary "
                         "condition.");
            }
        }
    }
}
 
const Array<OneD, const NekDouble> &DisContField::v_GetBndCondBwdWeight()
{
    return m_bndCondBndWeight;
}
 
void DisContField::v_SetBndCondBwdWeight(const int index, const NekDouble value)
{
    m_bndCondBndWeight[index] = value;
}
 
void DisContField::v_AddTraceQuadPhysToField(
    const Array<OneD, const NekDouble> &Fwd,
    const Array<OneD, const NekDouble> &Bwd, Array<OneD, NekDouble> &field)
{
    // Basis definition on each element
    LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
    if (basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
    {
        Array<OneD, NekDouble> tracevals(
            m_locTraceToTraceMap->GetNLocTracePts(), 0.0);
 
        Array<OneD, NekDouble> invals =
            tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
 
        m_locTraceToTraceMap->InterpTraceToLocTrace(1, Bwd, invals);
 
        m_locTraceToTraceMap->InterpTraceToLocTrace(0, Fwd, tracevals);
 
        m_locTraceToTraceMap->AddLocTracesToField(tracevals, field);
    }
    else
    {
        NEKERROR(ErrorUtil::efatal,
                 "v_AddTraceQuadPhysToField not coded for eGauss_Lagrange");
    }
}
 
void DisContField::v_ExtractTracePhys(Array<OneD, NekDouble> &outarray)
{
    ASSERTL1(m_physState == true, "local physical space is not true ");
    v_ExtractTracePhys(m_phys, outarray);
}
/**
 * @brief This method extracts the trace (verts in 1D) from
 * the field @a inarray and puts the values in @a outarray.
 *
 * It assumes the field is C0 continuous so that it can
 * overwrite the edge data when visited by the two adjacent
 * elements.
 *
 * @param inarray   An array containing the 1D data from which we wish
 *                  to extract the edge data.
 * @param outarray  The resulting edge information.
 *
 * @param gridVelocity Avoid performing parallel exchanges of the grid velocity
 *
 * This will not work for non-boundary expansions
 */
void DisContField::v_ExtractTracePhys(
    const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray, bool gridVelocity)
{
    LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
    if ((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
    {
        Vmath::Zero(outarray.size(), outarray, 1);
        Array<OneD, NekDouble> tracevals(
            m_locTraceToTraceMap->GetNFwdLocTracePts());
        m_locTraceToTraceMap->FwdLocTracesFromField(inarray, tracevals);
        m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, outarray);
        if (!gridVelocity)
        {
            m_traceMap->GetAssemblyCommDG()->PerformExchange(outarray,
                                                             outarray);
        }
    }
    else
    {
        // Loop over elemente and collect forward expansion
        int nexp = GetExpSize();
        int n, p, offset, phys_offset;
        Array<OneD, NekDouble> t_tmp;
        Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
            &elmtToTrace = m_traceMap->GetElmtToTrace();
        ASSERTL1(outarray.size() >= m_trace->GetNpoints(),
                 "input array is of insufficient length");
        for (n = 0; n < nexp; ++n)
        {
            phys_offset = GetPhys_Offset(n);
            for (p = 0; p < (*m_exp)[n]->GetNtraces(); ++p)
            {
                offset =
                    m_trace->GetPhys_Offset(elmtToTrace[n][p]->GetElmtId());
                (*m_exp)[n]->GetTracePhysVals(p, elmtToTrace[n][p],
                                              inarray + phys_offset,
                                              t_tmp = outarray + offset);
            }
        }
    }
}
/**
 * @brief Add trace contributions into elemental coefficient spaces.
 *
 * Given some quantity \f$ \vec{Fn} \f$, which conatins this
 * routine calculates the integral
 *
 * \f[
 * \int_{\Omega^e} \vec{Fn}, \mathrm{d}S
 * \f]
 *
 * and adds this to the coefficient space provided by outarray.
 *
 * @see Expansion3D::AddFaceNormBoundaryInt
 *
 * @param Fn        The trace quantities.
 * @param outarray  Resulting 3D coefficient space.
 *
 */
void DisContField::v_AddTraceIntegral(const Array<OneD, const NekDouble> &Fn,
                                      Array<OneD, NekDouble> &outarray)
{
    if (m_expType == e1D)
    {
        int n, offset, t_offset;
        Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
            &elmtToTrace               = m_traceMap->GetElmtToTrace();
        vector<bool> negatedFluxNormal = GetNegatedFluxNormal();
        for (n = 0; n < GetExpSize(); ++n)
        {
            // Number of coefficients on each element
            int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
            offset        = GetCoeff_Offset(n);
            // Implementation for every points except Gauss points
            if ((*m_exp)[n]->GetBasis(0)->GetBasisType() !=
                LibUtilities::eGauss_Lagrange)
            {
                t_offset =
                    GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
                if (negatedFluxNormal[2 * n])
                {
                    outarray[offset] -= Fn[t_offset];
                }
                else
                {
                    outarray[offset] += Fn[t_offset];
                }
                t_offset =
                    GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
                if (negatedFluxNormal[2 * n + 1])
                {
                    outarray[offset + (*m_exp)[n]->GetVertexMap(1)] -=
                        Fn[t_offset];
                }
                else
                {
                    outarray[offset + (*m_exp)[n]->GetVertexMap(1)] +=
                        Fn[t_offset];
                }
            }
            else
            {
                int j;
                static DNekMatSharedPtr m_Ixm, m_Ixp;
                static int sav_ncoeffs = 0;
                if (!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
                {
                    LibUtilities::BasisSharedPtr BASE;
                    const LibUtilities::PointsKey BS_p(
                        e_ncoeffs, LibUtilities::eGaussGaussLegendre);
                    const LibUtilities::BasisKey BS_k(
                        LibUtilities::eGauss_Lagrange, e_ncoeffs, BS_p);
                    BASE = LibUtilities::BasisManager()[BS_k];
                    Array<OneD, NekDouble> coords(1, 0.0);
                    coords[0]   = -1.0;
                    m_Ixm       = BASE->GetI(coords);
                    coords[0]   = 1.0;
                    m_Ixp       = BASE->GetI(coords);
                    sav_ncoeffs = e_ncoeffs;
                }
                t_offset =
                    GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
                if (negatedFluxNormal[2 * n])
                {
                    for (j = 0; j < e_ncoeffs; j++)
                    {
                        outarray[offset + j] -=
                            (m_Ixm->GetPtr())[j] * Fn[t_offset];
                    }
                }
                else
                {
                    for (j = 0; j < e_ncoeffs; j++)
                    {
                        outarray[offset + j] +=
                            (m_Ixm->GetPtr())[j] * Fn[t_offset];
                    }
                }
                t_offset =
                    GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
                if (negatedFluxNormal[2 * n + 1])
                {
                    for (j = 0; j < e_ncoeffs; j++)
                    {
                        outarray[offset + j] -=
                            (m_Ixp->GetPtr())[j] * Fn[t_offset];
                    }
                }
                else
                {
                    for (j = 0; j < e_ncoeffs; j++)
                    {
                        outarray[offset + j] +=
                            (m_Ixp->GetPtr())[j] * Fn[t_offset];
                    }
                }
            }
        }
    }
    else // other dimensions
    {
        // Basis definition on each element
        LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
        if ((m_expType != e1D) &&
            (basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
        {
            Array<OneD, NekDouble> Fcoeffs(m_trace->GetNcoeffs());
            m_trace->IProductWRTBase(Fn, Fcoeffs);
            m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(Fcoeffs,
                                                              outarray);
        }
        else
        {
            int e, n, offset, t_offset;
            Array<OneD, NekDouble> e_outarray;
            Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
                &elmtToTrace = m_traceMap->GetElmtToTrace();
            if (m_expType == e2D)
            {
                for (n = 0; n < GetExpSize(); ++n)
                {
                    offset = GetCoeff_Offset(n);
                    for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
                    {
                        t_offset = GetTrace()->GetPhys_Offset(
                            elmtToTrace[n][e]->GetElmtId());
                        (*m_exp)[n]->AddEdgeNormBoundaryInt(
                            e, elmtToTrace[n][e], Fn + t_offset,
                            e_outarray = outarray + offset);
                    }
                }
            }
            else if (m_expType == e3D)
            {
                for (n = 0; n < GetExpSize(); ++n)
                {
                    offset = GetCoeff_Offset(n);
                    for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
                    {
                        t_offset = GetTrace()->GetPhys_Offset(
                            elmtToTrace[n][e]->GetElmtId());
                        (*m_exp)[n]->AddFaceNormBoundaryInt(
                            e, elmtToTrace[n][e], Fn + t_offset,
                            e_outarray = outarray + offset);
                    }
                }
            }
        }
    }
}
/**
 * @brief Add trace contributions into elemental coefficient spaces.
 *
 * Given some quantity \f$ \vec{q} \f$, calculate the elemental integral
 *
 * \f[
 * \int_{\Omega^e} \vec{q}, \mathrm{d}S
 * \f]
 *
 * and adds this to the coefficient space provided by
 * outarray. The value of q is determined from the routine
 * IsLeftAdjacentTrace() which if true we use Fwd else we use
 * Bwd
 *
 * @param Fwd       The trace quantities associated with left (fwd)
 *                  adjancent elmt.
 * @param Bwd       The trace quantities associated with right (bwd)
 *                  adjacent elet.
 * @param outarray  Resulting  coefficient space.
 */
void DisContField::v_AddFwdBwdTraceIntegral(
    const Array<OneD, const NekDouble> &Fwd,
    const Array<OneD, const NekDouble> &Bwd, Array<OneD, NekDouble> &outarray)
{
    ASSERTL0(m_expType != e1D, "This method is not setup or "
                               "tested for 1D expansion");
    Array<OneD, NekDouble> Coeffs(m_trace->GetNcoeffs());
    m_trace->IProductWRTBase(Fwd, Coeffs);
    m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0, Coeffs, outarray);
    m_trace->IProductWRTBase(Bwd, Coeffs);
    m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1, Coeffs, outarray);
}
 
GlobalLinSysKey DisContField::v_HelmSolve(
    const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray, const StdRegions::ConstFactorMap &factors,
    const StdRegions::VarCoeffMap &varcoeff,
    [[maybe_unused]] const StdRegions::VarFactorsMap &varfactors,
    [[maybe_unused]] const Array<OneD, const NekDouble> &dirForcing,
    const bool PhysSpaceForcing)
{
    int i, n, cnt, nbndry;
    int nexp = GetExpSize();
    Array<OneD, NekDouble> f(m_ncoeffs);
    DNekVec F(m_ncoeffs, f, eWrapper);
    Array<OneD, NekDouble> e_f, e_l;
    //----------------------------------
    // Setup RHS Inner product if required
    //----------------------------------
    if (PhysSpaceForcing)
    {
        IProductWRTBase(inarray, f);
        Vmath::Neg(m_ncoeffs, f, 1);
    }
    else
    {
        Vmath::Smul(m_ncoeffs, -1.0, inarray, 1, f, 1);
    }
    //----------------------------------
    // Solve continuous Boundary System
    //----------------------------------
    int GloBndDofs = m_traceMap->GetNumGlobalBndCoeffs();
    int NumDirBCs  = m_traceMap->GetNumLocalDirBndCoeffs();
    int e_ncoeffs;
    GlobalMatrixKey HDGLamToUKey(StdRegions::eHybridDGLamToU,
                                 NullAssemblyMapSharedPtr, factors, varcoeff);
    const DNekScalBlkMatSharedPtr &HDGLamToU = GetBlockMatrix(HDGLamToUKey);
    // Retrieve number of local trace space coefficients N_{\lambda},
    // and set up local elemental trace solution \lambda^e.
    int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
    Array<OneD, NekDouble> bndrhs(LocBndCoeffs, 0.0);
    Array<OneD, NekDouble> loclambda(LocBndCoeffs, 0.0);
    DNekVec LocLambda(LocBndCoeffs, loclambda, eWrapper);
    //----------------------------------
    // Evaluate Trace Forcing
    // Kirby et al, 2010, P23, Step 5.
    //----------------------------------
    // Determing <u_lam,f> terms using HDGLamToU matrix
    for (cnt = n = 0; n < nexp; ++n)
    {
        nbndry    = (*m_exp)[n]->NumDGBndryCoeffs();
        e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
        e_f       = f + m_coeff_offset[n];
        e_l       = bndrhs + cnt;
        // use outarray as tmp space
        DNekVec Floc(nbndry, e_l, eWrapper);
        DNekVec ElmtFce(e_ncoeffs, e_f, eWrapper);
        Floc = Transpose(*(HDGLamToU->GetBlock(n, n))) * ElmtFce;
        cnt += nbndry;
    }
    Array<OneD, const int> bndCondMap =
        m_traceMap->GetBndCondCoeffsToLocalTraceMap();
    Array<OneD, const NekDouble> Sign = m_traceMap->GetLocalToGlobalBndSign();
    // Copy Dirichlet boundary conditions and weak forcing
    // into trace space
    int locid;
    cnt = 0;
    for (i = 0; i < m_bndCondExpansions.size(); ++i)
    {
        const Array<OneD, const NekDouble> bndcoeffs =
            m_bndCondExpansions[i]->GetCoeffs();
 
        if (m_bndConditions[i]->GetBoundaryConditionType() ==
            SpatialDomains::eDirichlet)
        {
            for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
            {
                locid            = bndCondMap[cnt + j];
                loclambda[locid] = Sign[locid] * bndcoeffs[j];
            }
        }
        else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                     SpatialDomains::eNeumann ||
                 m_bndConditions[i]->GetBoundaryConditionType() ==
                     SpatialDomains::eRobin)
        {
            // Add weak boundary condition to trace forcing
            for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
            {
                locid = bndCondMap[cnt + j];
                bndrhs[locid] += Sign[locid] * bndcoeffs[j];
            }
        }
        else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                 SpatialDomains::ePeriodic)
        {
            // skip increment of cnt if ePeriodic
            // because bndCondMap does not include ePeriodic
            continue;
        }
        cnt += (m_bndCondExpansions[i])->GetNcoeffs();
    }
 
    //----------------------------------
    // Solve trace problem: \Lambda = K^{-1} F
    // K is the HybridDGHelmBndLam matrix.
    //----------------------------------
    if (GloBndDofs - NumDirBCs > 0)
    {
        GlobalLinSysKey key(StdRegions::eHybridDGHelmBndLam, m_traceMap,
                            factors, varcoeff);
        GlobalLinSysSharedPtr LinSys = GetGlobalBndLinSys(key);
        LinSys->Solve(bndrhs, loclambda, m_traceMap);
 
        // For consistency with previous version put global
        // solution into m_trace->m_coeffs
        m_traceMap->LocalToGlobal(loclambda, m_trace->UpdateCoeffs());
    }
 
    //----------------------------------
    // Internal element solves
    //----------------------------------
    GlobalMatrixKey invHDGhelmkey(StdRegions::eInvHybridDGHelmholtz,
                                  NullAssemblyMapSharedPtr, factors, varcoeff);
 
    const DNekScalBlkMatSharedPtr &InvHDGHelm = GetBlockMatrix(invHDGhelmkey);
    DNekVec out(m_ncoeffs, outarray, eWrapper);
    Vmath::Zero(m_ncoeffs, outarray, 1);
 
    //  out =  u_f + u_lam = (*InvHDGHelm)*f + (LamtoU)*Lam
    out = (*InvHDGHelm) * F + (*HDGLamToU) * LocLambda;
 
    // Return empty GlobalLinSysKey
    return NullGlobalLinSysKey;
}
 
/* \brief This function evaluates the boundary conditions at a certain
 * time-level.
 *
 * Based on the boundary condition \f$g(\boldsymbol{x},t)\f$ evaluated
 * at a given time-level \a t, this function transforms the boundary
 * conditions onto the coefficients of the (one-dimensional) boundary
 * expansion. Depending on the type of boundary conditions, these
 * expansion coefficients are calculated in different ways:
 * - <b>Dirichlet boundary conditions</b><BR>
 *   In order to ensure global \f$C^0\f$ continuity of the spectral/hp
 *   approximation, the Dirichlet boundary conditions are projected onto
 *   the boundary expansion by means of a modified \f$C^0\f$ continuous
 *   Galerkin projection. This projection can be viewed as a collocation
 *   projection at the vertices, followed by an \f$L^2\f$ projection on
 *   the interior modes of the edges. The resulting coefficients
 *   \f$\boldsymbol{\hat{u}}^{\mathcal{D}}\f$ will be stored for the
 *   boundary expansion.
 * - <b>Neumann boundary conditions</b>
 *   In the discrete Galerkin formulation of the problem to be solved,
 *   the Neumann boundary conditions appear as the set of surface
 *   integrals: \f[\boldsymbol{\hat{g}}=\int_{\Gamma}
 *   \phi^e_n(\boldsymbol{x})g(\boldsymbol{x})d(\boldsymbol{x})\quad
 *   \forall n \f]
 *   As a result, it are the coefficients \f$\boldsymbol{\hat{g}}\f$
 *   that will be stored in the boundary expansion
 *
 * @param   time        The time at which the boundary conditions
 *                      should be evaluated.
 * @param   bndCondExpansions   List of boundary conditions.
 * @param   bndConditions   Information about the boundary conditions.
 *
 * This will only be undertaken for time dependent
 * boundary conditions unless time == 0.0 which is the
 * case when the method is called from the constructor.
 */
void DisContField::v_EvaluateBoundaryConditions(const NekDouble time,
                                                const std::string varName,
                                                const NekDouble x2_in,
                                                const NekDouble x3_in)
{
    int i;
    int npoints;
 
    MultiRegions::ExpListSharedPtr locExpList;
 
    for (i = 0; i < m_bndCondExpansions.size(); ++i)
    {
        if (m_bndCondExpansions[i] &&
            (time == 0.0 || m_bndConditions[i]->IsTimeDependent()))
        {
            m_bndCondBndWeight[i] = 1.0;
            locExpList            = m_bndCondExpansions[i];
 
            npoints = locExpList->GetNpoints();
            Array<OneD, NekDouble> x0(npoints, 0.0);
            Array<OneD, NekDouble> x1(npoints, 0.0);
            Array<OneD, NekDouble> x2(npoints, 0.0);
 
            locExpList->GetCoords(x0, x1, x2);
 
            if (x2_in != NekConstants::kNekUnsetDouble &&
                x3_in != NekConstants::kNekUnsetDouble)
            {
                Vmath::Fill(npoints, x2_in, x1, 1);
                Vmath::Fill(npoints, x3_in, x2, 1);
            }
            else if (x2_in != NekConstants::kNekUnsetDouble)
            {
                Vmath::Fill(npoints, x2_in, x2, 1);
            }
 
            // treat 1D expansions separately since we only
            // require an evaluation at a point rather than
            // any projections or inner products that are not
            // available in a PointExp
            if (m_expType == e1D)
            {
                if (m_bndConditions[i]->GetBoundaryConditionType() ==
                    SpatialDomains::eDirichlet)
                {
                    for (int n = 0; n < npoints; ++n)
                    {
                        m_bndCondExpansions[i]->SetCoeff(
                            n, (std::static_pointer_cast<
                                    SpatialDomains::DirichletBoundaryCondition>(
                                    m_bndConditions[i])
                                    ->m_dirichletCondition)
                                   ->Evaluate(x0[n], x1[n], x2[n], time));
                        m_bndCondExpansions[i]->SetPhys(
                            n, m_bndCondExpansions[i]->GetCoeff(n));
                    }
                }
                else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                         SpatialDomains::eNeumann)
                {
                    for (int n = 0; n < npoints; ++n)
                    {
                        m_bndCondExpansions[i]->SetCoeff(
                            n, (std::static_pointer_cast<
                                    SpatialDomains::NeumannBoundaryCondition>(
                                    m_bndConditions[i])
                                    ->m_neumannCondition)
                                   ->Evaluate(x0[n], x1[n], x2[n], time));
                    }
                }
                else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                         SpatialDomains::eRobin)
                {
                    for (int n = 0; n < npoints; ++n)
                    {
                        m_bndCondExpansions[i]->SetCoeff(
                            n, (std::static_pointer_cast<
                                    SpatialDomains::RobinBoundaryCondition>(
                                    m_bndConditions[i])
                                    ->m_robinFunction)
                                   ->Evaluate(x0[n], x1[n], x2[n], time));
                    }
                }
                else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                         SpatialDomains::ePeriodic)
                {
                    continue;
                }
                else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                         SpatialDomains::eNotDefined)
                {
                }
                else
                {
                    NEKERROR(ErrorUtil::efatal,
                             "This type of BC not implemented yet");
                }
            }
            else // 2D and 3D versions
            {
                if (m_bndConditions[i]->GetBoundaryConditionType() ==
                    SpatialDomains::eDirichlet)
                {
                    SpatialDomains::DirichletBCShPtr bcPtr =
                        std::static_pointer_cast<
                            SpatialDomains::DirichletBoundaryCondition>(
                            m_bndConditions[i]);
 
                    Array<OneD, NekDouble> valuesFile(npoints, 1.0),
                        valuesExp(npoints, 1.0);
 
                    string bcfilename = bcPtr->m_filename;
                    string exprbcs    = bcPtr->m_expr;
 
                    if (bcfilename != "")
                    {
                        locExpList->ExtractCoeffsFromFile(
                            bcfilename, bcPtr->GetComm(), varName,
                            locExpList->UpdateCoeffs());
                        locExpList->BwdTrans(locExpList->GetCoeffs(),
                                             locExpList->UpdatePhys());
                        valuesFile = locExpList->GetPhys();
                    }
 
                    if (exprbcs != "")
                    {
                        LibUtilities::EquationSharedPtr condition =
                            std::static_pointer_cast<
                                SpatialDomains::DirichletBoundaryCondition>(
                                m_bndConditions[i])
                                ->m_dirichletCondition;
 
                        condition->Evaluate(x0, x1, x2, time, valuesExp);
                    }
 
                    Vmath::Vmul(npoints, valuesExp, 1, valuesFile, 1,
                                locExpList->UpdatePhys(), 1);
                    locExpList->FwdTransBndConstrained(
                        locExpList->GetPhys(), locExpList->UpdateCoeffs());
                }
                else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                         SpatialDomains::eNeumann)
                {
                    SpatialDomains::NeumannBCShPtr bcPtr =
                        std::static_pointer_cast<
                            SpatialDomains::NeumannBoundaryCondition>(
                            m_bndConditions[i]);
                    string bcfilename = bcPtr->m_filename;
 
                    if (bcfilename != "")
                    {
                        locExpList->ExtractCoeffsFromFile(
                            bcfilename, bcPtr->GetComm(), varName,
                            locExpList->UpdateCoeffs());
                        locExpList->BwdTrans(locExpList->GetCoeffs(),
                                             locExpList->UpdatePhys());
                    }
                    else
                    {
                        LibUtilities::EquationSharedPtr condition =
                            std::static_pointer_cast<
                                SpatialDomains::NeumannBoundaryCondition>(
                                m_bndConditions[i])
                                ->m_neumannCondition;
                        condition->Evaluate(x0, x1, x2, time,
                                            locExpList->UpdatePhys());
                    }
 
                    locExpList->IProductWRTBase(locExpList->GetPhys(),
                                                locExpList->UpdateCoeffs());
                }
                else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                         SpatialDomains::eRobin)
                {
                    SpatialDomains::RobinBCShPtr bcPtr =
                        std::static_pointer_cast<
                            SpatialDomains::RobinBoundaryCondition>(
                            m_bndConditions[i]);
 
                    string bcfilename = bcPtr->m_filename;
 
                    if (bcfilename != "")
                    {
                        locExpList->ExtractCoeffsFromFile(
                            bcfilename, bcPtr->GetComm(), varName,
                            locExpList->UpdateCoeffs());
                        locExpList->BwdTrans(locExpList->GetCoeffs(),
                                             locExpList->UpdatePhys());
                    }
                    else
                    {
                        LibUtilities::EquationSharedPtr condition =
                            std::static_pointer_cast<
                                SpatialDomains::RobinBoundaryCondition>(
                                m_bndConditions[i])
                                ->m_robinFunction;
                        condition->Evaluate(x0, x1, x2, time,
                                            locExpList->UpdatePhys());
                    }
 
                    locExpList->IProductWRTBase(locExpList->GetPhys(),
                                                locExpList->UpdateCoeffs());
                }
                else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                         SpatialDomains::ePeriodic)
                {
                    continue;
                }
                else
                {
                    NEKERROR(ErrorUtil::efatal,
                             "This type of BC not implemented yet");
                }
            }
        }
    }
}
 
void DisContField::v_SetBCsToHomogeneous(void)
{
    int i;
    int npoints;
 
    MultiRegions::ExpListSharedPtr locExpList;
 
    for (i = 0; i < m_bndCondExpansions.size(); ++i)
    {
        locExpList = m_bndCondExpansions[i];
 
        npoints = locExpList->GetNpoints();
 
        // treat 1D expansions separately since we only
        // require an evaluation at a point rather than
        // any projections or inner products that are not
        // available in a PointExp
        if (m_expType == e1D)
        {
            NEKERROR(ErrorUtil::efatal,
                     "Function needs setting up for 1D boundary expansion");
        }
        else // 2D and 3D versions
        {
            if (m_bndConditions[i]->GetBoundaryConditionType() ==
                SpatialDomains::eDirichlet)
            {
                SpatialDomains::DirichletBCShPtr bcPtr =
                    std::static_pointer_cast<
                        SpatialDomains::DirichletBoundaryCondition>(
                        m_bndConditions[i]);
 
                Array<OneD, NekDouble> valuesExp(npoints, 0.0);
                bcPtr->m_filename = "";
                bcPtr->m_expr     = "0.0";
 
                bcPtr->m_dirichletCondition =
                    std::make_shared<LibUtilities::Equation>(
                        m_session->GetInterpreter(), bcPtr->m_expr);
 
                Vmath::Zero(locExpList->GetNcoeffs(),
                            locExpList->UpdateCoeffs(), 1);
                Vmath::Zero(locExpList->GetNpoints(), locExpList->UpdatePhys(),
                            1);
            }
            else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                     SpatialDomains::eNeumann)
            {
                SpatialDomains::NeumannBCShPtr bcPtr = std::static_pointer_cast<
                    SpatialDomains::NeumannBoundaryCondition>(
                    m_bndConditions[i]);
 
                bcPtr->m_filename = "";
                std::string eqn   = "0.0";
 
                bcPtr->m_neumannCondition =
                    std::make_shared<LibUtilities::Equation>(
                        m_session->GetInterpreter(), eqn);
 
                Vmath::Zero(locExpList->GetNcoeffs(),
                            locExpList->UpdateCoeffs(), 1);
                Vmath::Zero(locExpList->GetNpoints(), locExpList->UpdatePhys(),
                            1);
            }
            else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                     SpatialDomains::eRobin)
            {
                SpatialDomains::RobinBCShPtr bcPtr = std::static_pointer_cast<
                    SpatialDomains::RobinBoundaryCondition>(m_bndConditions[i]);
 
                bcPtr->m_filename = "";
                std::string eqn   = "0.0";
 
                bcPtr->m_robinFunction =
                    std::make_shared<LibUtilities::Equation>(
                        m_session->GetInterpreter(), eqn);
 
                Vmath::Zero(locExpList->GetNcoeffs(),
                            locExpList->UpdateCoeffs(), 1);
                Vmath::Zero(locExpList->GetNpoints(), locExpList->UpdatePhys(),
                            1);
            }
            else if (m_bndConditions[i]->GetBoundaryConditionType() ==
                     SpatialDomains::ePeriodic)
            {
                continue;
            }
            else
            {
                NEKERROR(ErrorUtil::efatal,
                         "This type of BC not implemented yet");
            }
        }
    }
}
 
/**
 * @brief Fill the weight with m_bndCondBndWeight.
 */
void DisContField::v_FillBwdWithBwdWeight(Array<OneD, NekDouble> &weightave,
                                          Array<OneD, NekDouble> &weightjmp)
{
    int cnt;
    int npts;
    int e = 0;
 
    // Fill boundary conditions into missing elements
    int id2 = 0;
 
    for (int n = cnt = 0; n < m_bndCondExpansions.size(); ++n)
    {
 
        if (m_bndConditions[n]->GetBoundaryConditionType() ==
            SpatialDomains::eDirichlet)
        {
            for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
            {
                npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                id2  = m_trace->GetPhys_Offset(
                    m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
                Vmath::Fill(npts, m_bndCondBndWeight[n], &weightave[id2], 1);
                Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
            }
 
            cnt += e;
        }
        else if (m_bndConditions[n]->GetBoundaryConditionType() ==
                     SpatialDomains::eNeumann ||
                 m_bndConditions[n]->GetBoundaryConditionType() ==
                     SpatialDomains::eRobin)
        {
            for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
            {
                npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
                id2  = m_trace->GetPhys_Offset(
                    m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
                Vmath::Fill(npts, m_bndCondBndWeight[n], &weightave[id2], 1);
                Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
            }
 
            cnt += e;
        }
        else if (m_bndConditions[n]->GetBoundaryConditionType() !=
                 SpatialDomains::ePeriodic)
        {
            NEKERROR(ErrorUtil::efatal,
                     "Method not set up for this boundary condition.");
        }
    }
}
 
// Set up a list of element ids and trace ids that link to the
// boundary conditions
void DisContField::v_GetBoundaryToElmtMap(Array<OneD, int> &ElmtID,
                                          Array<OneD, int> &TraceID)
{
 
    if (m_BCtoElmMap.size() == 0)
    {
        switch (m_expType)
        {
            case e1D:
            {
                map<int, int> VertGID;
                int i, n, id;
                int bid, cnt, Vid;
                int nbcs = 0;
                // Determine number of boundary condition expansions.
                for (i = 0; i < m_bndConditions.size(); ++i)
                {
                    nbcs += m_bndCondExpansions[i]->GetExpSize();
                }
 
                // make sure arrays are of sufficient length
                m_BCtoElmMap   = Array<OneD, int>(nbcs, -1);
                m_BCtoTraceMap = Array<OneD, int>(nbcs);
 
                // setup map of all global ids along boundary
                cnt = 0;
                for (n = 0; n < m_bndCondExpansions.size(); ++n)
                {
                    for (i = 0; i < m_bndCondExpansions[n]->GetExpSize(); ++i)
                    {
                        Vid = m_bndCondExpansions[n]
                                  ->GetExp(i)
                                  ->GetGeom()
                                  ->GetVid(0);
                        VertGID[Vid] = cnt++;
                    }
                }
 
                // Loop over elements and find verts that match;
                LocalRegions::ExpansionSharedPtr exp;
                for (cnt = n = 0; n < GetExpSize(); ++n)
                {
                    exp = (*m_exp)[n];
                    for (i = 0; i < exp->GetNverts(); ++i)
                    {
                        id = exp->GetGeom()->GetVid(i);
 
                        if (VertGID.count(id) > 0)
                        {
                            bid = VertGID.find(id)->second;
                            ASSERTL1(m_BCtoElmMap[bid] == -1,
                                     "Edge already set");
                            m_BCtoElmMap[bid]   = n;
                            m_BCtoTraceMap[bid] = i;
                            cnt++;
                        }
                    }
                }
                ASSERTL1(cnt == nbcs,
                         "Failed to visit all boundary condtiions");
            }
            break;
            case e2D:
            {
                map<int, int> globalIdMap;
                int i, n;
                int cnt;
                int nbcs = 0;
 
                // Populate global ID map (takes global geometry
                // ID to local expansion list ID).
                for (i = 0; i < GetExpSize(); ++i)
                {
                    globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
                }
 
                // Determine number of boundary condition expansions.
                for (i = 0; i < m_bndConditions.size(); ++i)
                {
                    nbcs += m_bndCondExpansions[i]->GetExpSize();
                }
 
                // Initialize arrays
                m_BCtoElmMap   = Array<OneD, int>(nbcs);
                m_BCtoTraceMap = Array<OneD, int>(nbcs);
 
                LocalRegions::Expansion1DSharedPtr exp1d;
                cnt = 0;
                for (n = 0; n < m_bndCondExpansions.size(); ++n)
                {
                    for (i = 0; i < m_bndCondExpansions[n]->GetExpSize();
                         ++i, ++cnt)
                    {
                        exp1d = m_bndCondExpansions[n]
                                    ->GetExp(i)
                                    ->as<LocalRegions::Expansion1D>();
 
                        // Use edge to element map from MeshGraph.
                        SpatialDomains::GeometryLinkSharedPtr tmp =
                            m_graph->GetElementsFromEdge(exp1d->GetGeom1D());
 
                        m_BCtoElmMap[cnt] =
                            globalIdMap[(*tmp)[0].first->GetGlobalID()];
                        m_BCtoTraceMap[cnt] = (*tmp)[0].second;
                    }
                }
            }
            break;
            case e3D:
            {
                map<int, int> globalIdMap;
                int i, n;
                int cnt;
                int nbcs = 0;
 
                // Populate global ID map (takes global geometry ID to local
                // expansion list ID).
                LocalRegions::Expansion3DSharedPtr exp3d;
                for (i = 0; i < GetExpSize(); ++i)
                {
                    globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
                }
 
                // Determine number of boundary condition expansions.
                for (i = 0; i < m_bndConditions.size(); ++i)
                {
                    if (m_bndCondExpansions[i])
                    {
                        nbcs += m_bndCondExpansions[i]->GetExpSize();
                    }
                }
 
                // Initialize arrays
                m_BCtoElmMap   = Array<OneD, int>(nbcs);
                m_BCtoTraceMap = Array<OneD, int>(nbcs);
 
                LocalRegions::Expansion2DSharedPtr exp2d;
                for (cnt = n = 0; n < m_bndCondExpansions.size(); ++n)
                {
                    for (i = 0; i < m_bndCondExpansions[n]->GetExpSize();
                         ++i, ++cnt)
                    {
                        exp2d = m_bndCondExpansions[n]
                                    ->GetExp(i)
                                    ->as<LocalRegions::Expansion2D>();
 
                        SpatialDomains::GeometryLinkSharedPtr tmp =
                            m_graph->GetElementsFromFace(exp2d->GetGeom2D());
                        m_BCtoElmMap[cnt] =
                            globalIdMap[tmp->at(0).first->GetGlobalID()];
                        m_BCtoTraceMap[cnt] = tmp->at(0).second;
                    }
                }
            }
            break;
            default:
                ASSERTL1(false, "Not setup for this expansion");
                break;
        }
    }
 
    ElmtID  = m_BCtoElmMap;
    TraceID = m_BCtoTraceMap;
}
 
void DisContField::v_GetBndElmtExpansion(int i,
                                         std::shared_ptr<ExpList> &result,
                                         const bool DeclareCoeffPhysArrays)
{
    int n, cnt;
    std::vector<unsigned int> eIDs;
 
    Array<OneD, int> ElmtID, TraceID;
    GetBoundaryToElmtMap(ElmtID, TraceID);
 
    // Skip other boundary regions
    for (cnt = n = 0; n < i; ++n)
    {
        cnt += m_bndCondExpansions[n]->GetExpSize();
    }
 
    // Populate eIDs with information from BoundaryToElmtMap
    for (n = 0; n < m_bndCondExpansions[i]->GetExpSize(); ++n)
    {
        eIDs.push_back(ElmtID[cnt + n]);
    }
 
    // Create expansion list
    result = MemoryManager<ExpList>::AllocateSharedPtr(
        *this, eIDs, DeclareCoeffPhysArrays, Collections::eNoCollection);
    // Copy phys and coeffs to new explist
    if (DeclareCoeffPhysArrays)
    {
        int nq;
        int offsetOld, offsetNew;
        Array<OneD, NekDouble> tmp1, tmp2;
        for (n = 0; n < result->GetExpSize(); ++n)
        {
            nq        = GetExp(ElmtID[cnt + n])->GetTotPoints();
            offsetOld = GetPhys_Offset(ElmtID[cnt + n]);
            offsetNew = result->GetPhys_Offset(n);
            Vmath::Vcopy(nq, tmp1 = GetPhys() + offsetOld, 1,
                         tmp2 = result->UpdatePhys() + offsetNew, 1);
            nq        = GetExp(ElmtID[cnt + n])->GetNcoeffs();
            offsetOld = GetCoeff_Offset(ElmtID[cnt + n]);
            offsetNew = result->GetCoeff_Offset(n);
            Vmath::Vcopy(nq, tmp1 = GetCoeffs() + offsetOld, 1,
                         tmp2 = result->UpdateCoeffs() + offsetNew, 1);
        }
    }
}
 
/**
 * @brief Reset this field, so that geometry information can be updated.
 */
void DisContField::v_Reset()
{
    ExpList::v_Reset();
 
    // Reset boundary condition expansions.
    for (int n = 0; n < m_bndCondExpansions.size(); ++n)
    {
        m_bndCondExpansions[n]->Reset();
    }
}
 
/**
 * Search through the edge expansions and identify which ones
 * have Robin/Mixed type boundary conditions. If find a Robin
 * boundary then store the edge id of the boundary condition
 * and the array of points of the physical space boundary
 * condition which are hold the boundary condition primitive
 * variable coefficient at the quatrature points
 *
 * \return std map containing the robin boundary condition
 * info using a key of the element id
 *
 * There is a next member to allow for more than one Robin
 * boundary condition per element
 */
map<int, RobinBCInfoSharedPtr> DisContField::v_GetRobinBCInfo(void)
{
    int i, cnt;
    map<int, RobinBCInfoSharedPtr> returnval;
    Array<OneD, int> ElmtID, TraceID;
    GetBoundaryToElmtMap(ElmtID, TraceID);
 
    for (cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
    {
        MultiRegions::ExpListSharedPtr locExpList;
 
        if (m_bndConditions[i]->GetBoundaryConditionType() ==
            SpatialDomains::eRobin)
        {
            int e, elmtid;
            Array<OneD, NekDouble> Array_tmp;
 
            locExpList = m_bndCondExpansions[i];
 
            int npoints = locExpList->GetNpoints();
            Array<OneD, NekDouble> x0(npoints, 0.0);
            Array<OneD, NekDouble> x1(npoints, 0.0);
            Array<OneD, NekDouble> x2(npoints, 0.0);
            Array<OneD, NekDouble> coeffphys(npoints);
 
            locExpList->GetCoords(x0, x1, x2);
 
            LibUtilities::EquationSharedPtr coeffeqn =
                std::static_pointer_cast<
                    SpatialDomains::RobinBoundaryCondition>(m_bndConditions[i])
                    ->m_robinPrimitiveCoeff;
 
            // evalaute coefficient
            coeffeqn->Evaluate(x0, x1, x2, 0.0, coeffphys);
 
            for (e = 0; e < locExpList->GetExpSize(); ++e)
            {
                RobinBCInfoSharedPtr rInfo =
                    MemoryManager<RobinBCInfo>::AllocateSharedPtr(
                        TraceID[cnt + e],
                        Array_tmp = coeffphys + locExpList->GetPhys_Offset(e));
 
                elmtid = ElmtID[cnt + e];
                // make link list if necessary
                if (returnval.count(elmtid) != 0)
                {
                    rInfo->next = returnval.find(elmtid)->second;
                }
                returnval[elmtid] = rInfo;
            }
        }
        cnt += m_bndCondExpansions[i]->GetExpSize();
    }
 
    return returnval;
}
 
/**
 * @brief Calculate the \f$ L^2 \f$ error of the \f$ Q_{\rm dir} \f$
 * derivative using the consistent DG evaluation of \f$ Q_{\rm dir} \f$.
 *
 * The solution provided is of the primative variation at the quadrature
 * points and the derivative is compared to the discrete derivative at
 * these points, which is likely to be undesirable unless using a much
 * higher number of quadrature points than the polynomial order used to
 * evaluate \f$ Q_{\rm dir} \f$.
 */
NekDouble DisContField::L2_DGDeriv(const int dir,
                                   const Array<OneD, const NekDouble> &coeffs,
                                   const Array<OneD, const NekDouble> &soln)
{
    int i, e, ncoeff_edge;
    Array<OneD, const NekDouble> tmp_coeffs;
    Array<OneD, NekDouble> out_d(m_ncoeffs), out_tmp;
 
    Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
        m_traceMap->GetElmtToTrace();
 
    StdRegions::Orientation edgedir;
 
    int cnt;
    int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
    Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda;
    m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
 
    edge_lambda = loc_lambda;
 
    // Calculate Q using standard DG formulation.
    for (i = cnt = 0; i < GetExpSize(); ++i)
    {
        // Probably a better way of setting up lambda than this.
        // Note cannot use PutCoeffsInToElmts since lambda space
        // is mapped during the solve.
        int nEdges = (*m_exp)[i]->GetGeom()->GetNumEdges();
        Array<OneD, Array<OneD, NekDouble>> edgeCoeffs(nEdges);
 
        for (e = 0; e < nEdges; ++e)
        {
            edgedir       = (*m_exp)[i]->GetTraceOrient(e);
            ncoeff_edge   = elmtToTrace[i][e]->GetNcoeffs();
            edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge);
            Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
            elmtToTrace[i][e]->SetCoeffsToOrientation(edgedir, edgeCoeffs[e],
                                                      edgeCoeffs[e]);
            edge_lambda = edge_lambda + ncoeff_edge;
        }
 
        (*m_exp)[i]->DGDeriv(dir, tmp_coeffs = coeffs + m_coeff_offset[i],
                             elmtToTrace[i], edgeCoeffs, out_tmp = out_d + cnt);
        cnt += (*m_exp)[i]->GetNcoeffs();
    }
 
    Array<OneD, NekDouble> phys(m_npoints);
    BwdTrans(out_d, phys);
    Vmath::Vsub(m_npoints, phys, 1, soln, 1, phys, 1);
    return L2(phys);
}
 
/**
 * @brief Evaluate HDG post-processing to increase polynomial order of
 * solution.
 *
 * This function takes the solution (assumed to be one order lower) in
 * physical space, and postprocesses at the current polynomial order by
 * solving the system:
 *
 * \f[
 * \begin{aligned}
 *   (\nabla w, \nabla u^*) &= (\nabla w, u), \\
 *   \langle \nabla u^*, 1 \rangle &= \langle \nabla u, 1 \rangle
 * \end{aligned}
 * \f]
 *
 * where \f$ u \f$ corresponds with the current solution as stored
 * inside #m_coeffs.
 *
 * @param outarray  The resulting field \f$ u^* \f$.
 */
void DisContField::EvaluateHDGPostProcessing(
    const Array<OneD, const NekDouble> &coeffs,
    Array<OneD, NekDouble> &outarray)
{
    int i, cnt, e, ncoeff_trace;
    Array<OneD, NekDouble> force, out_tmp, qrhs, qrhs1;
    Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
        m_traceMap->GetElmtToTrace();
 
    StdRegions::Orientation edgedir;
 
    int nq_elmt, nm_elmt;
    int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
    Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), trace_lambda;
    Array<OneD, NekDouble> tmp_coeffs;
    m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
 
    trace_lambda = loc_lambda;
 
    int dim = (m_expType == e2D) ? 2 : 3;
 
    int num_points[] = {0, 0, 0};
    int num_modes[]  = {0, 0, 0};
 
    // Calculate Q using standard DG formulation.
    for (i = cnt = 0; i < GetExpSize(); ++i)
    {
        nq_elmt = (*m_exp)[i]->GetTotPoints();
        nm_elmt = (*m_exp)[i]->GetNcoeffs();
        qrhs    = Array<OneD, NekDouble>(nq_elmt);
        qrhs1   = Array<OneD, NekDouble>(nq_elmt);
        force   = Array<OneD, NekDouble>(2 * nm_elmt);
        out_tmp = force + nm_elmt;
        LocalRegions::ExpansionSharedPtr ppExp;
 
        for (int j = 0; j < dim; ++j)
        {
            num_points[j] = (*m_exp)[i]->GetBasis(j)->GetNumPoints();
            num_modes[j]  = (*m_exp)[i]->GetBasis(j)->GetNumModes();
        }
 
        // Probably a better way of setting up lambda than this.  Note
        // cannot use PutCoeffsInToElmts since lambda space is mapped
        // during the solve.
        int nTraces = (*m_exp)[i]->GetNtraces();
        Array<OneD, Array<OneD, NekDouble>> traceCoeffs(nTraces);
 
        for (e = 0; e < (*m_exp)[i]->GetNtraces(); ++e)
        {
            edgedir        = (*m_exp)[i]->GetTraceOrient(e);
            ncoeff_trace   = elmtToTrace[i][e]->GetNcoeffs();
            traceCoeffs[e] = Array<OneD, NekDouble>(ncoeff_trace);
            Vmath::Vcopy(ncoeff_trace, trace_lambda, 1, traceCoeffs[e], 1);
            if (dim == 2)
            {
                elmtToTrace[i][e]->SetCoeffsToOrientation(
                    edgedir, traceCoeffs[e], traceCoeffs[e]);
            }
            else
            {
                (*m_exp)[i]
                    ->as<LocalRegions::Expansion3D>()
                    ->SetFaceToGeomOrientation(e, traceCoeffs[e]);
            }
            trace_lambda = trace_lambda + ncoeff_trace;
        }
 
        // creating orthogonal expansion (checking if we have quads or
        // triangles)
        LibUtilities::ShapeType shape = (*m_exp)[i]->DetShapeType();
        switch (shape)
        {
            case LibUtilities::eQuadrilateral:
            {
                const LibUtilities::PointsKey PkeyQ1(
                    num_points[0], LibUtilities::eGaussLobattoLegendre);
                const LibUtilities::PointsKey PkeyQ2(
                    num_points[1], LibUtilities::eGaussLobattoLegendre);
                LibUtilities::BasisKey BkeyQ1(LibUtilities::eOrtho_A,
                                              num_modes[0], PkeyQ1);
                LibUtilities::BasisKey BkeyQ2(LibUtilities::eOrtho_A,
                                              num_modes[1], PkeyQ2);
                SpatialDomains::QuadGeom *qGeom =
                    dynamic_cast<SpatialDomains::QuadGeom *>(
                        (*m_exp)[i]->GetGeom());
                ppExp = MemoryManager<LocalRegions::QuadExp>::AllocateSharedPtr(
                    BkeyQ1, BkeyQ2, qGeom);
            }
            break;
            case LibUtilities::eTriangle:
            {
                const LibUtilities::PointsKey PkeyT1(
                    num_points[0], LibUtilities::eGaussLobattoLegendre);
                const LibUtilities::PointsKey PkeyT2(
                    num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
                LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A,
                                              num_modes[0], PkeyT1);
                LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B,
                                              num_modes[1], PkeyT2);
                SpatialDomains::TriGeom *tGeom =
                    dynamic_cast<SpatialDomains::TriGeom *>(
                        (*m_exp)[i]->GetGeom());
                ppExp = MemoryManager<LocalRegions::TriExp>::AllocateSharedPtr(
                    BkeyT1, BkeyT2, tGeom);
            }
            break;
            case LibUtilities::eHexahedron:
            {
                const LibUtilities::PointsKey PkeyH1(
                    num_points[0], LibUtilities::eGaussLobattoLegendre);
                const LibUtilities::PointsKey PkeyH2(
                    num_points[1], LibUtilities::eGaussLobattoLegendre);
                const LibUtilities::PointsKey PkeyH3(
                    num_points[2], LibUtilities::eGaussLobattoLegendre);
                LibUtilities::BasisKey BkeyH1(LibUtilities::eOrtho_A,
                                              num_modes[0], PkeyH1);
                LibUtilities::BasisKey BkeyH2(LibUtilities::eOrtho_A,
                                              num_modes[1], PkeyH2);
                LibUtilities::BasisKey BkeyH3(LibUtilities::eOrtho_A,
                                              num_modes[2], PkeyH3);
                SpatialDomains::HexGeom *hGeom =
                    dynamic_cast<SpatialDomains::HexGeom *>(
                        (*m_exp)[i]->GetGeom());
                ppExp = MemoryManager<LocalRegions::HexExp>::AllocateSharedPtr(
                    BkeyH1, BkeyH2, BkeyH3, hGeom);
            }
            break;
            case LibUtilities::eTetrahedron:
            {
                const LibUtilities::PointsKey PkeyT1(
                    num_points[0], LibUtilities::eGaussLobattoLegendre);
                const LibUtilities::PointsKey PkeyT2(
                    num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
                const LibUtilities::PointsKey PkeyT3(
                    num_points[2], LibUtilities::eGaussRadauMAlpha2Beta0);
                LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A,
                                              num_modes[0], PkeyT1);
                LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B,
                                              num_modes[1], PkeyT2);
                LibUtilities::BasisKey BkeyT3(LibUtilities::eOrtho_C,
                                              num_modes[2], PkeyT3);
                SpatialDomains::TetGeom *tGeom =
                    dynamic_cast<SpatialDomains::TetGeom *>(
                        (*m_exp)[i]->GetGeom());
                ppExp = MemoryManager<LocalRegions::TetExp>::AllocateSharedPtr(
                    BkeyT1, BkeyT2, BkeyT3, tGeom);
            }
            break;
            default:
                NEKERROR(ErrorUtil::efatal,
                         "Wrong shape type, HDG postprocessing is not "
                         "implemented");
        };
 
        // DGDeriv
        // (d/dx w, d/dx q_0)
        (*m_exp)[i]->DGDeriv(0, tmp_coeffs = coeffs + m_coeff_offset[i],
                             elmtToTrace[i], traceCoeffs, out_tmp);
        (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
        ppExp->IProductWRTDerivBase(0, qrhs, force);
 
        // + (d/dy w, d/dy q_1)
        (*m_exp)[i]->DGDeriv(1, tmp_coeffs = coeffs + m_coeff_offset[i],
                             elmtToTrace[i], traceCoeffs, out_tmp);
 
        (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
        ppExp->IProductWRTDerivBase(1, qrhs, out_tmp);
 
        Vmath::Vadd(nm_elmt, force, 1, out_tmp, 1, force, 1);
 
        // determine force[0] = (1,u)
        (*m_exp)[i]->BwdTrans(tmp_coeffs = coeffs + m_coeff_offset[i], qrhs);
        force[0] = (*m_exp)[i]->Integral(qrhs);
 
        // multiply by inverse Laplacian matrix
        // get matrix inverse
        LocalRegions::MatrixKey lapkey(StdRegions::eInvLaplacianWithUnityMean,
                                       ppExp->DetShapeType(), *ppExp);
        DNekScalMatSharedPtr lapsys = ppExp->GetLocMatrix(lapkey);
 
        NekVector<NekDouble> in(nm_elmt, force, eWrapper);
        NekVector<NekDouble> out(nm_elmt);
 
        out = (*lapsys) * in;
 
        // Transforming back to modified basis
        Array<OneD, NekDouble> work(nq_elmt);
        ppExp->BwdTrans(out.GetPtr(), work);
        (*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray + m_coeff_offset[i]);
    }
}
 
void DisContField::v_GetLocTraceFromTracePts(
    const Array<OneD, const NekDouble> &Fwd,
    const Array<OneD, const NekDouble> &Bwd,
    Array<OneD, NekDouble> &locTraceFwd, Array<OneD, NekDouble> &locTraceBwd)
{
    if (locTraceFwd != NullNekDouble1DArray)
    {
        m_locTraceToTraceMap->InterpLocTracesToTraceTranspose(0, Fwd,
                                                              locTraceFwd);
    }
 
    if (locTraceBwd != NullNekDouble1DArray)
    {
        m_locTraceToTraceMap->InterpLocTracesToTraceTranspose(1, Bwd,
                                                              locTraceBwd);
    }
}
 
void DisContField::v_AddTraceIntegralToOffDiag(
    const Array<OneD, const NekDouble> &FwdFlux,
    const Array<OneD, const NekDouble> &BwdFlux,
    Array<OneD, NekDouble> &outarray)
{
    Array<OneD, NekDouble> FCoeffs(m_trace->GetNcoeffs());
 
    m_trace->IProductWRTBase(FwdFlux, FCoeffs);
    m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1, FCoeffs, outarray);
    m_trace->IProductWRTBase(BwdFlux, FCoeffs);
    m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0, FCoeffs, outarray);
}
 
/// \brief Rotate the slice [offset, offset + npts) of the
/// 3D vector field in Bwd around the axis 'dir' by 'angle'.
void DisContField::Rotate(Array<OneD, Array<OneD, NekDouble>> &Bwd,
                          const int dir, const NekDouble angle,
                          const int offset, const int npts)
{
    // Precompute trigonometric values for improved precision
    const NekDouble cosA = cos(angle);
    const NekDouble sinA = sin(angle);
 
    switch (dir)
    {
        // ----------------------------------------------------------
        // Rotate around the x-axis by 'angle':
        //  y' =  y*cos(angle) - z*sin(angle)
        //  z' =  y*sin(angle) + z*cos(angle)
        //  x' =  x (unchanged)
        // ----------------------------------------------------------
        case 0:
        {
            for (int i = offset; i < offset + npts; ++i)
            {
                NekDouble tmpY = Bwd[2][i];
                NekDouble tmpZ = Bwd[3][i];
 
                NekDouble yrot = cosA * tmpY - sinA * tmpZ;
                NekDouble zrot = sinA * tmpY + cosA * tmpZ;
 
                Bwd[2][i] = yrot;
                Bwd[3][i] = zrot;
            }
        }
        break;
 
        // ----------------------------------------------------------
        // Rotate around the y-axis by 'angle':
        //  z' =  z*cos(angle) - x*sin(angle)
        //  x' =  z*sin(angle) + x*cos(angle)
        //  y' =  y (unchanged)
        // ----------------------------------------------------------
        case 1:
        {
            for (int i = offset; i < offset + npts; ++i)
            {
                NekDouble tmpX = Bwd[1][i];
                NekDouble tmpZ = Bwd[3][i];
 
                NekDouble zrot = cosA * tmpZ - sinA * tmpX;
                NekDouble xrot = sinA * tmpZ + cosA * tmpX;
 
                Bwd[1][i] = xrot;
                Bwd[3][i] = zrot;
            }
        }
        break;
 
        // ----------------------------------------------------------
        // Rotate around the z-axis by 'angle':
        //  x' =  x*cos(angle) - y*sin(angle)
        //  y' =  x*sin(angle) + y*cos(angle)
        //  z' =  z (unchanged)
        // ----------------------------------------------------------
        case 2:
        {
            for (int i = offset; i < offset + npts; ++i)
            {
                NekDouble tmpX = Bwd[1][i];
                NekDouble tmpY = Bwd[2][i];
 
                NekDouble xrot = cosA * tmpX - sinA * tmpY;
                NekDouble yrot = sinA * tmpX + cosA * tmpY;
 
                Bwd[1][i] = xrot;
                Bwd[2][i] = yrot;
            }
        }
        break;
 
        default:
            NEKERROR(ErrorUtil::efatal,
                     "Rotate() axis must be 0 (x), 1 (y), or 2 (z).");
            break;
    }
}
 
/// \brief Rotate the slice [offset, offset + npts) of the
/// 3D vector field in Bwd around the axis 'dir' by 'angle'.
void DisContField::DeriveRotate(TensorOfArray3D<NekDouble> &Bwd, const int dir,
                                const NekDouble angle, const int offset,
                                const int npts)
{
    // Precompute trigonometric values using correct angle input
    const NekDouble cosA = cos(angle);
    const NekDouble sinA = sin(angle);
 
    // Define rotation matrix R (3x3) using Array<OneD, NekDouble> in
    // **column-major order**
    Array<OneD, NekDouble> R(9, 0.0);
 
    if (dir == 0) // Rotation around X-axis
    {
        R[0] = 1.0;
        R[3] = 0.0;
        R[6] = 0.0;
        R[1] = 0.0;
        R[4] = cosA;
        R[7] = -sinA;
        R[2] = 0.0;
        R[5] = sinA;
        R[8] = cosA;
    }
    else if (dir == 1) // Rotation around Y-axis
    {
        R[0] = cosA;
        R[3] = 0.0;
        R[6] = sinA;
        R[1] = 0.0;
        R[4] = 1.0;
        R[7] = 0.0;
        R[2] = -sinA;
        R[5] = 0.0;
        R[8] = cosA;
    }
    else if (dir == 2) // Rotation around Z-axis
    {
        R[0] = cosA;
        R[3] = -sinA;
        R[6] = 0.0;
        R[1] = sinA;
        R[4] = cosA;
        R[7] = 0.0;
        R[2] = 0.0;
        R[5] = 0.0;
        R[8] = 1.0;
    }
    else
    {
        NEKERROR(
            ErrorUtil::efatal,
            "Invalid rotation axis for first-order derivative transformation.");
    }
 
    for (int p = offset; p < offset + npts; ++p)
    {
        // Allocate gradient tensor (3x3) in **column-major order**
        Array<OneD, NekDouble> gradU(9, 0.0);
 
        // Load original velocity gradient tensor into **column-major** format
        gradU[0] = Bwd[0][1][p]; // ∂u_x/∂x
        gradU[3] = Bwd[1][1][p]; // ∂u_x/∂y
        gradU[6] = Bwd[2][1][p]; // ∂u_x/∂z
 
        gradU[1] = Bwd[0][2][p]; // ∂u_y/∂x
        gradU[4] = Bwd[1][2][p]; // ∂u_y/∂y
        gradU[7] = Bwd[2][2][p]; // ∂u_y/∂z
 
        gradU[2] = Bwd[0][3][p]; // ∂u_z/∂x
        gradU[5] = Bwd[1][3][p]; // ∂u_z/∂y
        gradU[8] = Bwd[2][3][p]; // ∂u_z/∂z
 
        // Allocate intermediate matrix R * gradU in **column-major order**
        Array<OneD, NekDouble> R_gradU(9, 0.0);
 
        // Perform matrix multiplication R * gradU using Blas::Dgemm in
        // **column-major order**
        Blas::Dgemm('N', 'N', 3, 3, 3, 1.0, R.data(), 3, gradU.data(), 3, 0.0,
                    R_gradU.data(), 3);
 
        // Allocate final rotated gradient (R * gradU) * R^T in **column-major
        // order**
        Array<OneD, NekDouble> gradU_rotated(9, 0.0);
 
        // Perform matrix multiplication (R * gradU) * R^T using Blas::Dgemm in
        // **column-major order**
        Blas::Dgemm('N', 'T', 3, 3, 3, 1.0, R_gradU.data(), 3, R.data(), 3, 0.0,
                    gradU_rotated.data(), 3);
 
        // Store rotated gradients back in Bwd in **column-major order**
        Bwd[0][1][p] = gradU_rotated[0]; // ∂u_x'/∂x'
        Bwd[1][1][p] = gradU_rotated[3]; // ∂u_x'/∂y'
        Bwd[2][1][p] = gradU_rotated[6]; // ∂u_x'/∂z'
 
        Bwd[0][2][p] = gradU_rotated[1]; // ∂u_y'/∂x'
        Bwd[1][2][p] = gradU_rotated[4]; // ∂u_y'/∂y'
        Bwd[2][2][p] = gradU_rotated[7]; // ∂u_y'/∂z'
 
        Bwd[0][3][p] = gradU_rotated[2]; // ∂u_z'/∂x'
        Bwd[1][3][p] = gradU_rotated[5]; // ∂u_z'/∂y'
        Bwd[2][3][p] = gradU_rotated[8]; // ∂u_z'/∂z'
    }
}
 
} // namespace Nektar::MultiRegions