Expansion1D.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: Expansion1D.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: File for Expansion1D routines
//
///////////////////////////////////////////////////////////////////////////////
 
#include <LocalRegions/Expansion1D.h>
#include <LocalRegions/Expansion2D.h>
 
using namespace std;
 
namespace Nektar::LocalRegions
{
 
DNekMatSharedPtr Expansion1D::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
{
    DNekMatSharedPtr returnval;
 
    switch (mkey.GetMatrixType())
    {
        case StdRegions::eHybridDGHelmholtz:
        {
            ASSERTL1(IsBoundaryInteriorExpansion(),
                     "HybridDGHelmholtz matrix not set up "
                     "for non boundary-interior expansions");
            int i;
            NekDouble lambdaval =
                mkey.GetConstFactor(StdRegions::eFactorLambda);
            NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
            int ncoeffs   = GetNcoeffs();
 
            int coordim = GetCoordim();
 
            DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
            StdRegions::MatrixType DerivType[3] = {StdRegions::eWeakDeriv0,
                                                   StdRegions::eWeakDeriv1,
                                                   StdRegions::eWeakDeriv2};
            DNekMat LocMat(ncoeffs, ncoeffs);
 
            returnval =
                MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs, ncoeffs);
            DNekMat &Mat = *returnval;
 
            Vmath::Zero(ncoeffs * ncoeffs, Mat.GetPtr(), 1);
 
            for (i = 0; i < coordim; ++i)
            {
                DNekScalMat &Dmat = *GetLocMatrix(DerivType[i]);
 
                Mat = Mat + Dmat * invMass * Transpose(Dmat);
            }
 
            // Add end Mass Matrix Contribution
            DNekScalMat &Mass = *GetLocMatrix(StdRegions::eMass);
            Mat               = Mat + lambdaval * Mass;
 
            Array<OneD, unsigned int> bmap;
            GetBoundaryMap(bmap);
 
            // Add tau*F_e using elemental mass matrices
            for (i = 0; i < 2; ++i)
            {
                Mat(bmap[i], bmap[i]) = Mat(bmap[i], bmap[i]) + tau;
            }
        }
        break;
        case StdRegions::eHybridDGLamToU:
        {
            int j, k;
            int nbndry  = NumDGBndryCoeffs();
            int ncoeffs = GetNcoeffs();
            StdRegions::ConstFactorMap factors;
            factors[StdRegions::eFactorLambda] =
                mkey.GetConstFactor(StdRegions::eFactorLambda);
            factors[StdRegions::eFactorTau] =
                mkey.GetConstFactor(StdRegions::eFactorTau);
 
            Array<OneD, NekDouble> lambda(nbndry);
            DNekVec Lambda(nbndry, lambda, eWrapper);
            Array<OneD, NekDouble> ulam(ncoeffs);
            DNekVec Ulam(ncoeffs, ulam, eWrapper);
            Array<OneD, NekDouble> f(ncoeffs);
            DNekVec F(ncoeffs, f, eWrapper);
 
            // declare matrix space
            returnval =
                MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs, nbndry);
            DNekMat &Umat = *returnval;
 
            // Helmholtz matrix
            DNekScalMat &invHmat =
                *GetLocMatrix(StdRegions::eInvHybridDGHelmholtz, factors);
 
            // for each degree of freedom of the lambda space
            // calculate Umat entry
            // Generate Lambda to U_lambda matrix
            for (j = 0; j < nbndry; ++j)
            {
                Vmath::Zero(nbndry, &lambda[0], 1);
                Vmath::Zero(ncoeffs, &f[0], 1);
                lambda[j] = 1.0;
 
                AddHDGHelmholtzTraceTerms(factors[StdRegions::eFactorTau],
                                          lambda, f);
 
                Ulam = invHmat * F; // generate Ulam from lambda
 
                // fill column of matrix
                for (k = 0; k < ncoeffs; ++k)
                {
                    Umat(k, j) = Ulam[k];
                }
            }
        }
        break;
        case StdRegions::eHybridDGLamToQ0:
        case StdRegions::eHybridDGLamToQ1:
        case StdRegions::eHybridDGLamToQ2:
        {
            int j       = 0;
            int k       = 0;
            int dir     = 0;
            int nbndry  = NumDGBndryCoeffs();
            int ncoeffs = GetNcoeffs();
 
            Array<OneD, NekDouble> lambda(nbndry);
            DNekVec Lambda(nbndry, lambda, eWrapper);
 
            Array<OneD, NekDouble> ulam(ncoeffs);
            DNekVec Ulam(ncoeffs, ulam, eWrapper);
            Array<OneD, NekDouble> f(ncoeffs);
            DNekVec F(ncoeffs, f, eWrapper);
            StdRegions::ConstFactorMap factors;
            factors[StdRegions::eFactorLambda] =
                mkey.GetConstFactor(StdRegions::eFactorLambda);
            factors[StdRegions::eFactorTau] =
                mkey.GetConstFactor(StdRegions::eFactorTau);
 
            // declare matrix space
            returnval =
                MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs, nbndry);
            DNekMat &Qmat = *returnval;
 
