LocTraceToTraceMap.h

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: LocTraceToTraceMap.h
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Local Trace to general trace map information
//
///////////////////////////////////////////////////////////////////////////////
 
#ifndef MULTIREGIONS_LOCTRACETOTRACEMAP_H
#define MULTIREGIONS_LOCTRACETOTRACEMAP_H
 
#include <tuple>
 
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Foundations/Points.h>
#include <LibUtilities/LinearAlgebra/NekTypeDefs.hpp>
#include <MultiRegions/MultiRegionsDeclspec.h>
 
namespace Nektar
{
namespace LocalRegions
{
class Expansion;
typedef std::shared_ptr<Expansion> ExpansionSharedPtr;
} // namespace LocalRegions
 
namespace MultiRegions
{
class ExpList;
typedef std::shared_ptr<ExpList> ExpListSharedPtr;
 
enum InterpLocTraceToTrace
{
    eNoInterp,
    eInterpDir0,
    eInterpEndPtDir0,
    eInterpDir1,
    eInterpEndPtDir1,
    eInterpBothDirs,
    eInterpEndPtDir0InterpDir1
};
 
/**
 * @brief Map holding points distributions required for interpolation of local
 * traces onto global trace in two and three dimensions.
 *
 * The tuple contains four PointsKey entries. All four are used in 3D, but for
 * 2D only two are used. The first two denote the local elemental points
 * distributions, the latter the target trace distributions. In the 2D case,
 * unused trace directions are set to LibUtilities::eNoPointsType.
 */
typedef std::tuple<LibUtilities::PointsKey, LibUtilities::PointsKey,
                   LibUtilities::PointsKey, LibUtilities::PointsKey>
    TraceInterpPoints;
 
struct cmpop
{
    bool operator()(TraceInterpPoints const &a,
                    TraceInterpPoints const &b) const
    {
        if (std::get<0>(a) < std::get<0>(b))
        {
            return true;
        }
 
        if (std::get<0>(b) < std::get<0>(a))
        {
            return false;
        }
 
        if (std::get<1>(a) < std::get<1>(b))
        {
            return true;
        }
        if (std::get<1>(b) < std::get<1>(a))
        {
            return false;
        }
 
        if (std::get<2>(a) < std::get<2>(b))
        {
            return true;
        }
 
        if (std::get<2>(b) < std::get<2>(a))
        {
            return false;
        }
 
        if (std::get<3>(a) < std::get<3>(b))
        {
            return true;
        }
 
        return false;
    }
};
 
struct TraceInterpEssential
{
    /// useful for allocate workspace
    int m_maxTraceSize;
    /// All execept m_interpTraceIndex are shared for collections
    /// subscript of m_interpTrace, referenced by element and trace id
    /// the storage is ordered by (e,n)
    /// if width > 1 then this array should better be ordered by (e, g, k)
    /// The argument is whether we need to reorder the element trace to group
    /// them by interpolation type... The current decision is to keep them
    /// in natrual order.
    Array<OneD, Array<OneD, int>> m_interpTraceIndex;
    /// A mapping holding the type of interpolation needed for each local trace.
    /// Dimension 0 holds forward traces, dimension 1 backward.
    Array<OneD, Array<OneD, InterpLocTraceToTrace>> m_interpTrace;
    /// Interpolation matrices for either 2D edges or first coordinate of 3D
    /// face.
    Array<OneD, Array<OneD, DNekMatSharedPtr>> m_interpTraceI0;
    /// Interpolation matrices for the second coordinate of 3D face, not used in
    /// 2D.
    Array<OneD, Array<OneD, DNekMatSharedPtr>> m_interpTraceI1;
    /// Interpolation matrices for either 2D edges or first coordinate
    /// of 3D face using going "from' to 'to' points (i.e. the reverse
    /// of other techniques)
    Array<OneD, Array<OneD, DNekMatSharedPtr>> m_interpFromTraceI0;
    /// Interpolation matrices for either 2D edges or first coordinate
    /// of 3D face using going "from' to 'to' points (i.e. the reverse
    /// of other techniques)
    Array<OneD, Array<OneD, DNekMatSharedPtr>> m_interpFromTraceI1;
    /// Interpolation points key distributions to each of the local to global
    /// mappings.
    Array<OneD, Array<OneD, TraceInterpPoints>> m_interpPoints;
    /// Mapping to hold first coordinate direction endpoint interpolation, which
    /// can be more optimal if using Gauss-Radau distribution for triangles
    Array<OneD, Array<OneD, Array<OneD, NekDouble>>> m_interpEndPtI0;
    /// Mapping to hold second coordinate direction endpoint interpolation,
    /// which can be more optimal if using Gauss-Radau distribution for
    /// triangles
    Array<OneD, Array<OneD, Array<OneD, NekDouble>>> m_interpEndPtI1;
    /// Number of edges/faces on a 2D/3D element that require interpolation.
    Array<OneD, Array<OneD, int>> m_interpNtraces;
};
 
struct TraceFieldMapEssential
{
    /// stores the maps to extract local trace from element
    /// each collection contains a set of different maps;
    /// the storage is ordered is (e) and will not be changed once created
    Array<OneD, Array<OneD, int>> m_locTracePhysToElmtMaps;
    /// stores the maps to reorder locTrace To Trace (in trace size)
    /// the storage will not change once created
    Array<OneD, Array<OneD, int>> m_orientationMaps;
    /// The orientation + type of each local element trace, just store the int
    /// each collection contains a different array of size (pCollExp * Ntraces)
    /// the storage is ordered is (e,n) and if width > 1 then this array should
    /// be ordered by (e, g, v). The padding entry will have default value of 0
    Array<OneD, Array<OneD, int>> m_orientationIds;
    /// Returns the start pos of trace phys by giving the local elmt trace id
    /// each collection contains a different array of size (pCollExp * Ntraces)
    /// originally in TraceToFieldMap, now moves to this class
    /// can be created on a single Array by offset
    /// the storage is ordered is (e,n) and if width > 1 then this array should
    /// be ordered by (e, g, v). The padding entry will have default value of
    /// ntracePts - implying the trace array must have bigger size to hold the
    /// unused data
    Array<OneD, Array<OneD, int>> m_locToTracePhysOffset;
    Array<OneD, Array<OneD, int>> m_locToTraceId;
    /// bool to indicate if this loc trace is left adjacent to trace
    /// each collection contains a different array of size (pCollExp * Ntraces)
    /// can be created by offset LeftAdjacents (which should be moved from
    /// DisContField to this field?) the storage is ordered is (e,n) and if
    /// width > 1 then this array should be ordered by (e, g, v). The padding
    /// entry will have default value of true
    Array<OneD, Array<OneD, bool>> m_isLocTraceLeftAdjacent;
};
 
/**
 * @brief A helper class to deal with trace operations in the discontinuous
 * Galerkin code.
 *
 * This class sets up a number of mappings to deal with operations that take the
 * "local trace" of an expansion list -- i.e. the concatenation of all elemental
 * facets -- to the "global trace" -- where the duplicate facets between
 * connected elements have been removed.
 *
 * <pre>
 * Elements:      Local trace:              Global trace:
 * +----+----+    + +---+ +   + +---+ +     + +---+ + +---+ +
 * |    |    |    |       |   |       |     |       |       |
 * |    |    |    |       |   |       |     |       |       |
 * +----+----+    + +---+ +   + +---+ +     + +---+ + +---+ +
 * </pre>
 *
 * There are a number of mappings that are required that this class provides
 * maps for:
 *
 *   - Extracting the local trace from the elements in physical space
 *   - Dealing with cases that need interpolation (i.e. where the polynomial
 *     order varies between two elements)
 *   - Adding global trace contributions back into the elements.
 *
 * These are documented in the member variables and class functions.
 */
class LocTraceToTraceMap
{
public:
    // Constructor
    MULTI_REGIONS_EXPORT LocTraceToTraceMap(
        const ExpList &locExp, const ExpListSharedPtr &trace,
        const Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
            &elmtToTrace,
        const std::vector<bool> &LeftAdjacents);
 
