QuadExp.cpp 93 KB
Newer Older
Mike Kirby's avatar
Mike Kirby committed
1 2
///////////////////////////////////////////////////////////////////////////////
//
3
// File QuadExp.cpp
Mike Kirby's avatar
Mike Kirby committed
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 31
//
// 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.
//
32
// Description: Expansion for quadrilateral elements.
Mike Kirby's avatar
Mike Kirby committed
33 34
//
///////////////////////////////////////////////////////////////////////////////
35
#include <LibUtilities/Foundations/InterpCoeff.h>
Mike Kirby's avatar
Mike Kirby committed
36
#include <LocalRegions/QuadExp.h>
37
#include <LocalRegions/Expansion3D.h>
38 39
#include <LibUtilities/BasicUtils/VmathArray.hpp>
#include <LibUtilities/BasicUtils/Vmath.hpp>
40
#include <LibUtilities/Foundations/Interp.h>
41
#include <LocalRegions/SegExp.h>
Mike Kirby's avatar
Mike Kirby committed
42

43
using namespace std;
Mike Kirby's avatar
Mike Kirby committed
44 45 46

namespace Nektar
{
47 48 49 50
    namespace LocalRegions
    {
        QuadExp::QuadExp(const LibUtilities::BasisKey &Ba,
                         const LibUtilities::BasisKey &Bb,
51
                         const SpatialDomains::QuadGeomSharedPtr &geom):
52 53
             StdExpansion  (Ba.GetNumModes()*Bb.GetNumModes(),2,Ba,Bb),
             StdExpansion2D(Ba.GetNumModes()*Bb.GetNumModes(),Ba,Bb),
54
             StdQuadExp    (Ba,Bb),
55 56
             Expansion     (geom),
             Expansion2D   (geom),
57
             m_matrixManager(
58 59
                 std::bind(&QuadExp::CreateMatrix, this, std::placeholders::_1),
                 std::string("QuadExpMatrix")),
60
             m_staticCondMatrixManager(
61 62
                 std::bind(&QuadExp::CreateStaticCondMatrix, this, std::placeholders::_1),
                 std::string("QuadExpStaticCondMatrix"))
63
        {
64 65
        }

66

67
        QuadExp::QuadExp(const QuadExp &T):
68 69 70
            StdExpansion(T),
            StdExpansion2D(T),
            StdQuadExp(T),
71 72
            Expansion   (T),
            Expansion2D (T),
73 74
            m_matrixManager(T.m_matrixManager),
            m_staticCondMatrixManager(T.m_staticCondMatrixManager)
75 76
        {
        }
77

78

79 80 81
        QuadExp::~QuadExp()
        {
        }
82 83


Daniele de Grazia's avatar
Daniele de Grazia committed
84 85 86
        NekDouble QuadExp::v_Integral(
            const Array<OneD,
            const NekDouble> &inarray)
87
        {
88 89
            int    nquad0 = m_base[0]->GetNumPoints();
            int    nquad1 = m_base[1]->GetNumPoints();
90
            Array<OneD, const NekDouble> jac = m_metricinfo->GetJac(GetPointsKeys());
91
            NekDouble ival;
Peter Vos's avatar
Peter Vos committed
92
            Array<OneD,NekDouble> tmp(nquad0*nquad1);
93

94
            // multiply inarray with Jacobian
Daniele de Grazia's avatar
Daniele de Grazia committed
95
            if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
96
            {
Peter Vos's avatar
Peter Vos committed
97
                Vmath::Vmul(nquad0*nquad1, jac, 1, inarray, 1, tmp, 1);
98 99 100
            }
            else
            {
Peter Vos's avatar
Peter Vos committed
101
                Vmath::Smul(nquad0*nquad1, jac[0], inarray, 1, tmp, 1);
102
            }
103

104
            // call StdQuadExp version;
Andrew Comerford's avatar
Andrew Comerford committed
105
            ival = StdQuadExp::v_Integral(tmp);
106
            return  ival;
107
        }
108

Daniele de Grazia's avatar
Daniele de Grazia committed
109 110 111 112 113
        void QuadExp::v_PhysDeriv(
            const Array<OneD, const NekDouble> & inarray,
                  Array<OneD,NekDouble> &out_d0,
                  Array<OneD,NekDouble> &out_d1,
                  Array<OneD,NekDouble> &out_d2)
114
        {
115 116 117
            int    nquad0 = m_base[0]->GetNumPoints();
            int    nquad1 = m_base[1]->GetNumPoints();
            int     nqtot = nquad0*nquad1;
118
            const Array<TwoD, const NekDouble>& df = m_metricinfo->GetDerivFactors(GetPointsKeys());
119 120
            Array<OneD,NekDouble> diff0(2*nqtot);
            Array<OneD,NekDouble> diff1(diff0+nqtot);
121

122
            StdQuadExp::v_PhysDeriv(inarray, diff0, diff1);
123

Daniele de Grazia's avatar
Daniele de Grazia committed
124
            if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
125
            {
Daniele de Grazia's avatar
Daniele de Grazia committed
126
                if (out_d0.num_elements())
127
                {
128 129 130
                    Vmath::Vmul  (nqtot, df[0], 1, diff0, 1, out_d0, 1);
                    Vmath::Vvtvp (nqtot, df[1], 1, diff1, 1, out_d0, 1, 
                    					 out_d0,1);
131 132 133 134
                }

                if(out_d1.num_elements())
                {
135 136
                    Vmath::Vmul  (nqtot,df[2],1,diff0,1, out_d1, 1);
                    Vmath::Vvtvp (nqtot,df[3],1,diff1,1, out_d1, 1, out_d1,1);
137 138
                }

Daniele de Grazia's avatar
Daniele de Grazia committed
139
                if (out_d2.num_elements())
140
                {
141 142
                    Vmath::Vmul  (nqtot,df[4],1,diff0,1, out_d2, 1);
                    Vmath::Vvtvp (nqtot,df[5],1,diff1,1, out_d2, 1, out_d2,1);
143
                }
144
            }
145
            else // regular geometry
146
            {
Daniele de Grazia's avatar
Daniele de Grazia committed
147
                if (out_d0.num_elements())
148
                {
149 150
                    Vmath::Smul (nqtot, df[0][0], diff0, 1, out_d0, 1);
                    Blas::Daxpy (nqtot, df[1][0], diff1, 1, out_d0, 1);
151
                }
152

Daniele de Grazia's avatar
Daniele de Grazia committed
153
                if (out_d1.num_elements())
154
                {
155 156
                    Vmath::Smul (nqtot, df[2][0], diff0, 1, out_d1, 1);
                    Blas::Daxpy (nqtot, df[3][0], diff1, 1, out_d1, 1);
157
                }
158

Daniele de Grazia's avatar
Daniele de Grazia committed
159
                if (out_d2.num_elements())
160
                {
161 162
                    Vmath::Smul (nqtot, df[4][0], diff0, 1, out_d2, 1);
                    Blas::Daxpy (nqtot, df[5][0], diff1, 1, out_d2, 1);
163 164
                }
            }
165
        }
166

167

Daniele de Grazia's avatar
Daniele de Grazia committed
168 169 170 171
        void QuadExp::v_PhysDeriv(
            const int dir,
            const Array<OneD, const NekDouble>& inarray,
                  Array<OneD, NekDouble> &outarray)
172
        {
Daniele de Grazia's avatar
Daniele de Grazia committed
173
            switch (dir)
174 175 176
            {
            case 0:
                {
Daniele de Grazia's avatar
Daniele de Grazia committed
177 178
                    v_PhysDeriv(inarray, outarray, NullNekDouble1DArray,
                                NullNekDouble1DArray);
179 180 181 182
                }
                break;
            case 1:
                {
Daniele de Grazia's avatar
Daniele de Grazia committed
183 184
                    v_PhysDeriv(inarray, NullNekDouble1DArray, outarray,
                                NullNekDouble1DArray);
185 186 187 188
                }
                break;
            case 2:
                {
Daniele de Grazia's avatar
Daniele de Grazia committed
189 190
                    v_PhysDeriv(inarray, NullNekDouble1DArray,
                                NullNekDouble1DArray, outarray);
191 192 193 194 195 196 197
                }
                break;
            default:
                {
                    ASSERTL1(false,"input dir is out of range");
                }
                break;
198
            }
199
        }
200

201

Daniele de Grazia's avatar
Daniele de Grazia committed
202 203 204 205
        void QuadExp::v_PhysDirectionalDeriv(
            const Array<OneD, const NekDouble> & inarray,
            const Array<OneD, const NekDouble>& direction,
                  Array<OneD,NekDouble> &out)
206
        {
207 208 209
            int    nquad0 = m_base[0]->GetNumPoints();
            int    nquad1 = m_base[1]->GetNumPoints();
            int    nqtot = nquad0*nquad1;
210

211
            const Array<TwoD, const NekDouble>& df = m_metricinfo->GetDerivFactors(GetPointsKeys());
212

213 214 215 216 217
            Array<OneD,NekDouble> diff0(2*nqtot);
            Array<OneD,NekDouble> diff1(diff0+nqtot);

            StdQuadExp::v_PhysDeriv(inarray, diff0, diff1);

Daniele de Grazia's avatar
Daniele de Grazia committed
218
            if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
219
            {
220 221 222 223 224 225 226 227
                Array<OneD, Array<OneD, NekDouble> > tangmat(2);

                // d/dx_v^s = v_x*ds/dx + v_y*ds/dy + v_z*dx/dz
                for (int i=0; i< 2; ++i)
                {
                    tangmat[i] = Array<OneD, NekDouble>(nqtot,0.0);
                    for (int k=0; k<(m_geom->GetCoordim()); ++k)
                    {
Daniele de Grazia's avatar
Daniele de Grazia committed
228
                        Vmath::Vvtvp(nqtot,
229
                                     &df[2*k+i][0], 1,
Daniele de Grazia's avatar
Daniele de Grazia committed
230 231 232
                                     &direction[k*nqtot], 1,
                                     &tangmat[i][0], 1,
                                     &tangmat[i][0], 1);
233 234 235 236
                    }
                }

                /// D_v = d/dx_v^s + d/dx_v^r
Daniele de Grazia's avatar
Daniele de Grazia committed
237
                if (out.num_elements())
238
                {
Daniele de Grazia's avatar
Daniele de Grazia committed
239 240 241 242 243 244 245 246 247
                    Vmath::Vmul  (nqtot,
                                  &tangmat[0][0], 1,
                                  &diff0[0], 1,
                                  &out[0], 1);
                    Vmath::Vvtvp (nqtot,
                                  &tangmat[1][0], 1,
                                  &diff1[0], 1,
                                  &out[0], 1,
                                  &out[0], 1);
248
                }
249

250 251
            }
            else
252
            {
Daniele de Grazia's avatar
Daniele de Grazia committed
253 254
                ASSERTL1(m_metricinfo->GetGtype() ==
                         SpatialDomains::eDeformed,"Wrong route");
255 256
            }
        }
257

258

Daniele de Grazia's avatar
Daniele de Grazia committed
259 260
        void QuadExp::v_FwdTrans(
            const Array<OneD, const NekDouble> & inarray,
261
            Array<OneD,NekDouble> &outarray)
262
        {
Daniele de Grazia's avatar
Daniele de Grazia committed
263
            if ((m_base[0]->Collocation())&&(m_base[1]->Collocation()))
264
            {
265
                Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
266 267 268 269
            }
            else
            {
                IProductWRTBase(inarray,outarray);
270

271 272
                // get Mass matrix inverse
                MatrixKey             masskey(StdRegions::eInvMass,
273
                                              DetShapeType(),*this);
274
                DNekScalMatSharedPtr  matsys = m_matrixManager[masskey];
275 276 277 278 279 280 281 282 283

                // copy inarray in case inarray == outarray
                NekVector<NekDouble> in(m_ncoeffs,outarray,eCopy);
                NekVector<NekDouble> out(m_ncoeffs,outarray,eWrapper);

                out = (*matsys)*in;
            }
        }

284

Daniele de Grazia's avatar
Daniele de Grazia committed
285 286 287
        void QuadExp::v_FwdTrans_BndConstrained(
            const Array<OneD, const NekDouble>& inarray,
                  Array<OneD, NekDouble> &outarray)
288
        {
Daniele de Grazia's avatar
Daniele de Grazia committed
289
            if ((m_base[0]->Collocation())&&(m_base[1]->Collocation()))
290 291
            {
                Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
292 293 294
            }
            else
            {
295 296 297 298 299
                int i,j;
                int npoints[2] = {m_base[0]->GetNumPoints(),
                                  m_base[1]->GetNumPoints()};
                int nmodes[2]  = {m_base[0]->GetNumModes(),
                                  m_base[1]->GetNumModes()};
300

301
                fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 );
302

303 304 305 306 307 308
                if(nmodes[0] == 1 && nmodes[1] == 1)
                {
                    outarray[0] = inarray[0];
                    return;
                }

309 310
                Array<OneD, NekDouble> physEdge[4];
                Array<OneD, NekDouble> coeffEdge[4];
311
                StdRegions::Orientation orient[4];
Daniele de Grazia's avatar
Daniele de Grazia committed
312
                for (i = 0; i < 4; i++)
313
                {
314 315 316
                    physEdge[i]  = Array<OneD, NekDouble>(npoints[i%2]);
                    coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i%2]);
                    orient[i]    = GetEorient(i);
317
                }
318

Daniele de Grazia's avatar
Daniele de Grazia committed
319
                for (i = 0; i < npoints[0]; i++)
320
                {
321
                    physEdge[0][i] = inarray[i];
Dave Moxey's avatar
Dave Moxey committed
322
                    physEdge[2][i] = inarray[npoints[0]*(npoints[1]-1)+i];
323 324
                }

Daniele de Grazia's avatar
Daniele de Grazia committed
325
                for (i = 0; i < npoints[1]; i++)
326
                {
Daniele de Grazia's avatar
Daniele de Grazia committed
327 328 329
                    physEdge[1][i] =
                        inarray[npoints[0]-1+i*npoints[0]];
                    physEdge[3][i] =
Dave Moxey's avatar
Dave Moxey committed
330
                        inarray[i*npoints[0]];
331
                }
332

Daniele de Grazia's avatar
Daniele de Grazia committed
333
                for (i = 0; i < 4; i++)
334
                {
Daniele de Grazia's avatar
Daniele de Grazia committed
335
                    if ( orient[i] == StdRegions::eBackwards )
336
                    {
Daniele de Grazia's avatar
Daniele de Grazia committed
337 338
                        reverse((physEdge[i]).get(),
                                (physEdge[i]).get() + npoints[i%2] );
339
                    }
340
                }
341

342
                SegExpSharedPtr segexp[4];
Daniele de Grazia's avatar
Daniele de Grazia committed
343
                for (i = 0; i < 4; i++)
344
                {
Daniele de Grazia's avatar
Daniele de Grazia committed
345 346 347
                    segexp[i] = MemoryManager<LocalRegions::SegExp>::
                        AllocateSharedPtr(
                            m_base[i%2]->GetBasisKey(),GetGeom2D()->GetEdge(i));
348
                }
349

350 351 352
                Array<OneD, unsigned int> mapArray;
                Array<OneD, int>          signArray;
                NekDouble sign;
353

Daniele de Grazia's avatar
Daniele de Grazia committed
354
                for (i = 0; i < 4; i++)
355
                {
Daniele de Grazia's avatar
Daniele de Grazia committed
356 357
                    segexp[i%2]->FwdTrans_BndConstrained(
                        physEdge[i],coeffEdge[i]);
358

359
                    GetEdgeToElementMap(i,orient[i],mapArray,signArray);
Daniele de Grazia's avatar
Daniele de Grazia committed
360
                    for (j=0; j < nmodes[i%2]; j++)
361
                    {
362 363
                        sign = (NekDouble) signArray[j];
                        outarray[ mapArray[j] ] = sign * coeffEdge[i][j];
364
                    }
365
                }
366

367 368 369 370
                int nBoundaryDofs = NumBndryCoeffs();
                int nInteriorDofs = m_ncoeffs - nBoundaryDofs;