            // Lambda to U matrix
            DNekScalMat &lamToU =
                *GetLocMatrix(StdRegions::eHybridDGLamToU, factors);
 
            // Inverse mass matrix
            DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
 
            // Weak Derivative matrix
            DNekScalMatSharedPtr Dmat;
            switch (mkey.GetMatrixType())
            {
                case StdRegions::eHybridDGLamToQ0:
                    dir  = 0;
                    Dmat = GetLocMatrix(StdRegions::eWeakDeriv0);
                    break;
                case StdRegions::eHybridDGLamToQ1:
                    dir  = 1;
                    Dmat = GetLocMatrix(StdRegions::eWeakDeriv1);
                    break;
                case StdRegions::eHybridDGLamToQ2:
                    dir  = 2;
                    Dmat = GetLocMatrix(StdRegions::eWeakDeriv2);
                    break;
                default:
                    ASSERTL0(false, "Direction not known");
                    break;
            }
 
            // for each degree of freedom of the lambda space
            // calculate Qmat entry
            // Generate Lambda to Q_lambda matrix
            for (j = 0; j < nbndry; ++j)
            {
                Vmath::Zero(nbndry, &lambda[0], 1);
                lambda[j] = 1.0;
 
                // for lambda[j] = 1 this is the solution to ulam
                for (k = 0; k < ncoeffs; ++k)
                {
                    Ulam[k] = lamToU(k, j);
                }
 
                // -D^T ulam
                Vmath::Neg(ncoeffs, &ulam[0], 1);
                F = Transpose(*Dmat) * Ulam;
 
                // + \tilde{G} \lambda
                AddNormTraceInt(dir, lambda, f);
 
                // multiply by inverse mass matrix
                Ulam = invMass * F;
 
                // fill column of matrix (Qmat is in column major format)
                Vmath::Vcopy(ncoeffs, &ulam[0], 1,
                             &(Qmat.GetPtr())[0] + j * ncoeffs, 1);
            }
        }
        break;
        case StdRegions::eHybridDGHelmBndLam:
        {
            int j;
            int nbndry = NumBndryCoeffs();
 
            StdRegions::ConstFactorMap factors;
            factors[StdRegions::eFactorLambda] =
                mkey.GetConstFactor(StdRegions::eFactorLambda);
            factors[StdRegions::eFactorTau] =
                mkey.GetConstFactor(StdRegions::eFactorTau);
 
            Array<OneD, unsigned int> bmap;
            Array<OneD, NekDouble> lam(2);
            GetBoundaryMap(bmap);
 
            // declare matrix space
            returnval =
                MemoryManager<DNekMat>::AllocateSharedPtr(nbndry, nbndry);
            DNekMat &BndMat = *returnval;
 
            // Matrix to map Lambda to U
            DNekScalMat &LamToU =
                *GetLocMatrix(StdRegions::eHybridDGLamToU, factors);
 
            // Matrix to map Lambda to Q
            DNekScalMat &LamToQ =
                *GetLocMatrix(StdRegions::eHybridDGLamToQ0, factors);
 
            lam[0] = 1.0;
            lam[1] = 0.0;
            for (j = 0; j < nbndry; ++j)
            {
                BndMat(0, j) =
                    -LamToQ(bmap[0], j) - factors[StdRegions::eFactorTau] *
                                              (LamToU(bmap[0], j) - lam[j]);
            }
 
            lam[0] = 0.0;
            lam[1] = 1.0;
            for (j = 0; j < nbndry; ++j)
            {
                BndMat(1, j) =
                    LamToQ(bmap[1], j) - factors[StdRegions::eFactorTau] *
                                             (LamToU(bmap[1], j) - lam[j]);
            }
        }
        break;
        default:
            ASSERTL0(false,
                     "This matrix type cannot be generated from this class");
            break;
    }
 
    return returnval;
}
 