    // Destructor
    MULTI_REGIONS_EXPORT virtual ~LocTraceToTraceMap();
 
    MULTI_REGIONS_EXPORT void LocTracesFromField(
        const Array<OneD, const NekDouble> &field,
        Array<OneD, NekDouble> faces);
 
    MULTI_REGIONS_EXPORT void FwdLocTracesFromField(
        const Array<OneD, const NekDouble> &field,
        Array<OneD, NekDouble> faces);
    MULTI_REGIONS_EXPORT void AddLocTracesToField(
        const Array<OneD, const NekDouble> &faces,
        Array<OneD, NekDouble> &field);
 
    MULTI_REGIONS_EXPORT void ReshuffleLocTracesForInterp(
        const int dir, const Array<OneD, const NekDouble> &loctraces,
        Array<OneD, NekDouble> reshuffle);
 
    MULTI_REGIONS_EXPORT void UnshuffleLocTraces(
        const int dir, const Array<OneD, const NekDouble> &loctraces,
        Array<OneD, NekDouble> unshuffle);
 
    MULTI_REGIONS_EXPORT void InterpLocTracesToTrace(
        const int dir, const Array<OneD, const NekDouble> &loctraces,
        Array<OneD, NekDouble> &traces);
 
    MULTI_REGIONS_EXPORT inline void InterpLocEdgesToTrace(
        const int dir, const Array<OneD, const NekDouble> &locfaces,
        Array<OneD, NekDouble> &edges);
 
