PhysGalerkinProject.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: PhysGalerkinProject.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: Definition of Physical Space Galerkin Projection methods
//
///////////////////////////////////////////////////////////////////////////////
#include <LibUtilities/BasicUtils/Vmath.hpp>
#include <LibUtilities/BasicUtils/VmathArray.hpp>
#include <LibUtilities/Foundations/PhysGalerkinProject.h>
 
#include <LibUtilities/Foundations/Basis.h>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <LibUtilities/Foundations/Points.h>
#include <LibUtilities/LinearAlgebra/NekTypeDefs.hpp>
 
namespace Nektar
{
namespace LibUtilities
{
 
// Physical Space Interpolation methods
 
// 1D Interpolation
void PhysGalerkinProject1D(const BasisKey &fbasis0,
                           const Array<OneD, const NekDouble> &from,
                           const BasisKey &tbasis0, Array<OneD, NekDouble> &to)
{
    PhysGalerkinProject1D(fbasis0.GetPointsKey(), from, tbasis0.GetPointsKey(),
                          to);
}
 
void PhysGalerkinProject1D(const PointsKey &fpoints0,
                           const Array<OneD, const NekDouble> &from,
                           const PointsKey &tpoints0,
                           Array<OneD, NekDouble> &to)
{
    if (fpoints0 == tpoints0) // check to see if the same
    {
        Vmath::Vcopy(fpoints0.GetNumPoints(), from, 1, to, 1);
    }
    else // interpolate
    {
        DNekMatSharedPtr GP0;
 
        GP0 = PointsManager()[tpoints0]->GetGalerkinProjection(fpoints0);
 
        NekVector<NekDouble> in(fpoints0.GetNumPoints(), from, eWrapper);
        NekVector<NekDouble> out(tpoints0.GetNumPoints(), to, eWrapper);
 
        GP0->Transpose();
        out = (*GP0) * in;
    }
}
 
void PhysGalerkinProject1D(const BasisKey &fbasis0, const NekDouble *from,
                           const BasisKey &tbasis0, NekDouble *to)
{
    PhysGalerkinProject1D(fbasis0.GetPointsKey(), from, tbasis0.GetPointsKey(),
                          to);
}
 
void PhysGalerkinProject1D(const PointsKey &fpoints0, const NekDouble *from,
                           const PointsKey &tpoints0, NekDouble *to)
{
    if (fpoints0 == tpoints0) // check to see if the same
    {
        Vmath::Vcopy(fpoints0.GetNumPoints(), from, 1, to, 1);
    }
    else // interpolate
    {
 
        DNekMatSharedPtr GP0;
 
        GP0 = PointsManager()[tpoints0]->GetGalerkinProjection(fpoints0);
 
        Blas::Dgemv('T', tpoints0.GetNumPoints(), fpoints0.GetNumPoints(), 1.0,
                    GP0->GetPtr().get(), tpoints0.GetNumPoints(), from, 1, 0.0,
                    to, 1);
    }
}
 
// 2D Galerkin Projection
void PhysGalerkinProject2D(const BasisKey &fbasis0, const BasisKey &fbasis1,
                           const Array<OneD, const NekDouble> &from,
                           const BasisKey &tbasis0, const BasisKey &tbasis1,
                           Array<OneD, NekDouble> &to)
{
    PhysGalerkinProject2D(fbasis0.GetPointsKey(), fbasis1.GetPointsKey(),
                          from.data(), tbasis0.GetPointsKey(),
                          tbasis1.GetPointsKey(), to.data());
}
 
void PhysGalerkinProject2D(const PointsKey &fpoints0, const PointsKey &fpoints1,
                           const Array<OneD, const NekDouble> &from,
                           const PointsKey &tpoints0, const PointsKey &tpoints1,
                           Array<OneD, NekDouble> &to)
{
    PhysGalerkinProject2D(fpoints0, fpoints1, from.data(), tpoints0, tpoints1,
                          to.data());
}
 
void PhysGalerkinProject2D(const PointsKey &fpoints0, const PointsKey &fpoints1,
                           const NekDouble *from, const PointsKey &tpoints0,
                           const PointsKey &tpoints1, NekDouble *to)
{
    DNekMatSharedPtr GP0, GP1;
    Array<OneD, NekDouble> wsp(tpoints1.GetNumPoints() *
                               fpoints0.GetNumPoints()); // fnp0*tnp1
 
    int fnp0 = fpoints0.GetNumPoints();
    int fnp1 = fpoints1.GetNumPoints();
    int tnp0 = tpoints0.GetNumPoints();
    int tnp1 = tpoints1.GetNumPoints();
 
    if (fpoints1 == tpoints1)
    {
        Vmath::Vcopy(fnp0 * tnp1, from, 1, wsp.get(), 1);
    }
    else
    {
        GP1 = PointsManager()[tpoints1]->GetGalerkinProjection(fpoints1);
        Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0, from, fnp0,
                    GP1->GetPtr().get(), tnp1, 0.0, wsp.get(), fnp0);
    }
 