void Expansion1D::v_PhysDeriv(const int dir,
                              const Array<OneD, const NekDouble> &inarray,
                              Array<OneD, NekDouble> &outarray)
{
    switch (dir)
    {
        case 0:
        {
            v_PhysDeriv(inarray, outarray, NullNekDouble1DArray,
                        NullNekDouble1DArray);
            break;
        }
 
        case 1:
        {
            v_PhysDeriv(inarray, NullNekDouble1DArray, outarray,
                        NullNekDouble1DArray);
            break;
        }
 
        case 2:
        {
            v_PhysDeriv(inarray, NullNekDouble1DArray, NullNekDouble1DArray,
                        outarray);
            break;
        }
 
        default:
        {
            ASSERTL1(false, "input dir is out of range");
        }
        break;
    }
}
 
/**
   \brief  Inner product of \a inarray over region with respect to
   expansion basis \a base and return in \a outarray
 
   Calculate \f$ I[p] = \int^{1}_{-1} \phi_p(\xi_1) u(\xi_1) d\xi_1
   = \sum_{i=0}^{nq-1} \phi_p(\xi_{1i}) u(\xi_{1i}) w_i \f$ where
   \f$ outarray[p] = I[p], inarray[i] = u(\xi_{1i}), base[p*nq+i] =
   \phi_p(\xi_{1i}) \f$.
 
   Inputs: \n
   - \a inarray: physical point array of function to be integrated
   \f$ u(\xi_1) \f$
 
   Output: \n
 
   - \a outarray: array of coefficients representing the inner
   product of function with ever  mode in the exapnsion
 
**/
void Expansion1D::v_IProductWRTBase(const Array<OneD, const NekDouble> &inarray,
                                    Array<OneD, NekDouble> &outarray)
{
    const Array<OneD, const NekDouble> &jac = m_geomFactors->GetJac();
    bool Deformed = (m_geomFactors->GetGtype() == SpatialDomains::eDeformed);
 
    if (v_IsCollocatedBasis())
    {
        int nqtot = GetTotPoints();
        if (Deformed)
        {
            Vmath::Vmul(nqtot, jac, 1, inarray, 1, outarray, 1);
        }
        else
        {
            Vmath::Smul(nqtot, jac[0], inarray, 1, outarray, 1);
        }
        v_MultiplyByStdQuadratureMetric(outarray, outarray);
    }
    else
    {
        v_IProductWRTBaseKernel(m_base[0]->GetBdata(), inarray, outarray, jac,
                                Deformed);
    }
}
 