    MULTI_REGIONS_EXPORT void InterpLocFacesToTrace(
        const int dir, const Array<OneD, const NekDouble> &locfaces,
        Array<OneD, NekDouble> faces);
 
    MULTI_REGIONS_EXPORT void InterpLocTracesToTraceTranspose(
        const int dir, const Array<OneD, const NekDouble> &traces,
        Array<OneD, NekDouble> &loctraces);
 
    /// Transpose of interp local edges to Trace methods
    MULTI_REGIONS_EXPORT void InterpLocEdgesToTraceTranspose(
        const int dir, const Array<OneD, const NekDouble> &edges,
        Array<OneD, NekDouble> &locedges);
 
    /// Transpose of interp local faces to Trace methods
    MULTI_REGIONS_EXPORT void InterpLocFacesToTraceTranspose(
        const int dir, const Array<OneD, const NekDouble> &traces,
        Array<OneD, NekDouble> &loctraces);
 
    MULTI_REGIONS_EXPORT void InterpTraceToLocTrace(
        const int dir, const Array<OneD, NekDouble> &traces,
        Array<OneD, NekDouble> &loctraces);
 
    MULTI_REGIONS_EXPORT inline void InterpTraceToLocEdges(
        const int dir, const Array<OneD, const NekDouble> &edges,
        Array<OneD, NekDouble> &locedges);
 
    MULTI_REGIONS_EXPORT void InterpTraceToLocFaces(
        const int dir, const Array<OneD, const NekDouble> &faces,
        Array<OneD, NekDouble> &locfaces);
 
    MULTI_REGIONS_EXPORT void AddTraceCoeffsToFieldCoeffs(
        const Array<OneD, const NekDouble> &trace,
        Array<OneD, NekDouble> &field);
 
    MULTI_REGIONS_EXPORT void AddTraceCoeffsToFieldCoeffs(
        const int dir, const Array<OneD, const NekDouble> &race,
        Array<OneD, NekDouble> &field);
 
    MULTI_REGIONS_EXPORT void CalcLocTracePhysToTraceIDMap(
        const ExpListSharedPtr &tracelist, const int ndim);
 
    MULTI_REGIONS_EXPORT void CalcLocTracePhysToTraceIDMap_2D(
        const ExpListSharedPtr &tracelist);
 
    MULTI_REGIONS_EXPORT void CalcLocTracePhysToTraceIDMap_3D(
        const ExpListSharedPtr &tracelist);
 
    /**
     * @brief Return the number of `forward' local trace points.
     */
    MULTI_REGIONS_EXPORT inline int GetNFwdLocTracePts()
    {
        return m_nFwdLocTracePts;
    }
 
    /**
     * @brief Return the number of local trace points.
     */
    MULTI_REGIONS_EXPORT inline int GetNLocTracePts()
    {
        return m_nLocTracePts;
    }
 
    MULTI_REGIONS_EXPORT inline const Array<OneD, const Array<OneD, bool>> &
    GetLeftRightAdjacentExpFlag() const
    {
        return m_leftRightAdjacentExpFlag;
    }
 
    MULTI_REGIONS_EXPORT inline const Array<OneD, const Array<OneD, int>> &
    GetLeftRightAdjacentExpId() const
    {
        return m_leftRightAdjacentExpId;
    }
 
