StdSegExp.cpp 30.5 KB
Newer Older
Mike Kirby's avatar
Mike Kirby committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
///////////////////////////////////////////////////////////////////////////////
//
// File StdSegExp.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.
mike turner's avatar
mike turner committed
31
//
Mike Kirby's avatar
Mike Kirby committed
32 33 34 35
// Description: Routines within Standard Segment Expansions
//
///////////////////////////////////////////////////////////////////////////////

Daniele de Grazia's avatar
Daniele de Grazia committed
36 37
#include <StdRegions/StdSegExp.h>
#include <LibUtilities/Foundations/InterpCoeff.h>
mike turner's avatar
mike turner committed
38

Mike Kirby's avatar
Mike Kirby committed
39

40 41
using namespace std;

Mike Kirby's avatar
Mike Kirby committed
42 43
namespace Nektar
{
44
    namespace StdRegions
Mike Kirby's avatar
Mike Kirby committed
45 46
    {

47 48 49 50 51 52 53
        /** \brief Default constructor */

        StdSegExp::StdSegExp()
        {
        }


mike turner's avatar
mike turner committed
54
        /** \brief Constructor using BasisKey class for quadrature points and
55 56
         *  order definition
         *
mike turner's avatar
mike turner committed
57
         *  \param Ba BasisKey class definition containing order and quadrature
58 59 60
         *  points.
         */

61
        StdSegExp::StdSegExp(const LibUtilities::BasisKey &Ba):
Pavel Burovskiy's avatar
Pavel Burovskiy committed
62 63
        StdExpansion(Ba.GetNumModes(), 1, Ba),
        StdExpansion1D(Ba.GetNumModes(),Ba)
64
        {
Peter Vos's avatar
Peter Vos committed
65
        }
66 67 68 69


        /** \brief Copy Constructor */

Peter Vos's avatar
Peter Vos committed
70
        StdSegExp::StdSegExp(const StdSegExp &T):
Pavel Burovskiy's avatar
Pavel Burovskiy committed
71 72
                StdExpansion(T),
                StdExpansion1D(T)
Peter Vos's avatar
Peter Vos committed
73 74
        {
        }
75 76


Peter Vos's avatar
Peter Vos committed
77
        StdSegExp::~StdSegExp()
78
        {
Peter Vos's avatar
Peter Vos committed
79 80
        }

Pavel Burovskiy's avatar
Pavel Burovskiy committed
81 82 83
        /** \brief Return Shape of region, using  ShapeType enum list.
         *  i.e. Segment
         */
84
        LibUtilities::ShapeType StdSegExp::v_DetShapeType() const
Pavel Burovskiy's avatar
Pavel Burovskiy committed
85
        {
86
            return LibUtilities::eSegment;
Pavel Burovskiy's avatar
Pavel Burovskiy committed
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
        }

        bool StdSegExp::v_IsBoundaryInteriorExpansion()
        {

            bool returnval = false;

            if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
            {
                returnval = true;
            }

            if(m_base[0]->GetBasisType() == LibUtilities::eGLL_Lagrange)
            {
                returnval = true;
            }

            return returnval;
        }



109 110

        //---------------------------------------------------------------------
Peter Vos's avatar
Peter Vos committed
111
        // Integration Methods
112 113
        //---------------------------------------------------------------------

mike turner's avatar
mike turner committed
114
        /** \brief Integrate the physical point list \a inarray over region
115 116
         *  and return the value
         *
mike turner's avatar
mike turner committed
117 118
         *  \param inarray definition of function to be integrated evauluated at
         *  quadrature point of expansion.
119 120 121
         *  \return returns \f$\int^1_{-1} u(\xi_1)d \xi_1 \f$ where \f$inarray[i]
         *  = u(\xi_{1i}) \f$
         */
Pavel Burovskiy's avatar
Pavel Burovskiy committed
122
        NekDouble StdSegExp::v_Integral(const Array<OneD, const NekDouble>& inarray )
Peter Vos's avatar
Peter Vos committed
123 124 125
        {
            NekDouble Int = 0.0;
            int    nquad0 = m_base[0]->GetNumPoints();
Peter Vos's avatar
Peter Vos committed
126
            Array<OneD, NekDouble> tmp(nquad0);
Peter Vos's avatar
Peter Vos committed
127 128
            Array<OneD, const NekDouble> z  = m_base[0]->GetZ();
            Array<OneD, const NekDouble> w0 = m_base[0]->GetW();
129

mike turner's avatar
mike turner committed
130
            // multiply by integration constants
Peter Vos's avatar
Peter Vos committed
131
            Vmath::Vmul(nquad0, inarray, 1, w0, 1, tmp, 1);
132

Peter Vos's avatar
Peter Vos committed
133
            Int = Vmath::Vsum(nquad0, tmp, 1);
134

Peter Vos's avatar
Peter Vos committed
135 136
            return Int;
        }
Peter Vos's avatar
Peter Vos committed
137

138 139 140 141 142 143 144 145



        //---------------------------------------------------------------------
        // Differentiation Methods
        //---------------------------------------------------------------------


mike turner's avatar
mike turner committed
146
        /** \brief Evaluate the derivative \f$ d/d{\xi_1} \f$ at the physical
147 148 149 150 151
         *  quadrature points given by \a inarray and return in \a outarray.
         *
         *  This is a wrapper around StdExpansion1D::Tensor_Deriv
         *  \param inarray array of a function evaluated at the quadrature points
         *  \param  outarray the resulting array of the derivative \f$
mike turner's avatar
mike turner committed
152
         *  du/d_{\xi_1}|_{\xi_{1i}} \f$ will be stored in the array \a outarra
153
         */
Pavel Burovskiy's avatar
Pavel Burovskiy committed
154

155
        void StdSegExp::v_PhysDeriv(const Array<OneD, const NekDouble>& inarray,
Pavel Burovskiy's avatar
Pavel Burovskiy committed
156 157 158
                Array<OneD, NekDouble> &out_d0,
                Array<OneD, NekDouble> &out_d1,
                Array<OneD, NekDouble> &out_d2)
159
        {
160
            PhysTensorDeriv(inarray,out_d0);
161
        }
162 163 164


        void StdSegExp::v_PhysDeriv(const int dir,
Pavel Burovskiy's avatar
Pavel Burovskiy committed
165 166
                const Array<OneD, const NekDouble>& inarray,
                Array<OneD, NekDouble> &outarray)
Peter Vos's avatar
Peter Vos committed
167
        {
168
            ASSERTL1(dir==0,"input dir is out of range");
Pavel Burovskiy's avatar
Pavel Burovskiy committed
169
            PhysTensorDeriv(inarray,outarray);
170
            // PhysDeriv(inarray, outarray);
Peter Vos's avatar
Peter Vos committed
171
        }
172

173
        void StdSegExp::v_StdPhysDeriv(
Pavel Burovskiy's avatar
Pavel Burovskiy committed
174 175 176 177
                const Array<OneD, const NekDouble>& inarray,
                Array<OneD, NekDouble> &out_d0,
                Array<OneD, NekDouble> &out_d1,
                Array<OneD, NekDouble> &out_d2)
178
        {
Pavel Burovskiy's avatar
Pavel Burovskiy committed
179
            PhysTensorDeriv(inarray,out_d0);
180
            // PhysDeriv(inarray, out_d0);
181
        }
182

183
        void StdSegExp::v_StdPhysDeriv(
Pavel Burovskiy's avatar
Pavel Burovskiy committed
184 185 186
                const int dir,
                const Array<OneD, const NekDouble>& inarray,
                Array<OneD, NekDouble> &outarray)
187
        {
188
            ASSERTL1(dir==0,"input dir is out of range");
Pavel Burovskiy's avatar
Pavel Burovskiy committed
189
            PhysTensorDeriv(inarray,outarray);
190
            // PhysDeriv(inarray, outarray);
191 192 193 194
        }



195 196 197 198 199 200 201 202 203
        //---------------------------------------------------------------------
        // Transforms
        //---------------------------------------------------------------------

        /** \brief Backward transform from coefficient space given
         *  in \a inarray and evaluate at the physical quadrature
         *  points \a outarray
         *
         *  Operation can be evaluated as \f$ u(\xi_{1i}) =
mike turner's avatar
mike turner committed
204
         *  \sum_{p=0}^{order-1} \hat{u}_p \phi_p(\xi_{1i}) \f$ or equivalently
205 206 207 208 209 210 211
         *  \f$ {\bf u} = {\bf B}^T {\bf \hat{u}} \f$ where
         *  \f${\bf B}[i][j] = \phi_i(\xi_{1j}), \mbox{\_coeffs}[p] = {\bf
         *  \hat{u}}[p] \f$
         *
         *  The function takes the coefficient array \a inarray as
         *  input for the transformation
         *
mike turner's avatar
mike turner committed
212
         *  \param inarray: the coeffficients of the expansion
213
         *
mike turner's avatar
mike turner committed
214
         *  \param outarray: the resulting array of the values of the function at
215 216 217 218
         *  the physical quadrature points will be stored in the array \a outarray
         */

        void StdSegExp::v_BwdTrans(
Pavel Burovskiy's avatar
Pavel Burovskiy committed
219 220
                const Array<OneD, const NekDouble>& inarray,
                Array<OneD, NekDouble> &outarray)
221
        {
Peter Vos's avatar
Peter Vos committed
222
            int  nquad = m_base[0]->GetNumPoints();
223

Peter Vos's avatar
Peter Vos committed
224 225
            if(m_base[0]->Collocation())
            {
Peter Vos's avatar
Peter Vos committed
226
                Vmath::Vcopy(nquad, inarray, 1, outarray, 1);
Peter Vos's avatar
Peter Vos committed
227 228 229
            }
            else
            {
230

231 232 233
            Blas::Dgemv('N',nquad,m_base[0]->GetNumModes(),1.0,
                        (m_base[0]->GetBdata()).get(),
                        nquad,&inarray[0],1,0.0,&outarray[0],1);
234

Peter Vos's avatar
Peter Vos committed
235
            }
236
        }
237

238
        void StdSegExp::v_ReduceOrderCoeffs(
Daniele de Grazia's avatar
Daniele de Grazia committed
239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265
            int                                 numMin,
            const Array<OneD, const NekDouble> &inarray,
                  Array<OneD,       NekDouble> &outarray)
        {
            int n_coeffs = inarray.num_elements();

            Array<OneD, NekDouble> coeff(n_coeffs);
            Array<OneD, NekDouble> coeff_tmp(n_coeffs,0.0);
            Array<OneD, NekDouble> tmp;
            Array<OneD, NekDouble> tmp2;

            int       nmodes0 = m_base[0]->GetNumModes();

            Vmath::Vcopy(n_coeffs,inarray,1,coeff_tmp,1);

            const LibUtilities::PointsKey Pkey0(
                nmodes0, LibUtilities::eGaussLobattoLegendre);

            LibUtilities::BasisKey b0(m_base[0]->GetBasisType(),nmodes0,Pkey0);

            LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A,nmodes0,Pkey0);

            LibUtilities::InterpCoeff1D(
                b0, coeff_tmp, bortho0, coeff);

            Vmath::Zero(n_coeffs,coeff_tmp,1);

266 267 268
            Vmath::Vcopy(numMin,
                         tmp  = coeff,1,
                         tmp2 = coeff_tmp,1);
Daniele de Grazia's avatar
Daniele de Grazia committed
269

270
            LibUtilities::InterpCoeff1D(
Daniele de Grazia's avatar
Daniele de Grazia committed
271 272
                bortho0, coeff_tmp, b0, outarray);
        }
273

274 275 276 277
        /**
         * \brief Forward transform from physical quadrature space stored in \a
         * inarray and evaluate the expansion coefficients and store in \a
         * outarray
278
         *
279 280 281 282 283 284
         * Perform a forward transform using a Galerkin projection by taking the
         * inner product of the physical points and multiplying by the inverse
         * of the mass matrix using the Solve method of the standard matrix
         * container holding the local mass matrix, i.e. \f$ {\bf \hat{u}} =
         * {\bf M}^{-1} {\bf I} \f$ where \f$ {\bf I}[p] = \int^1_{-1}
         * \phi_p(\xi_1) u(\xi_1) d\xi_1 \f$
285
         *
286 287
         * This function stores the expansion coefficients calculated by the
         * transformation in the coefficient space array \a outarray
288
         *
289 290 291
         * \param inarray: array of physical quadrature points to be transformed
         * \param outarray: the coeffficients of the expansion
         */
292
        void StdSegExp::v_FwdTrans(const Array<OneD, const NekDouble>& inarray,
Pavel Burovskiy's avatar
Pavel Burovskiy committed
293
                Array<OneD, NekDouble> &outarray)
Peter Vos's avatar
Peter Vos committed
294 295 296
        {
            if(m_base[0]->Collocation())
            {
Peter Vos's avatar
Peter Vos committed
297
                Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
Peter Vos's avatar
Peter Vos committed
298 299 300
            }
            else
            {
301 302
                v_IProductWRTBase(inarray,outarray);

Peter Vos's avatar
Peter Vos committed
303
                // get Mass matrix inverse
304
                StdMatrixKey      masskey(eInvMass,v_DetShapeType(),*this);
305
                DNekMatSharedPtr  matsys = GetStdMatrix(masskey);
306

307
                NekVector<NekDouble> in(m_ncoeffs,outarray,eCopy);
Peter Vos's avatar
Peter Vos committed
308
                NekVector<NekDouble> out(m_ncoeffs,outarray,eWrapper);
309

Peter Vos's avatar
Peter Vos committed
310 311
                out = (*matsys)*in;
            }
312
        }
313

314
        void StdSegExp::v_FwdTrans_BndConstrained(
Pavel Burovskiy's avatar
Pavel Burovskiy committed
315 316
                const Array<OneD, const NekDouble>& inarray,
                Array<OneD, NekDouble> &outarray)
317 318 319 320 321 322 323 324 325
        {
            if(m_base[0]->Collocation())
            {
                Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
            }
            else
            {
                int nInteriorDofs = m_ncoeffs-2;
                int offset;
326

327 328 329 330 331 332 333
                switch(m_base[0]->GetBasisType())
                {
                case LibUtilities::eGLL_Lagrange:
                    {
                        offset = 1;
                    }
                    break;
Daniele de Grazia's avatar
Daniele de Grazia committed
334 335
                case LibUtilities::eGauss_Lagrange:
                {
336
                    nInteriorDofs = m_ncoeffs;
Daniele de Grazia's avatar
Daniele de Grazia committed
337 338 339
                    offset = 0;
                }
                break;
340
                case LibUtilities::eModified_A:
341
                case LibUtilities::eModified_B:
342 343 344 345 346 347
                    {
                        offset = 2;
                    }
                    break;
                default:
                    ASSERTL0(false,"This type of FwdTrans is not defined for this expansion type");
348
                }
349 350

                fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 );
mike turner's avatar
mike turner committed
351

Daniele de Grazia's avatar
Daniele de Grazia committed
352
                if(m_base[0]->GetBasisType() != LibUtilities::eGauss_Lagrange)
Peter Vos's avatar
Peter Vos committed
353
                {
Daniele de Grazia's avatar
Daniele de Grazia committed
354 355 356 357 358 359 360 361 362
                    outarray[GetVertexMap(0)] = inarray[0];
                    outarray[GetVertexMap(1)] = inarray[m_base[0]->GetNumPoints()-1];

                    if(m_ncoeffs>2)
                    {
                        // ideally, we would like to have tmp0 to be replaced by
                        // outarray (currently MassMatrixOp does not allow aliasing)
                        Array<OneD, NekDouble> tmp0(m_ncoeffs);
                        Array<OneD, NekDouble> tmp1(m_ncoeffs);
363

Daniele de Grazia's avatar
Daniele de Grazia committed
364 365 366
                        StdMatrixKey      masskey(eMass,v_DetShapeType(),*this);
                        MassMatrixOp(outarray,tmp0,masskey);
                        v_IProductWRTBase(inarray,tmp1);
367

Daniele de Grazia's avatar
Daniele de Grazia committed
368
                        Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
369

Daniele de Grazia's avatar
Daniele de Grazia committed
370 371 372
                        // get Mass matrix inverse (only of interior DOF)
                        DNekMatSharedPtr  matsys =
                        (m_stdStaticCondMatrixManager[masskey])-> GetBlock(1,1);
373

Daniele de Grazia's avatar
Daniele de Grazia committed
374 375 376 377 378 379 380 381
                        Blas::Dgemv('N',nInteriorDofs,nInteriorDofs,1.0,
                                &(matsys->GetPtr())[0],nInteriorDofs,tmp1.get()
                                +offset,1,0.0,outarray.get()+offset,1);
                    }
                }
                else
                {
                    StdSegExp::v_FwdTrans(inarray, outarray);
Peter Vos's avatar
Peter Vos committed
382
                }
383
            }
Pavel Burovskiy's avatar
Pavel Burovskiy committed
384

385
        }