/** \brief Evaluate the derivative \f$ d/d{\xi_1} \f$ at the
    physical quadrature points given by \a inarray and return in \a
    outarray.
 
    This is a wrapper around StdExpansion1D::Tensor_Deriv
 
    Input:\n
 
    - \a n: number of derivatives to be evaluated where \f$ n \leq  dim\f$
 
    - \a inarray: array of function evaluated at the quadrature points
 
    Output: \n
 
    - \a outarray: array of the derivatives \f$
    du/d_{\xi_1}|_{\xi_{1i}} d\xi_1/dx,
    du/d_{\xi_1}|_{\xi_{1i}} d\xi_1/dy,
    du/d_{\xi_1}|_{\xi_{1i}} d\xi_1/dz,
    \f$ depending on value of \a dim
*/
void Expansion1D::v_PhysDeriv(const Array<OneD, const NekDouble> &inarray,
                              Array<OneD, NekDouble> &out_d0,
                              Array<OneD, NekDouble> &out_d1,
                              Array<OneD, NekDouble> &out_d2)
{
    int nquad0                        = m_base[0]->GetNumPoints();
    Array<TwoD, const NekDouble> gmat = m_geomFactors->GetDerivFactors();
    Array<OneD, NekDouble> diff(nquad0);
 
    PhysTensorDeriv(inarray, diff);
    if (m_geomFactors->GetGtype() == SpatialDomains::eDeformed)
    {
        if (out_d0.size())
        {
            Vmath::Vmul(nquad0, &gmat[0][0], 1, &diff[0], 1, &out_d0[0], 1);
        }
 
        if (out_d1.size())
        {
            Vmath::Vmul(nquad0, &gmat[1][0], 1, &diff[0], 1, &out_d1[0], 1);
        }
 
        if (out_d2.size())
        {
            Vmath::Vmul(nquad0, &gmat[2][0], 1, &diff[0], 1, &out_d2[0], 1);
        }
    }
    else
    {
        if (out_d0.size())
        {
            Vmath::Smul(nquad0, gmat[0][0], diff, 1, out_d0, 1);
        }
 
        if (out_d1.size())
        {
            Vmath::Smul(nquad0, gmat[1][0], diff, 1, out_d1, 1);
        }
 
        if (out_d2.size())
        {
            Vmath::Smul(nquad0, gmat[2][0], diff, 1, out_d2, 1);
        }
    }
}
 
void Expansion1D::AddNormTraceInt([[maybe_unused]] const int dir,
                                  Array<OneD, const NekDouble> &inarray,
                                  Array<OneD, NekDouble> &outarray)
{
    int k;
    int nbndry                                = NumBndryCoeffs();
    int nquad                                 = GetNumPoints(0);
    const Array<OneD, const NekDouble> &Basis = GetBasis(0)->GetBdata();
    Array<OneD, unsigned int> vmap;
 
    GetBoundaryMap(vmap);
 
    // add G \lambda term (can assume G is diagonal since one
    // of the basis is zero at boundary otherwise)
    for (k = 0; k < nbndry; ++k)
    {
        outarray[vmap[k]] += (Basis[(vmap[k] + 1) * nquad - 1] *
                                  Basis[(vmap[k] + 1) * nquad - 1] -
                              Basis[vmap[k] * nquad] * Basis[vmap[k] * nquad]) *
                             inarray[vmap[k]];
    }
}
 