    MULTI_REGIONS_EXPORT inline const Array<
        OneD, const Array<OneD, Array<OneD, int>>> &
    GetTraceCoeffToLeftRightExpCoeffMap() const
    {
        return m_traceCoeffToLeftRightExpCoeffMap;
    }
 
    MULTI_REGIONS_EXPORT inline const Array<
        OneD, const Array<OneD, Array<OneD, int>>> &
    GetTraceCoeffToLeftRightExpCoeffSign() const
    {
        return m_traceCoeffToLeftRightExpCoeffSign;
    }
 
    MULTI_REGIONS_EXPORT inline const Array<OneD, int> &GetElemNeighbsNumb()
        const
    {
        return m_ElemNeighbsNumb;
    }
 
    MULTI_REGIONS_EXPORT inline const Array<OneD, const Array<OneD, int>> &
    GetElemNeighbsId() const
    {
        return m_ElemNeighbsId;
    }
 
    MULTI_REGIONS_EXPORT inline const Array<OneD, const Array<OneD, int>> &
    GetLocTracephysToTraceIDMap() const
    {
        return m_LocTracephysToTraceIDMap;
    }
 
    MULTI_REGIONS_EXPORT inline void SetLocTracePhysToTraceIDMap(
        const Array<OneD, Array<OneD, int>> &inarray)
    {
        m_LocTracephysToTraceIDMap = inarray;
    }
 
    MULTI_REGIONS_EXPORT inline const Array<OneD, const int> &
    GetLocTraceToFieldMap() const
    {
        return m_locTraceToFieldMap;
    }
 
    MULTI_REGIONS_EXPORT void TraceLocToElmtLocCoeffMap(
        const ExpList &locExp, const ExpListSharedPtr &trace);
 
    MULTI_REGIONS_EXPORT const TraceInterpEssential &GetTraceInterpEssential(
        const int cid);
 
    MULTI_REGIONS_EXPORT const TraceFieldMapEssential &GetTraceFieldMapEssential(
        const int cid);
 
private:
    /// Expansion Dimension we have setup for trace mapping.
    int m_expdim;
    /// The number of forward trace points. A local trace element is `forward'
    /// if it is the side selected for the global trace.
    int m_nFwdLocTracePts;
    /// The number of local trace points.
    int m_nLocTracePts;
    /// The number of global trace points.
    int m_nTracePts;
    /// A mapping from the local elemental trace points, arranged as
    /// all forwards traces followed by backwards traces, to elemental
    /// storage.
    Array<OneD, int> m_locTraceToFieldMap;
    /// A mapping from the local elemental trace points, arranged as
    /// all forwards traces followed by backwards traces, to elemental
    /// storage.
    Array<OneD, Array<OneD, int>> m_locTraceToElmtTraceMap;
    /// A mapping from local trace points to the global trace. Dimension 0 holds
    /// forward traces, dimension 1 backward.
    Array<OneD, Array<OneD, int>> m_locInterpTraceToTraceMap;
    /// A mapping holding the type of interpolation needed for each local trace.
    /// Dimension 0 holds forward traces, dimension 1 backward.
    Array<OneD, Array<OneD, InterpLocTraceToTrace>> m_interpTrace;
    /// Interpolation matrices for either 2D edges or first coordinate of 3D
    /// face.
    Array<OneD, Array<OneD, DNekMatSharedPtr>> m_interpTraceI0;
    /// Interpolation matrices for the second coordinate of 3D face, not used in
    /// 2D.
    Array<OneD, Array<OneD, DNekMatSharedPtr>> m_interpTraceI1;
    /// Interpolation matrices for either 2D edges or first coordinate
    /// of 3D face using going "from' to 'to' points (i.e. the reverse
    /// of other techniques)
    Array<OneD, Array<OneD, DNekMatSharedPtr>> m_interpFromTraceI0;
    /// Interpolation matrices for either 2D edges or first coordinate
    /// of 3D face using going "from' to 'to' points (i.e. the reverse
    /// of other techniques)
    Array<OneD, Array<OneD, DNekMatSharedPtr>> m_interpFromTraceI1;
    /// Interpolation points key distributions to each of the local to global
    /// mappings.
    Array<OneD, Array<OneD, TraceInterpPoints>> m_interpPoints;
    /// Mapping to hold first coordinate direction endpoint interpolation, which
    /// can be more optimal if using Gauss-Radau distribution for triangles
    Array<OneD, Array<OneD, Array<OneD, NekDouble>>> m_interpEndPtI0;
    /// Mapping to hold second coordinate direction endpoint interpolation,
    /// which can be more optimal if using Gauss-Radau distribution for
    /// triangles
    Array<OneD, Array<OneD, Array<OneD, NekDouble>>> m_interpEndPtI1;
    /// Number of edges/faces on a 2D/3D element that require interpolation.
    Array<OneD, Array<OneD, int>> m_interpNtraces;
    /// Number of forwards/backwards trace coefficients.
    int m_nTraceCoeffs[2];
    /// Mapping from forwards/backwards trace coefficients to elemental
    /// coefficient storage.
    Array<OneD, Array<OneD, int>> m_traceCoeffsToElmtMap;
    /// Mapping from forwards/backwards trace coefficients to the position of
    /// the trace element in global storage.
    Array<OneD, Array<OneD, int>> m_traceCoeffsToElmtTrace;
    /// Sign array for mapping from forwards/backwards trace coefficients to
    /// local trace storage.
    Array<OneD, Array<OneD, int>> m_traceCoeffsToElmtSign;
    /// Flag indicates whether the expansion that are the left & right adjacent
    /// to current trace exists.
    Array<OneD, Array<OneD, bool>> m_leftRightAdjacentExpFlag;
    /// The expansion ID that are the left & right adjacent to current trace.
    Array<OneD, Array<OneD, int>> m_leftRightAdjacentExpId;
    /// The map of every coeff from current trace to the left & right adjacent
    /// expasion coeffs.
    Array<OneD, Array<OneD, Array<OneD, int>>>
        m_traceCoeffToLeftRightExpCoeffMap;
    /// The sign of every coeff from current trace to the left & right adjacent
    /// expasion coeffs.
    Array<OneD, Array<OneD, Array<OneD, int>>>
        m_traceCoeffToLeftRightExpCoeffSign;
    /// The map of every phys from current trace to the left & right adjacent
    /// expasion phys. This map is only used when no interpolation is needed in
    /// getting GetFwdBwdTracePhys. If interpolation is needed, it should be
    /// determined as the InnerProduct of m_locTraceToFieldMap matrix and
    /// interpolation matrix.
    Array<OneD, Array<OneD, Array<OneD, int>>> m_tracePhysToLeftRightExpPhysMap;
 