Peter Vos's avatar
Peter Vos committed
386

387 388

        void StdSegExp::v_BwdTrans_SumFac(const Array<OneD, const NekDouble>& inarray,
Pavel Burovskiy's avatar
Pavel Burovskiy committed
389
                Array<OneD, NekDouble> &outarray)
390 391 392 393 394 395 396 397 398 399 400 401
        {
            v_BwdTrans(inarray, outarray);
        }



        //---------------------------------------------------------------------
        // Inner product functions
        //---------------------------------------------------------------------



mike turner's avatar
mike turner committed
402 403
        /** \brief  Inner product of \a inarray over region with respect to
         *  expansion basis \a base and return in \a outarray
404 405 406 407 408 409
         *
         *  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$.
         *
mike turner's avatar
mike turner committed
410
         *  \param  base an array defining the local basis for the inner product
411 412 413
         *  usually passed from Basis->GetBdata() or Basis->GetDbdata()
         *  \param inarray: physical point array of function to be integrated
         *  \f$ u(\xi_1) \f$
mike turner's avatar
mike turner committed
414 415 416
         *  \param coll_check flag to identify when a Basis->Collocation() call
         *  should be performed to see if this is a GLL_Lagrange basis with a
         *  collocation property. (should be set to 0 if taking the inner
417
         *  product with respect to the derivative of basis)
mike turner's avatar
mike turner committed
418
         *  \param outarray  the values of the inner product with respect to
419 420 421
         *  each basis over region will be stored in the array \a outarray as
         *  output of the function
         */