void Expansion1D::AddHDGHelmholtzTraceTerms(
    const NekDouble tau, const Array<OneD, const NekDouble> &inarray,
    Array<OneD, NekDouble> &outarray)
{
    int i, n;
    int nbndry  = NumBndryCoeffs();
    int nquad   = GetNumPoints(0);
    int ncoeffs = GetNcoeffs();
    int coordim = GetCoordim();
    Array<OneD, unsigned int> vmap;
 
    ASSERTL0(&inarray[0] != &outarray[0],
             "Input and output arrays use the same memory");
 
    const Array<OneD, const NekDouble> &Basis = GetBasis(0)->GetBdata();
    DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
 
    GetBoundaryMap(vmap);
 
    // Add F = \tau <phi_i,phi_j> (note phi_i is zero if phi_j is non-zero)
    for (i = 0; i < nbndry; ++i)
    {
        outarray[vmap[i]] += tau * Basis[(vmap[i] + 1) * nquad - 1] *
                             Basis[(vmap[i] + 1) * nquad - 1] *
                             inarray[vmap[i]];
        outarray[vmap[i]] += tau * Basis[vmap[i] * nquad] *
                             Basis[vmap[i] * nquad] * inarray[vmap[i]];
    }
 
    //===============================================================
    // Add -\sum_i D_i^T M^{-1} G_i + E_i M^{-1} G_i =
    //                         \sum_i D_i M^{-1} G_i term
 
    StdRegions::MatrixType DerivType[3] = {StdRegions::eWeakDeriv0,
                                           StdRegions::eWeakDeriv1,
                                           StdRegions::eWeakDeriv2};
    Array<OneD, NekDouble> tmpcoeff(ncoeffs, 0.0);
    DNekVec Coeffs(ncoeffs, outarray, eWrapper);
    DNekVec Tmpcoeff(ncoeffs, tmpcoeff, eWrapper);
 
    for (n = 0; n < coordim; ++n)
    {
        // evaluate M^{-1} G
        for (i = 0; i < ncoeffs; ++i)
        {
            // lower boundary (negative normal)
            tmpcoeff[i] -= invMass(i, vmap[0]) * Basis[vmap[0] * nquad] *
                           Basis[vmap[0] * nquad] * inarray[vmap[0]];
 
            // upper boundary (positive normal)
            tmpcoeff[i] += invMass(i, vmap[1]) *
                           Basis[(vmap[1] + 1) * nquad - 1] *
                           Basis[(vmap[1] + 1) * nquad - 1] * inarray[vmap[1]];
        }
 
        DNekScalMat &Dmat = *GetLocMatrix(DerivType[n]);
        Coeffs            = Coeffs + Dmat * Tmpcoeff;
    }
}
 
void Expansion1D::v_AddRobinMassMatrix(
    const int vert, const Array<OneD, const NekDouble> &primCoeffs,
    DNekMatSharedPtr &inoutmat)
{
    ASSERTL0(IsBoundaryInteriorExpansion(),
             "Robin boundary conditions are only implemented for "
             "boundary-interior expanisons");
    ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
             "Assuming that input matrix was square");
 
    // Get local Element mapping for vertex point
    int map = GetVertexMap(vert);
 
    // Now need to identify a map which takes the local edge
    // mass matrix to the matrix stored in inoutmat;
    // This can currently be deduced from the size of the matrix
    // - if inoutmat.m_rows() == v_NCoeffs() it is a full
    //   matrix system
    // - if inoutmat.m_rows() == v_NumBndCoeffs() it is a
    //  boundary CG system
 
    int rows = inoutmat->GetRows();
 
    if (rows == GetNcoeffs())
    {
        // no need to do anything
    }
    else if (rows == NumBndryCoeffs()) // same as NumDGBndryCoeffs()
    {
        int i;
        Array<OneD, unsigned int> bmap;
        GetBoundaryMap(bmap);
 
        for (i = 0; i < 2; ++i)
        {
            if (map == bmap[i])
            {
                map = i;
                break;
            }
        }
        ASSERTL1(i != 2, "Did not find number in map");
    }
 
    // assumes end points have unit magnitude
    (*inoutmat)(map, map) += primCoeffs[0];
}
 
/**
 * Given an edge and vector of element coefficients:
 * - maps those elemental coefficients corresponding to the trace into
 *   an vector.
 * - update the element coefficients
 * - multiplies the edge vector by the edge mass matrix
 * - maps the edge coefficients back onto the elemental coefficients
 */