    // store the number of neighbor elements for each element
    Array<OneD, int> m_ElemNeighbsNumb;
    // store the id of neighbor elements for each element
    Array<OneD, Array<OneD, int>> m_ElemNeighbsId;
 
    Array<OneD, Array<OneD, int>> m_LocTracephysToTraceIDMap;
 
    void FindElmtNeighbors(const ExpList &locExp,
                           const ExpListSharedPtr &trace);
 
    //------- Below is new members to support cell-based trace operations ------
    /// start entry of each local trace in m_locTraceToFieldMap, referenced by
    /// element and trace id
    Array<OneD, Array<OneD, int>> m_locTracePtsEntry;
    /// start entry of each global trace in m_locInterpTraceToTraceMap,
    /// referenced by element and trace id
    Array<OneD, Array<OneD, int>> m_interpTracePtsEntry;
    /// subscript of m_interpTrace, referenced by element and trace id
    Array<OneD, Array<OneD, int>> m_interpTraceIndex;
    /// start entry of each global trace in m_traceCoeffsToElmtMap, referenced
    /// by element and trace id
    Array<OneD, Array<OneD, int>> m_traceCoeffsEntry;
 
    // stores start expid of each collection
    Array<OneD, int> m_collExpOffset;
 
    // From loctrace id to element trace phys offset
    Array<OneD, int> m_locToTracePhysOffset;
 
    // A series of structs, each is for the mapping of a collection
    Array<OneD, TraceFieldMapEssential> m_traceFieldMap;
    Array<OneD, TraceInterpEssential> m_traceInterp;
};
 
typedef std::shared_ptr<LocTraceToTraceMap> LocTraceToTraceMapSharedPtr;
 
static LocTraceToTraceMapSharedPtr NullLocTraceToTraceMapSharedPtr;
 
} // namespace MultiRegions
} // namespace Nektar
 
#endif