Pavel Burovskiy's avatar
Pavel Burovskiy committed
422

423 424 425 426 427
        void StdSegExp::v_IProductWRTBase(
            const Array<OneD, const NekDouble> &base,
            const Array<OneD, const NekDouble> &inarray,
                  Array<OneD,       NekDouble> &outarray,
            int coll_check)
428 429 430 431 432 433
        {
            int    nquad = m_base[0]->GetNumPoints();
            Array<OneD, NekDouble> tmp(nquad);
            Array<OneD, const NekDouble> w =  m_base[0]->GetW();

            Vmath::Vmul(nquad, inarray, 1, w, 1, tmp, 1);
mike turner's avatar
mike turner committed
434

435
            /* Comment below was a bug for collocated basis
436 437 438 439 440 441 442 443
            if(coll_check&&m_base[0]->Collocation())
            {
                Vmath::Vcopy(nquad, tmp, 1, outarray, 1);
            }
            else
            {
                Blas::Dgemv('T',nquad,m_ncoeffs,1.0,base.get(),nquad,
                            &tmp[0],1,0.0,outarray.get(),1);
444
            }*/
mike turner's avatar
mike turner committed
445

446 447 448
            // Correct implementation
            Blas::Dgemv('T',nquad,m_ncoeffs,1.0,base.get(),nquad,
                        &tmp[0],1,0.0,outarray.get(),1);
449 450 451 452 453 454 455 456
        }

        /** \brief Inner product of \a inarray over region with respect to the
         *  expansion basis (this)->m_base[0] and return in \a outarray
         *
         *  Wrapper call to StdSegExp::IProductWRTBase
         *  \param inarray array of function values evaluated at the physical
         *  collocation points
mike turner's avatar
mike turner committed
457
         *  \param outarray  the values of the inner product with respect to
458 459 460
         *  each basis over region will be stored in the array \a outarray as
         *  output of the function
         */