    if (fpoints0 == tpoints0)
    {
        Vmath::Vcopy(tnp0 * tnp1, wsp.get(), 1, to, 1);
    }
    else
    {
        GP0 = PointsManager()[tpoints0]->GetGalerkinProjection(fpoints0);
        Blas::Dgemm('N', 'N', tnp0, tnp1, fnp0, 1.0, GP0->GetPtr().get(), tnp0,
                    wsp.get(), fnp0, 0.0, to, tnp0);
    }
}
 
// 3D Galerkin Projection
void PhysGalerkinProject3D(const BasisKey &fbasis0, const BasisKey &fbasis1,
                           const BasisKey &fbasis2,
                           const Array<OneD, const NekDouble> &from,
                           const BasisKey &tbasis0, const BasisKey &tbasis1,
                           const BasisKey &tbasis2, Array<OneD, NekDouble> &to)
{
    PhysGalerkinProject3D(fbasis0.GetPointsKey(), fbasis1.GetPointsKey(),
                          fbasis2.GetPointsKey(), from.data(),
                          tbasis0.GetPointsKey(), tbasis1.GetPointsKey(),
                          tbasis2.GetPointsKey(), to.data());
}
 
void PhysGalerkinProject3D(const PointsKey &fpoints0, const PointsKey &fpoints1,
                           const PointsKey &fpoints2,
                           const Array<OneD, const NekDouble> &from,
                           const PointsKey &tpoints0, const PointsKey &tpoints1,
                           const PointsKey &tpoints2,
                           Array<OneD, NekDouble> &to)
{
    PhysGalerkinProject3D(fpoints0, fpoints1, fpoints2, from.data(), tpoints0,
                          tpoints1, tpoints2, to.data());
}
 
void PhysGalerkinProject3D(const PointsKey &fpoints0, const PointsKey &fpoints1,
                           const PointsKey &fpoints2, const NekDouble *from,
                           const PointsKey &tpoints0, const PointsKey &tpoints1,
                           const PointsKey &tpoints2, NekDouble *to)
{
    DNekMatSharedPtr GP0, GP1, GP2;
 
    int fnp0 = fpoints0.GetNumPoints();
    int fnp1 = fpoints1.GetNumPoints();
    int fnp2 = fpoints2.GetNumPoints();
    int tnp0 = tpoints0.GetNumPoints();
    int tnp1 = tpoints1.GetNumPoints();
    int tnp2 = tpoints2.GetNumPoints();
 
    Array<OneD, NekDouble> wsp1(fnp0 * tnp1 * tnp2);
    Array<OneD, NekDouble> wsp2(fnp0 * fnp1 * tnp2);
 
    GP2 = PointsManager()[tpoints2]->GetGalerkinProjection(fpoints2);
    Blas::Dgemm('N', 'T', fnp0 * fnp1, tnp2, fnp2, 1.0, from, fnp0 * fnp1,
                GP2->GetPtr().get(), tnp2, 0.0, wsp2.get(), fnp0 * fnp1);
 
    GP1 = PointsManager()[tpoints1]->GetGalerkinProjection(fpoints1);
    for (int i = 0; i < tnp2; i++)
    {
        Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0,
                    wsp2.get() + i * fnp0 * fnp1, fnp0, GP1->GetPtr().get(),
                    tnp1, 0.0, wsp1.get() + i * fnp0 * tnp1, fnp0);
    }
 
    GP0 = PointsManager()[tpoints0]->GetGalerkinProjection(fpoints0);
    Blas::Dgemm('N', 'N', tnp0, tnp1 * tnp2, fnp0, 1.0, GP0->GetPtr().get(),
                tnp0, wsp1.get(), fnp0, 0.0, to, tnp0);
}
 
} // namespace LibUtilities
} // namespace Nektar