void Expansion1D::v_AddRobinTraceContribution(
    const int vert, const Array<OneD, const NekDouble> &primCoeffs,
    const Array<OneD, NekDouble> &incoeffs, Array<OneD, NekDouble> &coeffs)
{
    ASSERTL1(IsBoundaryInteriorExpansion(),
             "Not set up for non boundary-interior expansions");
 
    int map = GetVertexMap(vert);
    coeffs[map] += primCoeffs[0] * incoeffs[map];
}
 
NekDouble Expansion1D::v_VectorFlux(
    const Array<OneD, Array<OneD, NekDouble>> &vec)
{
    const Array<OneD, const Array<OneD, NekDouble>> &normals =
        GetLeftAdjacentElementExp()->GetTraceNormal(
            GetLeftAdjacentElementTrace());
 
    int nq = m_base[0]->GetNumPoints();
    Array<OneD, NekDouble> Fn(nq);
    Vmath::Vmul(nq, &vec[0][0], 1, &normals[0][0], 1, &Fn[0], 1);
    Vmath::Vvtvp(nq, &vec[1][0], 1, &normals[1][0], 1, &Fn[0], 1, &Fn[0], 1);
 
    return Integral(Fn);
}
 
/** @brief: This method gets all of the factors which are
    required as part of the Gradient Jump Penalty
    stabilisation and involves the product of the normal and
    geometric factors along the element trace.
*/
void Expansion1D::v_NormalTraceDerivFactors(
    Array<OneD, Array<OneD, NekDouble>> &factors,
    [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &d0factors,
    [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &d1factors)
{
    int nquad                         = GetNumPoints(0);
    Array<TwoD, const NekDouble> gmat = m_geomFactors->GetDerivFactors();
 
    if (factors.size() <= 2)
    {
        factors    = Array<OneD, Array<OneD, NekDouble>>(2);
        factors[0] = Array<OneD, NekDouble>(1);
        factors[1] = Array<OneD, NekDouble>(1);
    }
 
    // Outwards normal
    const Array<OneD, const Array<OneD, NekDouble>> &normal_0 =
        GetTraceNormal(0);
    const Array<OneD, const Array<OneD, NekDouble>> &normal_1 =
        GetTraceNormal(1);
 
    if (m_geomFactors->GetGtype() == SpatialDomains::eDeformed)
    {
        factors[0][0] = gmat[0][nquad - 1] * normal_0[0][0];
        factors[1][0] = gmat[0][0] * normal_1[0][0];
 
        for (int n = 1; n < normal_0.size(); ++n)
        {
            factors[0][0] += gmat[n][0] * normal_0[n][0];
            factors[1][0] += gmat[n][nquad - 1] * normal_1[n][0];
        }
    }
    else
    {
        factors[0][0] = gmat[0][0] * normal_0[0][0];
        factors[1][0] = gmat[0][0] * normal_1[0][0];
 
        for (int n = 1; n < normal_0.size(); ++n)
        {
            factors[0][0] += gmat[n][0] * normal_0[n][0];
            factors[1][0] += gmat[n][0] * normal_1[n][0];
        }
    }
}
 
void Expansion1D::v_ReOrientTracePhysMap(
    [[maybe_unused]] const StdRegions::Orientation orient,
    Array<OneD, int> &idmap, [[maybe_unused]] const int nq0,
    [[maybe_unused]] const int nq1, [[maybe_unused]] bool Forwards)
{
    if (idmap.size() != 1)
    {
        idmap = Array<OneD, int>(1);
    }
 
    idmap[0] = 0;
}
 
void Expansion1D::v_TraceNormLen([[maybe_unused]] const int traceid,
                                 NekDouble &h, NekDouble &p)
{
    h = GetGeom()->GetVertex(1)->dist(*GetGeom()->GetVertex(0));
    p = m_ncoeffs - 1;
}
 
} // namespace Nektar::LocalRegions