mike turner's avatar
mike turner committed
461
        void StdSegExp::v_IProductWRTBase(const Array<OneD, const NekDouble>& inarray,
Pavel Burovskiy's avatar
Pavel Burovskiy committed
462
                Array<OneD, NekDouble> &outarray)
463 464 465 466
        {
            v_IProductWRTBase(m_base[0]->GetBdata(),inarray,outarray,1);
        }

Pavel Burovskiy's avatar
Pavel Burovskiy committed
467 468 469 470 471 472 473 474
        void StdSegExp::v_IProductWRTDerivBase(
                const int dir,
                const Array<OneD, const NekDouble>& inarray,
                Array<OneD, NekDouble> & outarray)
      {
            ASSERTL1(dir >= 0 && dir < 1,"input dir is out of range");
            v_IProductWRTBase(m_base[0]->GetDbdata(),inarray,outarray,1);
        }
475 476

        void StdSegExp::v_IProductWRTBase_SumFac(
Pavel Burovskiy's avatar
Pavel Burovskiy committed
477
                const Array<OneD, const NekDouble>& inarray,
478 479
                Array<OneD, NekDouble> &outarray,
                bool multiplybyweights)
480
        {
481 482 483 484
            int    nquad = m_base[0]->GetNumPoints();
            Array<OneD, NekDouble> tmp(nquad);
            Array<OneD, const NekDouble> w =  m_base[0]->GetW();
            Array<OneD, const NekDouble> base =  m_base[0]->GetBdata();
485

486 487 488 489 490 491 492 493 494 495 496 497
            if(multiplybyweights)
            {
                Vmath::Vmul(nquad, inarray, 1, w, 1, tmp, 1);

                Blas::Dgemv('T',nquad,m_ncoeffs,1.0,base.get(),nquad,
                                &tmp[0],1,0.0,outarray.get(),1);
            }
            else
            {
                Blas::Dgemv('T',nquad,m_ncoeffs,1.0,base.get(),nquad,
                            &inarray[0],1,0.0,outarray.get(),1);
            }
498 499 500 501 502 503
        }