                if (nInteriorDofs > 0) {
371 372
                    Array<OneD, NekDouble> tmp0(m_ncoeffs);
                    Array<OneD, NekDouble> tmp1(m_ncoeffs);
373

Daniele de Grazia's avatar
Daniele de Grazia committed
374 375
                    StdRegions::StdMatrixKey
                        stdmasskey(StdRegions::eMass,DetShapeType(),*this);
376 377
                    MassMatrixOp(outarray,tmp0,stdmasskey);
                    IProductWRTBase(inarray,tmp1);
378

379
                    Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
380

381 382 383
                    // get Mass matrix inverse (only of interior DOF)
                    // use block (1,1) of the static condensed system
                    // note: this block alreay contains the inverse matrix
Daniele de Grazia's avatar
Daniele de Grazia committed
384 385 386 387
                    MatrixKey
                        masskey(StdRegions::eMass,DetShapeType(),*this);
                    DNekScalMatSharedPtr matsys =
                        (m_staticCondMatrixManager[masskey])->GetBlock(1,1);
388

389 390 391 392 393
                    Array<OneD, NekDouble> rhs(nInteriorDofs);
                    Array<OneD, NekDouble> result(nInteriorDofs);

                    GetInteriorMap(mapArray);

Daniele de Grazia's avatar
Daniele de Grazia committed
394
                    for (i = 0; i < nInteriorDofs; i++)
395
                    {
396
                        rhs[i] = tmp1[ mapArray[i] ];
397
                    }
398

Daniele de Grazia's avatar
Daniele de Grazia committed
399 400 401
                    Blas::Dgemv('N', nInteriorDofs, nInteriorDofs,
                                matsys->Scale(),
                                &((matsys->GetOwnedMatrix())->GetPtr())[0],
402 403
                                nInteriorDofs,rhs.get(),1,0.0,result.get(),1);

Daniele de Grazia's avatar
Daniele de Grazia committed
404
                    for (i = 0; i < nInteriorDofs; i++)
405
                    {
406 407 408 409
                        outarray[ mapArray[i] ] = result[i];
                    }
                }
            }
410

411
        }
412

413

414
        void QuadExp::v_IProductWRTBase(
Daniele de Grazia's avatar
Daniele de Grazia committed
415 416
            const Array<OneD, const NekDouble>& inarray,
                  Array<OneD, NekDouble> &outarray)
417
        {
Daniele de Grazia's avatar
Daniele de Grazia committed
418
            if (m_base[0]->Collocation() && m_base[1]->Collocation())
419 420
            {
                MultiplyByQuadratureMetric(inarray,outarray);
421 422 423
            }
            else
            {
424
                IProductWRTBase_SumFac(inarray,outarray);
425 426 427
            }
        }

428

Daniele de Grazia's avatar
Daniele de Grazia committed
429 430 431 432
        void QuadExp::v_IProductWRTDerivBase(
            const int dir,
            const Array<OneD, const NekDouble>& inarray,
                  Array<OneD, NekDouble> & outarray)
433
        {
434 435
            IProductWRTDerivBase_SumFac(dir,inarray,outarray);
        }
436

437

Daniele de Grazia's avatar
Daniele de Grazia committed
438 439
        void QuadExp::v_IProductWRTBase_SumFac(
            const Array<OneD, const NekDouble>& inarray,
440 441
            Array<OneD, NekDouble> &outarray,
                                               bool multiplybyweights)