Pavel Burovskiy's avatar
Pavel Burovskiy committed
504 505 506
        //----------------------------
        // Evaluation
        //----------------------------
507 508 509 510 511 512 513 514 515 516

        void StdSegExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
        {
            int   nquad = m_base[0]->GetNumPoints();
            const NekDouble * base  = m_base[0]->GetBdata().get();

            ASSERTL2(mode <= m_ncoeffs,
             "calling argument mode is larger than total expansion order");

            Vmath::Vcopy(nquad,(NekDouble *)base+mode*nquad,1, &outarray[0],1);
Pavel Burovskiy's avatar
Pavel Burovskiy committed
517
        }
518

519 520

        NekDouble StdSegExp::v_PhysEvaluate(
521 522
                const Array<OneD, const NekDouble>& coords,
                const Array<OneD, const NekDouble>& physvals)
523
        {
524
            return  StdExpansion1D::v_PhysEvaluate(coords, physvals);
525 526
        }

527 528 529 530
        void StdSegExp::v_LaplacianMatrixOp(
            const Array<OneD, const NekDouble> &inarray,
                  Array<OneD,       NekDouble> &outarray,
            const StdMatrixKey                 &mkey)
531
        {
532 533 534 535 536 537 538 539 540 541 542 543 544
            int    nquad = m_base[0]->GetNumPoints();

            Array<OneD,NekDouble> physValues(nquad);
            Array<OneD,NekDouble> dPhysValuesdx(nquad);

            v_BwdTrans(inarray,physValues);

            // Laplacian matrix operation
            v_PhysDeriv(physValues,dPhysValuesdx);
            v_IProductWRTBase(m_base[0]->GetDbdata(),dPhysValuesdx,outarray,1);
        }