442 443 444 445
        {
            int    nquad0 = m_base[0]->GetNumPoints();
            int    nquad1 = m_base[1]->GetNumPoints();
            int    order0 = m_base[0]->GetNumModes();
446

447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464
            if(multiplybyweights)
            {
                Array<OneD,NekDouble> tmp(nquad0*nquad1+nquad1*order0);
                Array<OneD,NekDouble> wsp(tmp+nquad0*nquad1);
                
                MultiplyByQuadratureMetric(inarray,tmp);
                StdQuadExp::IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
                                                         m_base[1]->GetBdata(),
                                                         tmp,outarray,wsp,true,true);
            }
            else
            {
                Array<OneD,NekDouble> wsp(nquad1*order0);
                
                StdQuadExp::IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
                                                         m_base[1]->GetBdata(),
                                                         inarray,outarray,wsp,true,true);
            }
465
        }
466

467

Daniele de Grazia's avatar
Daniele de Grazia committed
468 469 470
        void QuadExp::v_IProductWRTBase_MatOp(
             const Array<OneD, const NekDouble>& inarray,
                   Array<OneD, NekDouble> &outarray)
471 472
        {
            int nq = GetTotPoints();
Daniele de Grazia's avatar
Daniele de Grazia committed
473 474
            MatrixKey
                iprodmatkey(StdRegions::eIProductWRTBase,DetShapeType(),*this);
475
            DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
476

Daniele de Grazia's avatar
Daniele de Grazia committed
477 478
            Blas::Dgemv('N',m_ncoeffs,nq,iprodmat->Scale(),
                        (iprodmat->GetOwnedMatrix())->GetPtr().get(),
479
                        m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
480 481 482

        }

483

Daniele de Grazia's avatar
Daniele de Grazia committed
484 485 486 487
        void QuadExp::v_IProductWRTDerivBase_SumFac(
            const int dir,
            const Array<OneD, const NekDouble>& inarray,
                  Array<OneD, NekDouble> & outarray)
488
        {
Daniele de Grazia's avatar
Daniele de Grazia committed
489 490 491 492
            ASSERTL1((dir==0) || (dir==1) || (dir==2),
                     "Invalid direction.");
            ASSERTL1((dir==2) ? (m_geom->GetCoordim() ==3):true,
                     "Invalid direction.");
493

494 495 496 497
            int    nquad0  = m_base[0]->GetNumPoints();
            int    nquad1  = m_base[1]->GetNumPoints();
            int    nqtot   = nquad0*nquad1;
            int    nmodes0 = m_base[0]->GetNumModes();
498

499
            const Array<TwoD, const NekDouble>& df = m_metricinfo->GetDerivFactors(GetPointsKeys());
500

501 502 503 504
            Array<OneD, NekDouble> tmp1(2*nqtot+m_ncoeffs+nmodes0*nquad1);
            Array<OneD, NekDouble> tmp2(tmp1 +   nqtot);
            Array<OneD, NekDouble> tmp3(tmp1 + 2*nqtot);
            Array<OneD, NekDouble> tmp4(tmp1 + 2*nqtot+m_ncoeffs);
505

Daniele de Grazia's avatar
Daniele de Grazia committed
506
            if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
507
            {
Daniele de Grazia's avatar
Daniele de Grazia committed
508
                Vmath::Vmul(nqtot,
509
                            &df[2*dir][0], 1,
Daniele de Grazia's avatar
Daniele de Grazia committed
510 511 512
                            inarray.get(), 1,
                            tmp1.get(), 1);
                Vmath::Vmul(nqtot,
513
                            &df[2*dir+1][0], 1,
Daniele de Grazia's avatar
Daniele de Grazia committed
514 515
                            inarray.get(), 1,
                            tmp2.get(),1);
516
            }
517
            else
518
            {
Daniele de Grazia's avatar
Daniele de Grazia committed
519
                Vmath::Smul(nqtot,
520
                            df[2*dir][0], inarray.get(), 1,
Daniele de Grazia's avatar
Daniele de Grazia committed
521 522
                            tmp1.get(), 1);
                Vmath::Smul(nqtot,
523
                            df[2*dir+1][0], inarray.get(), 1,
Daniele de Grazia's avatar
Daniele de Grazia committed
524
                            tmp2.get(), 1);
525
            }
526

527 528
            MultiplyByQuadratureMetric(tmp1,tmp1);
            MultiplyByQuadratureMetric(tmp2,tmp2);
529

Daniele de Grazia's avatar
Daniele de Grazia committed
530 531 532 533 534 535
            IProductWRTBase_SumFacKernel(
                m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
                tmp1, tmp3, tmp4, false, true);
            IProductWRTBase_SumFacKernel(
                m_base[0]->GetBdata() , m_base[1]->GetDbdata(),
                tmp2, outarray, tmp4, true, false);
536
            Vmath::Vadd(m_ncoeffs, tmp3, 1, outarray, 1, outarray, 1);
537
        }
538

539

Daniele de Grazia's avatar
Daniele de Grazia committed
540 541 542 543
        void QuadExp::v_IProductWRTDerivBase_MatOp(
            const int dir,
            const Array<OneD, const NekDouble>& inarray,
                  Array<OneD, NekDouble> &outarray)
544
        {
545
            int nq = GetTotPoints();
546
            StdRegions::MatrixType mtype = StdRegions::eIProductWRTDerivBase0;
547

Daniele de Grazia's avatar
Daniele de Grazia committed
548
            switch (dir)
549
            {
Daniele de Grazia's avatar
Daniele de Grazia committed
550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570
                case 0:
                    {
                        mtype = StdRegions::eIProductWRTDerivBase0;
                    }
                    break;
                case 1:
                    {
                        mtype = StdRegions::eIProductWRTDerivBase1;
                    }
                    break;
                case 2:
                    {
                        mtype = StdRegions::eIProductWRTDerivBase2;
                    }
                    break;
                default:
                    {
                        ASSERTL1(false,"input dir is out of range");
                    }
                    break;
            }   
571

572
            MatrixKey      iprodmatkey(mtype,DetShapeType(),*this);
573
            DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
574

Daniele de Grazia's avatar
Daniele de Grazia committed
575 576
            Blas::Dgemv('N', m_ncoeffs, nq, iprodmat->Scale(),
                        (iprodmat->GetOwnedMatrix())->GetPtr().get(),
577
                        m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
578
        }
579

580

581
        void QuadExp::v_NormVectorIProductWRTBase(
Dave Moxey's avatar
Dave Moxey committed
582 583 584 585
            const Array<OneD, const NekDouble> &Fx,
            const Array<OneD, const NekDouble> &Fy,
            const Array<OneD, const NekDouble> &Fz,
                  Array<OneD,       NekDouble> &outarray)
586
        {
587
            int nq = m_base[0]->GetNumPoints()*m_base[1]->GetNumPoints();
Dave Moxey's avatar
Dave Moxey committed
588 589 590 591 592 593 594 595
            Array<OneD, NekDouble> Fn(nq);
            
            const Array<OneD, const Array<OneD, NekDouble> > &normals = 
                GetLeftAdjacentElementExp()->GetFaceNormal(
                    GetLeftAdjacentElementFace());
            
            if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
            {
596 597 598
                Vmath::Vvtvvtp(nq,&normals[0][0],1,&Fx[0],1,
                                  &normals[1][0],1,&Fy[0],1,&Fn[0],1);
                Vmath::Vvtvp  (nq,&normals[2][0],1,&Fz[0],1,&Fn[0],1,&Fn[0],1);
Dave Moxey's avatar
Dave Moxey committed
599 600 601
            }
            else
            {
602 603 604
                Vmath::Svtsvtp(nq,normals[0][0],&Fx[0],1,
                                  normals[1][0],&Fy[0],1,&Fn[0],1);
                Vmath::Svtvp  (nq,normals[2][0],&Fz[0],1,&Fn[0],1,&Fn[0],1);
Dave Moxey's avatar
Dave Moxey committed
605
            }
606

607
            IProductWRTBase(Fn,outarray);
608
        }
609

Chris Cantwell's avatar
Chris Cantwell committed
610 611 612
        void QuadExp::v_NormVectorIProductWRTBase(
            const Array<OneD, const Array<OneD, NekDouble> > &Fvec,
                  Array<OneD,       NekDouble>               &outarray)
613 614 615 616
        {
            NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
        }

617 618 619 620 621 622
        StdRegions::StdExpansionSharedPtr QuadExp::v_GetStdExp(void) const
        {
            return MemoryManager<StdRegions::StdQuadExp>
                ::AllocateSharedPtr(m_base[0]->GetBasisKey(),
                                    m_base[1]->GetBasisKey());
        }
Sehun Chun's avatar
Sehun Chun committed
623

624 625 626 627 628 629 630 631 632 633 634 635

        StdRegions::StdExpansionSharedPtr QuadExp::v_GetLinStdExp(void) const
        {
            LibUtilities::BasisKey bkey0(m_base[0]->GetBasisType(),
                           2, m_base[0]->GetPointsKey());
            LibUtilities::BasisKey bkey1(m_base[1]->GetBasisType(),
                           2, m_base[1]->GetPointsKey());

            return MemoryManager<StdRegions::StdQuadExp>
                ::AllocateSharedPtr( bkey0, bkey1);
        }

Daniele de Grazia's avatar
Daniele de Grazia committed
636
        void QuadExp::v_GetCoords(
637 638 639
            Array<OneD, NekDouble> &coords_0,
            Array<OneD, NekDouble> &coords_1,
            Array<OneD, NekDouble> &coords_2)
640
        {
641
            Expansion::v_GetCoords(coords_0, coords_1, coords_2);
642
        }
643 644


645
        void QuadExp::v_GetCoord(const Array<OneD, const NekDouble> &Lcoords,
646
                               Array<OneD,NekDouble> &coords)
647
        {
648
            int  i;
649 650

            ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 &&
651 652
                     Lcoords[1] >= -1.0 && Lcoords[1]  <=1.0,
                     "Local coordinates are not in region [-1,1]");
653 654

            m_geom->FillGeom();
Daniele de Grazia's avatar
Daniele de Grazia committed
655
            for (i = 0; i < m_geom->GetCoordim(); ++i)
656 657 658
            {
                coords[i] = m_geom->GetCoord(i,Lcoords);
            }
659
        }
660

661

662

Chris Cantwell's avatar
Chris Cantwell committed
663
        /**
664 665 666 667 668
         * Given the local cartesian coordinate \a Lcoord evaluate the
         * value of physvals at this point by calling through to the
         * StdExpansion method
         */
        NekDouble QuadExp::v_StdPhysEvaluate(
Dave Moxey's avatar
Dave Moxey committed
669 670
            const Array<OneD, const NekDouble> &Lcoord,
            const Array<OneD, const NekDouble> &physvals)
671 672 673 674 675
        {
            // Evaluate point in local (eta) coordinates.
            return StdQuadExp::v_PhysEvaluate(Lcoord,physvals);
        }

Daniele de Grazia's avatar
Daniele de Grazia committed
676
        NekDouble QuadExp::v_PhysEvaluate(
677 678
            const Array<OneD, const NekDouble> &coord,
            const Array<OneD, const NekDouble> &physvals)
679
        {
680
            Array<OneD,NekDouble> Lcoord = Array<OneD, NekDouble>(2);
681

682 683
            ASSERTL0(m_geom,"m_geom not defined");
            m_geom->GetLocCoords(coord,Lcoord);
684

685 686
            return StdQuadExp::v_PhysEvaluate(Lcoord, physvals);
        }
687

688

689 690 691
        // Get edge values from the 2D Phys space along an edge
        // following a counter clockwise edge convention for definition
        // of edgedir, Note that point distribution is given by QuadExp.
Daniele de Grazia's avatar
Daniele de Grazia committed
692 693 694 695
        void QuadExp::v_GetEdgePhysVals(
            const int edge,
            const Array<OneD, const NekDouble> &inarray,
                  Array<OneD,NekDouble> &outarray)
696 697 698
        {
            int nquad0 = m_base[0]->GetNumPoints();
            int nquad1 = m_base[1]->GetNumPoints();
699

700
            StdRegions::Orientation edgedir = GetEorient(edge);
701 702
            switch(edge)
            {
703 704 705 706 707 708 709 710 711 712 713 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 740 741 742 743 744 745 746 747 748
            case 0:
                if (edgedir == StdRegions::eForwards)
                {
                    Vmath::Vcopy(nquad0,&(inarray[0]),1,&(outarray[0]),1);
                }
                else
                {
                    Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1),-1,
                                 &(outarray[0]),1);
                }
                break;
            case 1:
                if (edgedir == StdRegions::eForwards)
                {
                    Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0-1),nquad0,
                                 &(outarray[0]),1);
                }
                else
                {
                    Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0*nquad1-1),
                                 -nquad0, &(outarray[0]),1);
                }
                break;
            case 2:
                if (edgedir == StdRegions::eForwards)
                {
                    Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1),-1,
                                 &(outarray[0]),1);
                }
                else
                {
                    Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1),1,
                                 &(outarray[0]),1);
                }
                break;
            case 3:
                if (edgedir == StdRegions::eForwards)
                {
                    Vmath::Vcopy(nquad1,&(inarray[0]) + nquad0*(nquad1-1),
                                 -nquad0,&(outarray[0]),1);
                }
                else
                {
                    Vmath::Vcopy(nquad1,&(inarray[0]),nquad0,
                                 &(outarray[0]),1);
                }