Pavel Burovskiy's avatar
Pavel Burovskiy committed
545
        void StdSegExp::v_HelmholtzMatrixOp(
546 547 548
            const Array<OneD, const NekDouble> &inarray,
                  Array<OneD,       NekDouble> &outarray,
            const StdMatrixKey                 &mkey)
549 550 551 552 553 554 555 556 557 558 559 560 561 562 563
        {
            int    nquad = m_base[0]->GetNumPoints();

            Array<OneD,NekDouble> physValues(nquad);
            Array<OneD,NekDouble> dPhysValuesdx(nquad);
            Array<OneD,NekDouble> wsp(m_ncoeffs);

            v_BwdTrans(inarray,physValues);

            // mass matrix operation
            v_IProductWRTBase((m_base[0]->GetBdata()),physValues,wsp,1);

            // Laplacian matrix operation
            v_PhysDeriv(physValues,dPhysValuesdx);
            v_IProductWRTBase(m_base[0]->GetDbdata(),dPhysValuesdx,outarray,1);
564
            Blas::Daxpy(m_ncoeffs, mkey.GetConstFactor(eFactorLambda), wsp.get(), 1, outarray.get(), 1);
565
        }
566

Douglas Serson's avatar
Douglas Serson committed
567 568 569
        void StdSegExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
                                            const StdMatrixKey &mkey)
        {
570
            // Generate an orthogonal expansion
Douglas Serson's avatar
Douglas Serson committed
571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607
            int nq = m_base[0]->GetNumPoints();
            int nmodes = m_base[0]->GetNumModes();
            // Declare orthogonal basis.
            LibUtilities::PointsKey pKey(nq,m_base[0]->GetPointsType());

            LibUtilities::BasisKey B(LibUtilities::eOrtho_A, nmodes, pKey);
            StdSegExp OrthoExp(B);

            //SVV parameters loaded from the .xml case file
            NekDouble  SvvDiffCoeff  = mkey.GetConstFactor(eFactorSVVDiffCoeff);
            int cutoff = (int) (mkey.GetConstFactor(eFactorSVVCutoffRatio))*nmodes;

            Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());

            // project onto modal  space.
            OrthoExp.FwdTrans(array,orthocoeffs);

            //
            for(int j = 0; j < nmodes; ++j)
            {
                if(j >= cutoff)//to filter out only the "high-modes"
                {
                     orthocoeffs[j] *=
                         (SvvDiffCoeff*exp(-(j-nmodes)*(j-nmodes)/
                                            ((NekDouble)((j-cutoff+1)*
                                                 (j-cutoff+1)))));
                 }
                else
                {
                     orthocoeffs[j] *= 0.0;
                }
            }

            // backward transform to physical space
            OrthoExp.BwdTrans(orthocoeffs,array);
        }

608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640
        void StdSegExp::v_ExponentialFilter(
                                          Array<OneD, NekDouble> &array,
                                    const NekDouble        alpha,
                                    const NekDouble        exponent,
                                    const NekDouble        cutoff)
        {
            // Generate an orthogonal expansion
            int nq     = m_base[0]->GetNumPoints();
            int nmodes = m_base[0]->GetNumModes();
            int P  = nmodes - 1;
            // Declare orthogonal basis.
            LibUtilities::PointsKey pKey(nq,m_base[0]->GetPointsType());

            LibUtilities::BasisKey B(LibUtilities::eOrtho_A, nmodes, pKey);
            StdSegExp OrthoExp(B);

            // Cutoff
            int Pcut = cutoff*P;

            // Project onto orthogonal space.
            Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
            OrthoExp.FwdTrans(array,orthocoeffs);

            //
            NekDouble fac;
            for(int j = 0; j < nmodes; ++j)
            {
                //to filter out only the "high-modes"
                if(j > Pcut)
                {
                    fac = (NekDouble) (j - Pcut) / ( (NekDouble) (P - Pcut) );
                    fac = pow(fac, exponent);
                    orthocoeffs[j] *= exp(-alpha*fac);
Douglas Serson's avatar
Douglas Serson committed
641
                }
642 643 644 645 646 647
            }

            // backward transform to physical space
            OrthoExp.BwdTrans(orthocoeffs,array);
        }