749 750 751 752 753 754
                break;
            default:
                ASSERTL0(false,"edge value (< 3) is out of range");
                break;
            }
        }
755 756
        
        
757 758 759 760 761 762
        void QuadExp::v_GetTracePhysVals(
             const int edge,
             const StdRegions::StdExpansionSharedPtr &EdgeExp,
             const Array<OneD, const NekDouble> &inarray,
             Array<OneD,NekDouble> &outarray,
             StdRegions::Orientation  orient)
763 764 765 766 767
        {
            v_GetEdgePhysVals(edge,EdgeExp,inarray,outarray);
        }


Daniele de Grazia's avatar
Daniele de Grazia committed
768 769 770 771 772
        void QuadExp::v_GetEdgePhysVals(
            const int edge,
            const StdRegions::StdExpansionSharedPtr &EdgeExp,
            const Array<OneD, const NekDouble> &inarray,
                  Array<OneD,NekDouble> &outarray)
773 774 775
        {
            int nquad0 = m_base[0]->GetNumPoints();
            int nquad1 = m_base[1]->GetNumPoints();
Daniele de Grazia's avatar
Daniele de Grazia committed
776 777
            
            // Implementation for all the basis except Gauss points
Daniele de Grazia's avatar
Daniele de Grazia committed
778 779 780 781
            if (m_base[0]->GetPointsType() !=
                LibUtilities::eGaussGaussLegendre &&
                m_base[1]->GetPointsType() !=
                LibUtilities::eGaussGaussLegendre)
782
            {
Daniele de Grazia's avatar
Daniele de Grazia committed
783
                switch (edge)
Daniele de Grazia's avatar
Daniele de Grazia committed
784 785 786 787 788
                {
                    case 0:
                        Vmath::Vcopy(nquad0,&(inarray[0]),1,&(outarray[0]),1);
                        break;
                    case 1:
789 790
                        Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0-1),
                                     nquad0,&(outarray[0]),1);
Daniele de Grazia's avatar
Daniele de Grazia committed
791 792 793 794 795 796
                        break;
                    case 2:
                        Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1),1,
                                     &(outarray[0]),1);
                        break;
                    case 3:
797 798
                        Vmath::Vcopy(nquad1,&(inarray[0]),nquad0,
                                     &(outarray[0]),1);
Daniele de Grazia's avatar
Daniele de Grazia committed
799 800 801 802 803
                        break;
                    default:
                        ASSERTL0(false,"edge value (< 3) is out of range");
                        break;
                }
804
            }
Daniele de Grazia's avatar
Daniele de Grazia committed
805 806 807 808 809
            else
            {
                QuadExp::v_GetEdgeInterpVals(edge, inarray, outarray);
            }
            
810
            // Interpolate if required
Daniele de Grazia's avatar
Daniele de Grazia committed
811 812
            if (m_base[edge%2]->GetPointsKey() !=
                EdgeExp->GetBasis(0)->GetPointsKey())
813 814
            {
                Array<OneD,NekDouble> outtmp(max(nquad0,nquad1));
815
				
816
                outtmp = outarray;
817
				
Daniele de Grazia's avatar
Daniele de Grazia committed
818
                LibUtilities::Interp1D(
Gianmarco Mengaldo's avatar
Minor.  
Gianmarco Mengaldo committed
819 820
                    m_base[edge%2]->GetPointsKey(), outtmp,
                    EdgeExp->GetBasis(0)->GetPointsKey(), outarray);
821
            }