648 649 650 651
        //up to here
        void StdSegExp::v_MultiplyByStdQuadratureMetric(
            const Array<OneD, const NekDouble> &inarray,
                  Array<OneD,       NekDouble> &outarray)
mike turner's avatar
mike turner committed
652
        {
653
            int nquad0 = m_base[0]->GetNumPoints();
mike turner's avatar
mike turner committed
654

655 656 657 658 659
            const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();

            Vmath::Vmul(nquad0, inarray.get(),1,
                        w0.get(),1,outarray.get(),1);
        }
660 661

        void StdSegExp::v_GetCoords(
Pavel Burovskiy's avatar
Pavel Burovskiy committed
662 663 664
                Array<OneD, NekDouble> &coords_0,
                Array<OneD, NekDouble> &coords_1,
                Array<OneD, NekDouble> &coords_2)
665
        {
Pavel Burovskiy's avatar
Pavel Burovskiy committed
666 667
            Blas::Dcopy(GetNumPoints(0),(m_base[0]->GetZ()).get(),
                        1,&coords_0[0],1);
668
        }
669

Peter Vos's avatar
Peter Vos committed
670

671

672 673 674 675 676 677


        //---------------------------------------------------------------------
        // Helper functions
        //---------------------------------------------------------------------

Pavel Burovskiy's avatar
Pavel Burovskiy committed
678 679 680
        int StdSegExp::v_GetNverts() const
        {
            return 2;
mike turner's avatar
mike turner committed
681
        }
682

Pavel Burovskiy's avatar
Pavel Burovskiy committed
683
        int StdSegExp::v_NumBndryCoeffs() const
684
        {
Pavel Burovskiy's avatar
Pavel Burovskiy committed
685
            return 2;
mike turner's avatar
mike turner committed
686
        }
687

Pavel Burovskiy's avatar
Pavel Burovskiy committed
688 689 690
        int StdSegExp::v_NumDGBndryCoeffs() const
        {
            return 2;
mike turner's avatar
mike turner committed
691
        }
692

Pavel Burovskiy's avatar
Pavel Burovskiy committed
693 694 695 696 697 698
        int StdSegExp::v_CalcNumberOfCoefficients(
                const std::vector<unsigned int> &nummodes,
                int &modes_offset)
        {
            int nmodes = nummodes[modes_offset];
            modes_offset += 1;
699

Pavel Burovskiy's avatar
Pavel Burovskiy committed
700
            return nmodes;
701 702 703 704 705 706 707 708 709 710 711 712 713
        }

        //---------------------------------------------------------------------
        // Wrapper functions
        //---------------------------------------------------------------------

        DNekMatSharedPtr StdSegExp::v_GenMatrix(const StdMatrixKey &mkey)
        {
            DNekMatSharedPtr Mat;
            MatrixType mattype;

            switch(mattype = mkey.GetMatrixType())
            {
714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739
            case ePhysInterpToEquiSpaced:
                {
                    int nq = m_base[0]->GetNumPoints();

                    // take definition from key
                    if(mkey.ConstFactorExists(eFactorConst))
                    {
                        nq = (int) mkey.GetConstFactor(eFactorConst);
                    }

                    int neq = LibUtilities::StdSegData::
                                              getNumberOfCoefficients(nq);
                    Array<OneD, NekDouble >   coords (1);
                    DNekMatSharedPtr          I     ;
                    Mat                     = MemoryManager<DNekMat>::
                                                AllocateSharedPtr(neq, nq);

                    for(int i = 0; i < neq; ++i)
                    {
                        coords[0] = -1.0 + 2*i/(NekDouble)(neq-1);
                        I         = m_base[0]->GetI(coords);
                        Vmath::Vcopy(nq, I->GetRawPtr(), 1,
                                         Mat->GetRawPtr()+i,neq);
                    }
                }
                break;
740 741 742
            case eFwdTrans:
                {
                    Mat = MemoryManager<DNekMat>::AllocateSharedPtr(m_ncoeffs,m_ncoeffs);
743
                    StdMatrixKey iprodkey(eIProductWRTBase,v_DetShapeType(),*this);
744
                    DNekMat &Iprod = *GetStdMatrix(iprodkey);
745
                    StdMatrixKey imasskey(eInvMass,v_DetShapeType(),*this);
746
                    DNekMat &Imass = *GetStdMatrix(imasskey);
Mike Kirby's avatar
Mike Kirby committed
747

748 749 750 751 752 753
                    (*Mat) = Imass*Iprod;
                }
                break;
            default:
                {
                    Mat = StdExpansion::CreateGeneralMatrix(mkey);
Mike Kirby's avatar
Mike Kirby committed
754

755 756 757 758 759 760 761 762 763 764 765 766 767
                    if(mattype ==  eMass)
                    {
                        // For Fourier basis set the imaginary component
                        // of mean mode to have a unit diagonal component
                        // in mass matrix
                        if(m_base[0]->GetBasisType() == LibUtilities::eFourier)
                        {
                            (*Mat)(1,1) = 1.0;
                        }
                    }
                }
                break;
            }
Mike Kirby's avatar
Mike Kirby committed
768

769 770 771
            return Mat;
        }