Daniele de Grazia's avatar
Daniele de Grazia committed
822
            
823
            //Reverse data if necessary
824
            if(GetEorient(edge) == StdRegions::eBackwards)
825
            {
Gianmarco Mengaldo's avatar
Minor.  
Gianmarco Mengaldo committed
826 827
                Vmath::Reverse(EdgeExp->GetNumPoints(0),&outarray[0], 1,
                               &outarray[0], 1);
828 829
            }
        }
830
        
Daniele de Grazia's avatar
Daniele de Grazia committed
831 832 833 834 835 836 837
        void QuadExp::v_GetEdgeInterpVals(const int edge,
                    const Array<OneD, const NekDouble> &inarray,
                    Array<OneD,NekDouble> &outarray)
        {
             int i;
             int nq0 = m_base[0]->GetNumPoints();
             int nq1 = m_base[1]->GetNumPoints();
838 839 840 841 842

             StdRegions::ConstFactorMap factors;
             factors[StdRegions::eFactorGaussEdge] = edge;

             StdRegions::StdMatrixKey key(
843
                 StdRegions::eInterpGauss,
844 845 846 847
                 DetShapeType(),*this,factors);
             
             DNekScalMatSharedPtr mat_gauss = m_matrixManager[key];

Daniele de Grazia's avatar
Daniele de Grazia committed
848
             switch (edge)
Daniele de Grazia's avatar
Daniele de Grazia committed
849
             {
850
                 case 0:
Daniele de Grazia's avatar
Daniele de Grazia committed
851
                 {
852
                     for (i = 0; i < nq0; i++)
Daniele de Grazia's avatar
Daniele de Grazia committed
853
                     {
854 855 856
                         outarray[i] = Blas::Ddot(
                             nq1, mat_gauss->GetOwnedMatrix()->GetPtr().get(),
                             1, &inarray[i], nq0);
Daniele de Grazia's avatar
Daniele de Grazia committed
857
                     }
858
                     break;
Daniele de Grazia's avatar
Daniele de Grazia committed
859
                 }
860
                 case 1:
Daniele de Grazia's avatar
Daniele de Grazia committed
861
                 {
862
                     for (i = 0; i < nq1; i++)
Daniele de Grazia's avatar
Daniele de Grazia committed
863
                     {
864 865 866
                         outarray[i] =  Blas::Ddot(
                             nq0, mat_gauss->GetOwnedMatrix()->GetPtr().get(),
                             1, &inarray[i * nq0], 1);
Daniele de Grazia's avatar
Daniele de Grazia committed
867
                     }
868
                     break;
Daniele de Grazia's avatar
Daniele de Grazia committed
869
                 }
870
                 case 2:
Daniele de Grazia's avatar
Daniele de Grazia committed
871
                 {
872
                     for (i = 0; i < nq0; i++)
Daniele de Grazia's avatar
Daniele de Grazia committed
873
                     {
874 875 876
                         outarray[i] = Blas::Ddot(
                             nq1, mat_gauss->GetOwnedMatrix()->GetPtr().get(),
                             1, &inarray[i], nq0);
Daniele de Grazia's avatar
Daniele de Grazia committed
877
                     }
878
                     break;
Daniele de Grazia's avatar
Daniele de Grazia committed
879
                 }
880
                 case 3:
Daniele de Grazia's avatar
Daniele de Grazia committed
881
                 {
882
                     for (i = 0; i < nq1; i++)
Daniele de Grazia's avatar
Daniele de Grazia committed
883
                     {
884 885 886
                         outarray[i] =  Blas::Ddot(
                             nq0, mat_gauss->GetOwnedMatrix()->GetPtr().get(),
                             1, &inarray[i * nq0], 1);
Daniele de Grazia's avatar
Daniele de Grazia committed
887
                     }
888
                     break;
Daniele de Grazia's avatar
Daniele de Grazia committed
889
                 }
890
                 default:
Gianmarco Mengaldo's avatar
Minor.  
Gianmarco Mengaldo committed
891
                     ASSERTL0(false, "edge value (< 3) is out of range");
892
                     break;
Daniele de Grazia's avatar
Daniele de Grazia committed
893 894
             }
        }
895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917
        
        
        void QuadExp::v_GetEdgePhysMap(
            const int                edge,
            Array<OneD, int>        &outarray)
        {
            int nquad0 = m_base[0]->GetNumPoints();
            int nquad1 = m_base[1]->GetNumPoints();
            
            // Get points in Cartesian orientation
            switch (edge)
            {
                case 0:
                    outarray = Array<OneD, int>(nquad0);
                    for (int i = 0; i < nquad0; ++i)
                    {
                        outarray[i] = i;
                    }
                    break;
                case 1:
                    outarray = Array<OneD, int>(nquad1);
                    for (int i = 0; i < nquad1; ++i)
                    {
918
                        outarray[i] = (nquad0-1) + i*nquad0;
919 920 921 922 923 924 925 926 927 928 929 930 931
                    }
                    break;
                case 2:
                    outarray = Array<OneD, int>(nquad0);
                    for (int i = 0; i < nquad0; ++i)
                    {
                        outarray[i] = i + nquad0*(nquad1-1);
                    }
                    break;
                case 3:
                    outarray = Array<OneD, int>(nquad1);
                    for (int i = 0; i < nquad1; ++i)
                    {
932
                        outarray[i] = i*nquad0;
933 934 935 936 937 938 939 940
                    }
                    break;
                default:
                    ASSERTL0(false, "edge value (< 3) is out of range");
                    break;
            }
            
        }
Daniele de Grazia's avatar
Daniele de Grazia committed
941
    
942 943 944
        
        
        
945
        void QuadExp::v_GetEdgeQFactors(
Daniele de Grazia's avatar
Daniele de Grazia committed
946
                const int edge,
947 948 949 950 951
                Array<OneD, NekDouble> &outarray)
        {
            int i;
            int nquad0 = m_base[0]->GetNumPoints();
            int nquad1 = m_base[1]->GetNumPoints();
952

953 954 955
            LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
            const Array<OneD, const NekDouble>& jac = m_metricinfo->GetJac(ptsKeys);
            const Array<TwoD, const NekDouble>& df  = m_metricinfo->GetDerivFactors(ptsKeys);
956 957 958 959 960 961 962 963
            
            Array<OneD, NekDouble> j (max(nquad0, nquad1), 0.0);
            Array<OneD, NekDouble> g0(max(nquad0, nquad1), 0.0);
            Array<OneD, NekDouble> g1(max(nquad0, nquad1), 0.0);
            Array<OneD, NekDouble> g2(max(nquad0, nquad1), 0.0);
            Array<OneD, NekDouble> g3(max(nquad0, nquad1), 0.0);
            
            if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
Daniele de Grazia's avatar
Daniele de Grazia committed