Pavel Burovskiy's avatar
Pavel Burovskiy committed
772

773 774
        DNekMatSharedPtr StdSegExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
        {
775
            return v_GenMatrix(mkey);
776
        }
Mike Kirby's avatar
Mike Kirby committed
777

Pavel Burovskiy's avatar
Pavel Burovskiy committed
778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796
        //---------------------------------------------------------------------
        // Mappings
        //---------------------------------------------------------------------


        void StdSegExp::v_GetBoundaryMap(Array<OneD, unsigned int>& outarray)
        {
            if(outarray.num_elements() != NumBndryCoeffs())
            {
                outarray = Array<OneD, unsigned int>(NumBndryCoeffs());
            }
            const LibUtilities::BasisType Btype = GetBasisType(0);
            int nummodes = m_base[0]->GetNumModes();

            outarray[0] = 0;

            switch(Btype)
            {
            case LibUtilities::eGLL_Lagrange:
797
            case LibUtilities::eGauss_Lagrange:
Pavel Burovskiy's avatar
Pavel Burovskiy committed
798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823
            case LibUtilities::eChebyshev:
            case LibUtilities::eFourier:
                outarray[1]= nummodes-1;
                break;
            case LibUtilities::eModified_A:
            case LibUtilities::eModified_B:
                outarray[1] = 1;
                break;
            default:
                ASSERTL0(0,"Mapping array is not defined for this expansion");
                break;
            }
        }

        void StdSegExp::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
        {
            int i;
            if(outarray.num_elements()!=GetNcoeffs()-NumBndryCoeffs())
            {
                outarray = Array<OneD, unsigned int>(GetNcoeffs()-NumBndryCoeffs());
            }
            const LibUtilities::BasisType Btype = GetBasisType(0);

            switch(Btype)
            {
            case LibUtilities::eGLL_Lagrange:
824
            case LibUtilities::eGauss_Lagrange:
Pavel Burovskiy's avatar
Pavel Burovskiy committed
825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844
            case LibUtilities::eChebyshev:
            case LibUtilities::eFourier:
                for(i = 0 ; i < GetNcoeffs()-2;i++)
                {
                    outarray[i] = i+1;
                }
                break;
            case LibUtilities::eModified_A:
            case LibUtilities::eModified_B:
                for(i = 0 ; i < GetNcoeffs()-2;i++)
                {
                    outarray[i] = i+2;
                }
                break;
            default:
                ASSERTL0(0,"Mapping array is not defined for this expansion");
                break;
            }
        }

845
        int StdSegExp::v_GetVertexMap(int localVertexId,bool useCoeffPacking)
Pavel Burovskiy's avatar
Pavel Burovskiy committed
846 847 848 849 850 851 852 853 854 855 856 857 858
        {
            ASSERTL0((localVertexId==0)||(localVertexId==1),"local vertex id"
                     "must be between 0 or 1");

            int localDOF = localVertexId;

            if( (m_base[0]->GetBasisType()==LibUtilities::eGLL_Lagrange) &&
                (localVertexId==1) )
            {
                localDOF = m_base[0]->GetNumModes()-1;
            }
            return localDOF;
        }
859

860 861 862 863 864 865 866 867 868 869 870 871 872 873
        void StdSegExp::v_GetSimplexEquiSpacedConnectivity(
            Array<OneD, int> &conn,
            bool              standard)
        {
            int np = m_base[0]->GetNumPoints();

            conn     = Array<OneD, int>(2*(np-1));
            int cnt  = 0;
            for(int i = 0; i < np-1; ++i)
            {
                conn[cnt++] = i;
                conn[cnt++] = i+1;
            }
        }
874 875 876

    }//end namespace
}//end namespace