StdPyrExp.cpp 66.1 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 StdPyrExp.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
//
32
// Description: pyramidic routines built upon StdExpansion3D
Mike Kirby's avatar
Mike Kirby committed
33
34
35
36
//
///////////////////////////////////////////////////////////////////////////////

#include <StdRegions/StdPyrExp.h>
Dave Moxey's avatar
Dave Moxey committed
37
38
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <iomanip>
Dave Moxey's avatar
Dave Moxey committed
39

40
41
using namespace std;

Mike Kirby's avatar
Mike Kirby committed
42
43
namespace Nektar
{
Peter Vos's avatar
Peter Vos committed
44
45
    namespace StdRegions
    {
mike turner's avatar
mike turner committed
46
        StdPyrExp::StdPyrExp() // Deafult construct of standard expansion directly called.
47
48
        {
        }
mike turner's avatar
mike turner committed
49

50
51
        StdPyrExp::StdPyrExp(const LibUtilities::BasisKey &Ba,
                             const LibUtilities::BasisKey &Bb,
mike turner's avatar
mike turner committed
52
                             const LibUtilities::BasisKey &Bc)
53
54
55
56
            : StdExpansion  (LibUtilities::StdPyrData::getNumberOfCoefficients(
                                 Ba.GetNumModes(),
                                 Bb.GetNumModes(),
                                 Bc.GetNumModes()),
57
                             3, Ba, Bb, Bc),
58
59
60
61
              StdExpansion3D(LibUtilities::StdPyrData::getNumberOfCoefficients(
                                 Ba.GetNumModes(),
                                 Bb.GetNumModes(),
                                 Bc.GetNumModes()),
62
                             Ba, Bb, Bc)
63
        {
64

65
66
67
68
            ASSERTL0(Ba.GetNumModes() <= Bc.GetNumModes(), "order in 'a' direction is higher "
                     "than order in 'c' direction");
            ASSERTL0(Bb.GetNumModes() <= Bc.GetNumModes(), "order in 'b' direction is higher "
                     "than order in 'c' direction");
69
70
71
72
            ASSERTL1(Bc.GetBasisType() == LibUtilities::eModifiedPyr_C ||
                     Bc.GetBasisType() == LibUtilities::eOrthoPyr_C,
                     "Expected basis type in 'c' direction to be ModifiedPyr_C or OrthoPyr_C");
            
73
74
75
        }

        StdPyrExp::StdPyrExp(const StdPyrExp &T)
76
            : StdExpansion  (T),
77
              StdExpansion3D(T)
78
79
80
81
82
83
        {
        }


        // Destructor
        StdPyrExp::~StdPyrExp()
mike turner's avatar
mike turner committed
84
85
        {
        }
86
87


88
89
90
        //---------------------------------------
        // Differentiation/integration Methods
        //---------------------------------------
mike turner's avatar
mike turner committed
91

92
        /**
mike turner's avatar
mike turner committed
93
94
         * \brief Calculate the derivative of the physical points
         *
95
96
         * The derivative is evaluated at the nodal physical points.
         * Derivatives with respect to the local Cartesian coordinates.
mike turner's avatar
mike turner committed
97
         *
98
99
100
101
102
103
104
105
106
107
108
109
         * \f$\begin{Bmatrix} \frac {\partial} {\partial \xi_1} \\ \frac
         * {\partial} {\partial \xi_2} \\ \frac {\partial} {\partial \xi_3}
         * \end{Bmatrix} = \begin{Bmatrix} \frac 2 {(1-\eta_3)} \frac \partial
         * {\partial \bar \eta_1} \\ \frac {\partial} {\partial \xi_2} \ \
         * \frac {(1 + \bar \eta_1)} {(1 - \eta_3)} \frac \partial {\partial
         * \bar \eta_1} + \frac {\partial} {\partial \eta_3} \end{Bmatrix}\f$
         */
        void StdPyrExp::v_PhysDeriv(
            const Array<OneD, const NekDouble> &u_physical,
                  Array<OneD,       NekDouble> &out_dxi1,
                  Array<OneD,       NekDouble> &out_dxi2,
                  Array<OneD,       NekDouble> &out_dxi3)
110
        {
111
            // PhysDerivative implementation based on Spen's book page 152.
112
113
114
115
            int    Qx = m_base[0]->GetNumPoints();
            int    Qy = m_base[1]->GetNumPoints();
            int    Qz = m_base[2]->GetNumPoints();

116
117
118
119
            Array<OneD, NekDouble> dEta_bar1(Qx*Qy*Qz,0.0);
            Array<OneD, NekDouble> dXi2     (Qx*Qy*Qz,0.0);
            Array<OneD, NekDouble> dEta3    (Qx*Qy*Qz,0.0);
            PhysTensorDeriv(u_physical, dEta_bar1, dXi2, dEta3);
120

121
            Array<OneD, const NekDouble> eta_x, eta_y, eta_z;
Peter Vos's avatar
Peter Vos committed
122
123
124
            eta_x = m_base[0]->GetZ();
            eta_y = m_base[1]->GetZ();
            eta_z = m_base[2]->GetZ();
125

126
            int i, j, k, n;
127

128
129
130
131
132
133
134
135
136
137
138
139
140
            for (k = 0, n = 0; k < Qz; ++k)
            {
                for (j = 0; j < Qy; ++j)
                {
                    for (i = 0; i < Qx; ++i, ++n)
                    {
                        if (out_dxi1.num_elements() > 0)
                            out_dxi1[n] = 2.0/(1.0 - eta_z[k]) * dEta_bar1[n];
                        if (out_dxi2.num_elements() > 0)
                            out_dxi2[n] = 2.0/(1.0 - eta_z[k]) * dXi2[n];
                        if (out_dxi3.num_elements() > 0)
                            out_dxi3[n] = (1.0+eta_x[i])/(1.0-eta_z[k])*dEta_bar1[n] +
                                (1.0+eta_y[j])/(1.0-eta_z[k])*dXi2[n] + dEta3[n];
mike turner's avatar
mike turner committed
141
                    }
142
                }
143
            }
144
145
        }

146
147
148
149
150
151
152
153
154
155
156
        void StdPyrExp::v_PhysDeriv(const int dir,
                                    const Array<OneD, const NekDouble>& inarray,
                                          Array<OneD,       NekDouble>& outarray)
        {
            switch(dir)
            {
                case 0:
                {
                    v_PhysDeriv(inarray, outarray, NullNekDouble1DArray,
                                NullNekDouble1DArray);
                    break;
157
                }
mike turner's avatar
mike turner committed
158

159
160
161
162
163
164
                case 1:
                {
                    v_PhysDeriv(inarray, NullNekDouble1DArray, outarray,
                                NullNekDouble1DArray);
                    break;
                }
mike turner's avatar
mike turner committed
165

166
167
168
169
170
171
                case 2:
                {
                    v_PhysDeriv(inarray, NullNekDouble1DArray,
                                NullNekDouble1DArray, outarray);
                    break;
                }
mike turner's avatar
mike turner committed
172

173
174
175
                default:
                {
                    ASSERTL1(false,"input dir is out of range");
176
                }
177
                break;
178
            }
179
        }
180

mike turner's avatar
mike turner committed
181
        void StdPyrExp::v_StdPhysDeriv(const Array<OneD, const NekDouble> &inarray,
182
183
184
185
                                             Array<OneD,       NekDouble> &out_d0,
                                             Array<OneD,       NekDouble> &out_d1,
                                             Array<OneD,       NekDouble> &out_d2)
        {
186
            StdPyrExp::v_PhysDeriv(inarray, out_d0, out_d1, out_d2);
187
        }
188
189
190
191
192
193
194

        void StdPyrExp::v_StdPhysDeriv(const int dir,
                                       const Array<OneD, const NekDouble>& inarray,
                                             Array<OneD,       NekDouble>& outarray)
        {
            StdPyrExp::v_PhysDeriv(dir, inarray, outarray);
        }
mike turner's avatar
mike turner committed
195

196
197
198
        //---------------------------------------
        // Transforms
        //---------------------------------------
mike turner's avatar
mike turner committed
199
200

	/**
201
202
203
204
205
         * \brief Backward transformation is evaluated at the quadrature
         * points.
         *
         * \f$ u^{\delta} (\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{m(pqr)} \hat
         * u_{pqr} \phi_{pqr} (\xi_{1i}, \xi_{2j}, \xi_{3k})\f$
mike turner's avatar
mike turner committed
206
         *
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
         * Backward transformation is three dimensional tensorial expansion
         *
         * \f$ u (\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{p=0}^{Q_x} \psi_p^a
	 *  (\xi_{1i}) \lbrace { \sum_{q=0}^{Q_y} \psi_{q}^a (\xi_{2j})
	 *  \lbrace { \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{pqr}^c (\xi_{3k})
	 *  \rbrace} \rbrace}. \f$
	 *
         * And sumfactorizing step of the form is as:\ \ \f$ f_{pqr}
         * (\xi_{3k}) = \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{pqr}^c
         * (\xi_{3k}),\\ g_{p} (\xi_{2j}, \xi_{3k}) = \sum_{r=0}^{Q_y}
         * \psi_{p}^a (\xi_{2j}) f_{pqr} (\xi_{3k}),\\ u(\xi_{1i}, \xi_{2j},
         * \xi_{3k}) = \sum_{p=0}^{Q_x} \psi_{p}^a (\xi_{1i}) g_{p}
         * (\xi_{2j}, \xi_{3k}).  \f$
         **/
        void StdPyrExp::v_BwdTrans(const Array<OneD, const NekDouble> &inarray,
                                         Array<OneD,       NekDouble> &outarray)
        {
mike turner's avatar
mike turner committed
224
            if (m_base[0]->Collocation() &&
225
226
227
228
229
230
231
232
233
234
235
236
237
                m_base[1]->Collocation() &&
                m_base[2]->Collocation())
            {
                Vmath::Vcopy(m_base[0]->GetNumPoints()*
                             m_base[1]->GetNumPoints()*
                             m_base[2]->GetNumPoints(),
                             inarray, 1, outarray, 1);
            }
            else
            {
                StdPyrExp::v_BwdTrans_SumFac(inarray,outarray);
            }
        }
238

239
240
241
242
243
244
245
        /**
         * Sum-factorisation implementation of the BwdTrans operation.
         */
        void StdPyrExp::v_BwdTrans_SumFac(
            const Array<OneD, const NekDouble>& inarray,
                  Array<OneD,       NekDouble>& outarray)
        {
246
            int  nquad0 = m_base[0]->GetNumPoints();
247
248
249
250
251
252
            int  nquad1 = m_base[1]->GetNumPoints();
            int  nquad2 = m_base[2]->GetNumPoints();
            int  order0 = m_base[0]->GetNumModes();
            int  order1 = m_base[1]->GetNumModes();

            Array<OneD, NekDouble> wsp(nquad2*order0*order1+
253
                                       nquad2*nquad1*nquad0);
254

255
256
257
258
259
260
            v_BwdTrans_SumFacKernel(m_base[0]->GetBdata(),
                                    m_base[1]->GetBdata(),
                                    m_base[2]->GetBdata(),
                                    inarray,outarray,wsp,
                                    true,true,true);
        }
261

262
263
264
265
266
267
268
269
270
271
272
        void StdPyrExp::v_BwdTrans_SumFacKernel(
            const Array<OneD, const NekDouble> &base0,
            const Array<OneD, const NekDouble> &base1,
            const Array<OneD, const NekDouble> &base2,
            const Array<OneD, const NekDouble> &inarray,
                  Array<OneD,       NekDouble> &outarray,
                  Array<OneD,       NekDouble> &wsp,
            bool                                doCheckCollDir0,
            bool                                doCheckCollDir1,
            bool                                doCheckCollDir2)
        {
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
            int  nquad0 = m_base[0]->GetNumPoints();
            int  nquad1 = m_base[1]->GetNumPoints();
            int  nquad2 = m_base[2]->GetNumPoints();

            int  order0 = m_base[0]->GetNumModes();
            int  order1 = m_base[1]->GetNumModes();
            int  order2 = m_base[2]->GetNumModes();

            Array<OneD, NekDouble > tmp  = wsp;
            Array<OneD, NekDouble > tmp1 = tmp + nquad2*order0*order1;

            int i, j, mode, mode1, cnt;

            // Perform summation over '2' direction
            mode = mode1 = cnt = 0;
            for(i = 0; i < order0; ++i)
            {
                for(j = 0; j < order1; ++j, ++cnt)
                {
                    int ijmax = max(i,j);
                    Blas::Dgemv('N', nquad2, order2-ijmax,
                                1.0, base2.get()+mode*nquad2, nquad2,
                                     inarray.get()+mode1,     1,
                                0.0, tmp.get()+cnt*nquad2,    1);
                    mode  += order2-ijmax;
                    mode1 += order2-ijmax;
                }
                //increment mode in case order1!=order2
                for(j = order1; j < order2-i; ++j)
                {
                    int ijmax = max(i,j);
                    mode += order2-ijmax;
                }
            }

            // fix for modified basis by adding split of top singular
            // vertex mode - currently (1+c)/2 x (1-b)/2 x (1-a)/2
            // component is evaluated
            if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
            {

                // Not sure why we could not use basis as 1.0 
                // top singular vertex - (1+c)/2 x (1+b)/2 x (1-a)/2 component
                Blas::Daxpy(nquad2,inarray[1],base2.get()+nquad2,1,
                            &tmp[0]+nquad2,1);

                // top singular vertex - (1+c)/2 x (1-b)/2 x (1+a)/2 component
                Blas::Daxpy(nquad2,inarray[1],base2.get()+nquad2,1,
                            &tmp[0]+order1*nquad2,1);

                // top singular vertex - (1+c)/2 x (1+b)/2 x (1+a)/2 component
                Blas::Daxpy(nquad2,inarray[1],base2.get()+nquad2,1,
                            &tmp[0]+order1*nquad2+nquad2,1);
}

            // Perform summation over '1' direction
            mode = 0;
            for(i = 0; i < order0; ++i)
            {
                Blas::Dgemm('N', 'T', nquad1, nquad2, order1,
                            1.0, base1.get(),      nquad1,
                            tmp.get()+mode*nquad2, nquad2,
                            0.0, tmp1.get()+i*nquad1*nquad2, nquad1);
                mode  += order1;
            }

            // Perform summation over '0' direction
            Blas::Dgemm('N', 'T', nquad0, nquad1*nquad2, order0,
                        1.0, base0.get(),    nquad0,
                             tmp1.get(),     nquad1*nquad2,
                        0.0, outarray.get(), nquad0);
            
345
346
        }

Sophia Han's avatar
Sophia Han committed
347
348
	/** \brief Forward transform from physical quadrature space
            stored in \a inarray and evaluate the expansion coefficients and
349
            store in \a outarray
mike turner's avatar
mike turner committed
350

Sophia Han's avatar
Sophia Han committed
351
            Inputs:\n
mike turner's avatar
mike turner committed
352

Sophia Han's avatar
Sophia Han committed
353
            - \a inarray: array of physical quadrature points to be transformed
mike turner's avatar
mike turner committed
354

Sophia Han's avatar
Sophia Han committed
355
            Outputs:\n
mike turner's avatar
mike turner committed
356
357
358
359

            - \a outarray: updated array of expansion coefficients.

        */
360
361
        void StdPyrExp::v_FwdTrans(const Array<OneD, const NekDouble> &inarray,
                                         Array<OneD,       NekDouble> &outarray)
362
        {
363
            v_IProductWRTBase(inarray,outarray);
364

365
366
367
            // get Mass matrix inverse
            StdMatrixKey      imasskey(eInvMass,DetShapeType(),*this);
            DNekMatSharedPtr  imatsys = GetStdMatrix(imasskey);
368

369
            
370
371
372
373
            // copy inarray in case inarray == outarray
            DNekVec in (m_ncoeffs, outarray);
            DNekVec out(m_ncoeffs, outarray, eWrapper);

374
            out = (*imatsys)*in;
375
        }
mike turner's avatar
mike turner committed
376
377


378
379
380
381
        //---------------------------------------
        // Inner product functions
        //---------------------------------------

mike turner's avatar
mike turner committed
382
383
384
        /** \brief  Inner product of \a inarray over region with respect to the
            expansion basis m_base[0]->GetBdata(),m_base[1]->GetBdata(), m_base[2]->GetBdata() and return in \a outarray

385
            Wrapper call to StdPyrExp::IProductWRTBase
mike turner's avatar
mike turner committed
386

387
            Input:\n
mike turner's avatar
mike turner committed
388

389
            - \a inarray: array of function evaluated at the physical collocation points
mike turner's avatar
mike turner committed
390

391
            Output:\n
mike turner's avatar
mike turner committed
392

393
            - \a outarray: array of inner product with respect to each basis over region
mike turner's avatar
mike turner committed
394

395
396
        */
        void StdPyrExp::v_IProductWRTBase(
mike turner's avatar
mike turner committed
397
            const Array<OneD, const NekDouble> &inarray,
398
                  Array<OneD,       NekDouble> &outarray)
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
        {
            if (m_base[0]->Collocation() &&
                m_base[1]->Collocation() &&
                m_base[2]->Collocation())
            {
                MultiplyByStdQuadratureMetric(inarray, outarray);
            }
            else
            {
                StdPyrExp::v_IProductWRTBase_SumFac(inarray,outarray);
            }
        }

        void StdPyrExp::v_IProductWRTBase_SumFac(
            const Array<OneD, const NekDouble>& inarray,
414
415
                  Array<OneD,       NekDouble>& outarray,
            bool                                multiplybyweights)
416
        {
417
418
419
420
421
422
423
424

            int  nquad1 = m_base[1]->GetNumPoints();
            int  nquad2 = m_base[2]->GetNumPoints();
            int  order0 = m_base[0]->GetNumModes();
            int  order1 = m_base[1]->GetNumModes();

            Array<OneD, NekDouble> wsp(nquad1*nquad2*order0 +
                                       nquad2*order0*order1);
425

426
427
428
429
430
            if(multiplybyweights)
            {
                Array<OneD, NekDouble> tmp(inarray.num_elements());

                v_MultiplyByStdQuadratureMetric(inarray, tmp);
431

432
433
434
435
436
437
438
439
440
441
442
443
444
445
                v_IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
                                               m_base[1]->GetBdata(),
                                               m_base[2]->GetBdata(),
                                               tmp,outarray,wsp,
                                               true,true,true);
            }
            else
            {
                v_IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
                                               m_base[1]->GetBdata(),
                                               m_base[2]->GetBdata(),
                                               inarray,outarray,wsp,
                                               true,true,true);
            }
446
447
448
449
450
451
452
453
454
455
456
457
        }

        void StdPyrExp::v_IProductWRTBase_SumFacKernel(
            const Array<OneD, const NekDouble> &base0,
            const Array<OneD, const NekDouble> &base1,
            const Array<OneD, const NekDouble> &base2,
            const Array<OneD, const NekDouble> &inarray,
                  Array<OneD,       NekDouble> &outarray,
                  Array<OneD,       NekDouble> &wsp,
            bool                                doCheckCollDir0,
            bool                                doCheckCollDir1,
            bool                                doCheckCollDir2)
458
        {
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
            int  nquad0 = m_base[0]->GetNumPoints();
            int  nquad1 = m_base[1]->GetNumPoints();
            int  nquad2 = m_base[2]->GetNumPoints();

            int  order0 = m_base[0]->GetNumModes();
            int  order1 = m_base[1]->GetNumModes();
            int  order2 = m_base[2]->GetNumModes();

            Array<OneD, NekDouble > tmp1 = wsp;
            Array<OneD, NekDouble > tmp2 = wsp + nquad1*nquad2*order0;

            int i,j, mode,mode1, cnt;

            // Inner product with respect to the '0' direction
            Blas::Dgemm('T', 'N', nquad1*nquad2, order0, nquad0,
                        1.0, inarray.get(), nquad0,
                             base0.get(),   nquad0,
                        0.0, tmp1.get(),    nquad1*nquad2);

            // Inner product with respect to the '1' direction
            for(mode=i=0; i < order0; ++i)
            {
                Blas::Dgemm('T', 'N', nquad2, order1, nquad1,
                            1.0, tmp1.get()+i*nquad1*nquad2, nquad1,
                                 base1.get(),    nquad1,
                            0.0, tmp2.get()+mode*nquad2,     nquad2);
                mode  += order1;
            }


            // Inner product with respect to the '2' direction
            mode = mode1 = cnt = 0;
            for(i = 0; i < order0; ++i)
            {
                for(j = 0; j < order1; ++j, ++cnt)
                {
                    int ijmax = max(i,j);

                    Blas::Dgemv('T', nquad2, order2-ijmax,
                                1.0, base2.get()+mode*nquad2, nquad2,
                                     tmp2.get()+cnt*nquad2,   1,
                                0.0, outarray.get()+mode1,    1);
                    mode  += order2-ijmax;
                    mode1 += order2-ijmax;
                }
                
                //increment mode in case order1!=order2
                for(j = order1; j < order2; ++j)
                {
                    int ijmax = max(i,j);
                    mode += order2-ijmax;
                }
            }

            // fix for modified basis for top singular vertex component
            // Already have evaluated (1+c)/2 (1-b)/2 (1-a)/2
            if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
            {
                // add in (1+c)/2 (1+b)/2 (1-a)/2  component
                outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
                                          &tmp2[nquad2],1);

                // add in (1+c)/2 (1-b)/2 (1+a)/2 component
                outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
                                          &tmp2[nquad2*order1],1);

                // add in (1+c)/2 (1+b)/2 (1+a)/2 component
                outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
                                          &tmp2[nquad2*order1+nquad2],1);
            }
529
        }
530

531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
        void StdPyrExp::v_IProductWRTDerivBase(
            const int                           dir,
            const Array<OneD, const NekDouble>& inarray,
                  Array<OneD,       NekDouble>& outarray)
        {
            StdPyrExp::v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
        }

        /**
         * @param   inarray     Function evaluated at physical collocation
         *                      points.
         * @param   outarray    Inner product with respect to each basis
         *                      function over the element.
         */
        void StdPyrExp::v_IProductWRTDerivBase_SumFac(
            const int dir,
            const Array<OneD, const NekDouble>& inarray,
                  Array<OneD,       NekDouble>& outarray)
        {
            int i;
            int nquad0  = m_base[0]->GetNumPoints();
            int nquad1  = m_base[1]->GetNumPoints();
            int nquad2  = m_base[2]->GetNumPoints();
            int nqtot   = nquad0*nquad1*nquad2;

            Array<OneD, NekDouble> gfac0(nquad0);
            Array<OneD, NekDouble> gfac1(nquad1);
            Array<OneD, NekDouble> gfac2(nquad2);
            Array<OneD, NekDouble> tmp0 (nqtot);
560
561
562
563
564
565
566


            int  order0 = m_base[0]->GetNumModes();
            int  order1 = m_base[1]->GetNumModes();

            Array<OneD, NekDouble> wsp(nquad1*nquad2*order0 +
                                       nquad2*order0*order1);
567
568
569
570
571
572
573
574
575
576
577
578
579
580

            const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
            const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
            const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();

            // set up geometric factor: (1+z0)/2
            for(i = 0; i < nquad0; ++i)
            {
                gfac0[i] = 0.5*(1+z0[i]);
            }

            // set up geometric factor: (1+z1)/2
            for(i = 0; i < nquad1; ++i)
            {
581
                gfac1[i] = 0.5*(1+z1[i]);
582
583
584
585
586
587
588
589
            }

            // Set up geometric factor: 2/(1-z2)
            for(i = 0; i < nquad2; ++i)
            {
            	gfac2[i] = 2.0/(1-z2[i]);
            }

Dave Moxey's avatar
Dave Moxey committed
590
            // Derivative in first/second direction is always scaled as follows
591
592
593
594
595
596
597
598
            const int nq01 = nquad0*nquad1;
            for(i = 0; i < nquad2; ++i)
            {
                Vmath::Smul(nq01, gfac2[i],
                            &inarray[0] + i*nq01, 1,
                            &tmp0   [0] + i*nq01, 1);
            }

Dave Moxey's avatar
Dave Moxey committed
599
            MultiplyByStdQuadratureMetric(tmp0, tmp0);
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617

            switch(dir)
            {
                case 0:
                {
                    IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
                                                 m_base[1]->GetBdata (),
                                                 m_base[2]->GetBdata (),
                                                 tmp0, outarray, wsp,
                                                 false, true, true);
                    break;
                }
                case 1:
                {
                    IProductWRTBase_SumFacKernel(m_base[0]->GetBdata (),
                                                 m_base[1]->GetDbdata(),
                                                 m_base[2]->GetBdata (),
                                                 tmp0, outarray, wsp,
Dave Moxey's avatar
Dave Moxey committed
618
                                                 true, false, true);
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
                    break;
                }
                case 2:
                {
                    Array<OneD, NekDouble> tmp3(m_ncoeffs);
                    Array<OneD, NekDouble> tmp4(m_ncoeffs);

                    // Scale eta_1 derivative by gfac0
                    for(i = 0; i < nquad1*nquad2; ++i)
                    {
                        Vmath::Vmul(nquad0, tmp0 .get() + i*nquad0, 1,
                                            gfac0.get(),            1,
                                            tmp0 .get() + i*nquad0, 1);
                    }
                    IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
                                                 m_base[1]->GetBdata(),
                                                 m_base[2]->GetBdata(),
                                                 tmp0,  tmp3, wsp,
                                                 false, true, true);

Dave Moxey's avatar
Dave Moxey committed
639
                    // Scale eta_2 derivative by gfac1*gfac2
640
641
642
643
644
645
646
647
648
649
650
651
652
                    for(i = 0; i < nquad2; ++i)
                    {
                        Vmath::Smul(nq01, gfac2[i],
                                    &inarray[0] + i*nq01, 1,
                                    &tmp0   [0] + i*nq01, 1);
                    }
                    for(i = 0; i < nquad1*nquad2; ++i)
                    {
                        Vmath::Smul(nquad0, gfac1[i%nquad1],
                                            &tmp0[0] + i*nquad0, 1,
                                            &tmp0[0] + i*nquad0, 1);
                    }

Dave Moxey's avatar
Dave Moxey committed
653
                    MultiplyByStdQuadratureMetric(tmp0, tmp0);
654
655
656
657
658
659
                    IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
                                                 m_base[1]->GetDbdata(),
                                                 m_base[2]->GetBdata(),
                                                 tmp0, tmp4,  wsp,
                                                 true, false, true);

Dave Moxey's avatar
Dave Moxey committed
660
                    MultiplyByStdQuadratureMetric(inarray,tmp0);
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
                    IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
                                                 m_base[1]->GetBdata(),
                                                 m_base[2]->GetDbdata(),
                                                 tmp0,outarray,wsp,
                                                 true, true, false);

                    Vmath::Vadd(m_ncoeffs,&tmp3[0],1,&outarray[0],1,&outarray[0],1);
                    Vmath::Vadd(m_ncoeffs,&tmp4[0],1,&outarray[0],1,&outarray[0],1);
                    break;
                }
                default:
                {
                    ASSERTL1(false, "input dir is out of range");
                    break;
                }
            }
        }

679
680
681
        //---------------------------------------
        // Evaluation functions
        //---------------------------------------
mike turner's avatar
mike turner committed
682

Dave Moxey's avatar
Dave Moxey committed
683
        void StdPyrExp::v_LocCoordToLocCollapsed(
684
            const Array<OneD, const NekDouble>& xi,
685
                  Array<OneD,       NekDouble>& eta)
686
        {
687
            if (fabs(xi[2]-1.0) < NekConstants::kNekZeroTol)
688
689
690
691
692
693
            {
                // Very top point of the pyramid
                eta[0] = -1.0;
                eta[1] = -1.0;
                eta[2] = xi[2];
            }
mike turner's avatar
mike turner committed
694
            else
695
            {
696
                // Below the line-singularity -- Common case
697
                eta[2] = xi[2]; // eta_z = xi_z
mike turner's avatar
mike turner committed
698
                eta[1] = 2.0*(1.0 + xi[1])/(1.0 - xi[2]) - 1.0;
699
                eta[0] = 2.0*(1.0 + xi[0])/(1.0 - xi[2]) - 1.0;
mike turner's avatar
mike turner committed
700
            }
701
        }
702

mike turner's avatar
mike turner committed
703
        void StdPyrExp::v_GetCoords(Array<OneD, NekDouble> &xi_x,
704
705
                                    Array<OneD, NekDouble> &xi_y,
                                    Array<OneD, NekDouble> &xi_z)
706
        {
707
            Array<OneD, const NekDouble> etaBar_x = m_base[0]->GetZ();
708
709
            Array<OneD, const NekDouble> eta_y    = m_base[1]->GetZ();
            Array<OneD, const NekDouble> eta_z    = m_base[2]->GetZ();
710
711
712
713
714
            int Qx = GetNumPoints(0);
            int Qy = GetNumPoints(1);
            int Qz = GetNumPoints(2);

            // Convert collapsed coordinates into cartesian coordinates: eta --> xi
715
716
            for (int k = 0; k < Qz; ++k )
            {
mike turner's avatar
mike turner committed
717
                for (int j = 0; j < Qy; ++j)
718
                {
mike turner's avatar
mike turner committed
719
                    for (int i = 0; i < Qx; ++i)
720
                    {
721
722
                        int s = i + Qx*(j + Qy*k);

723
                        xi_z[s] = eta_z[k];
724
                        xi_y[s] = (1.0 + eta_y[j]) * (1.0 - eta_z[k]) / 2.0  -  1.0;
725
                        xi_x[s] = (1.0 + etaBar_x[i]) * (1.0 - eta_z[k]) / 2.0  -  1.0;
726
727
728
729
730
731
732
                    }
                }
            }
        }

        void StdPyrExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
        {
733
734
735
            Array<OneD, NekDouble> tmp(m_ncoeffs, 0.0);
            tmp[mode] = 1.0;
            v_BwdTrans(tmp, outarray);
736
        }
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751

        void StdPyrExp::v_GetFaceNumModes(
                    const int                  fid,
                    const Orientation          faceOrient,
                    int &numModes0,
                    int &numModes1)
        {
            int nummodes [3] = {m_base[0]->GetNumModes(),
                                m_base[1]->GetNumModes(),
                                m_base[2]->GetNumModes()};
            switch(fid)
            {
            // quad
            case 0:
                {
752
753
                    numModes0 = nummodes[0];
                    numModes1 = nummodes[1];
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
                }
                break;
            case 1:
            case 3:
                {
                    numModes0 = nummodes[0];
                    numModes1 = nummodes[2];
                }
                break;
            case 2:
            case 4:
                {
                    numModes0 = nummodes[1];
                    numModes1 = nummodes[2];
                }
                break;
            }
771
772
773
774
775

            if ( faceOrient >= 9 )
            {
                std::swap(numModes0, numModes1);
            }
776
777
        }

778
779
780
781
782
783
784
785
786
787
788
789
790
        //---------------------------------------
        // Helper functions
        //---------------------------------------

        int StdPyrExp::v_GetNverts() const
        {
            return 5;
        }

        int StdPyrExp::v_GetNedges() const
        {
            return 8;
        }
791

792
        int StdPyrExp::v_GetNfaces() const
793
        {
794
795
            return 5;
        }
mike turner's avatar
mike turner committed
796

797
        LibUtilities::ShapeType StdPyrExp::v_DetShapeType() const
798
        {
799
            return LibUtilities::ePyramid;
800
801
802
803
804
805
806
807
808
809
        }

        int StdPyrExp::v_NumBndryCoeffs() const
        {
            ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
                     GetBasisType(0) == LibUtilities::eGLL_Lagrange,
                     "BasisType is not a boundary interior form");
            ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
                     GetBasisType(1) == LibUtilities::eGLL_Lagrange,
                     "BasisType is not a boundary interior form");
810
            ASSERTL1(GetBasisType(2) == LibUtilities::eModifiedPyr_C ||
811
812
                     GetBasisType(2) == LibUtilities::eGLL_Lagrange,
                     "BasisType is not a boundary interior form");
mike turner's avatar
mike turner committed
813

814
815
816
            int P = m_base[0]->GetNumModes();
            int Q = m_base[1]->GetNumModes();
            int R = m_base[2]->GetNumModes();
mike turner's avatar
mike turner committed
817

818
819
            return LibUtilities::StdPyrData::
                                    getNumberOfBndCoefficients(P, Q, R);
820
821
822
823
824
        }

        int StdPyrExp::v_GetEdgeNcoeffs(const int i) const
        {
            ASSERTL2(i >= 0 && i <= 7, "edge id is out of range");
mike turner's avatar
mike turner committed
825

826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
            if (i == 0 || i == 2)
            {
                return GetBasisNumModes(0);
            }
            else if (i == 1 || i == 3)
            {
                return GetBasisNumModes(1);
            }
            else
            {
                return GetBasisNumModes(2);
            }
        }

        int StdPyrExp::v_GetFaceNcoeffs(const int i) const
        {
            ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
mike turner's avatar
mike turner committed
843

844
845
846
847
848
849
            if (i == 0)
            {
                return GetBasisNumModes(0)*GetBasisNumModes(1);
            }
            else if (i == 1 || i == 3)
            {
850
851
                int P = GetBasisNumModes(0)-1, Q = GetBasisNumModes(2)-1;
                return Q+1 + (P*(1 + 2*Q - P))/2;
852
853
854
            }
            else
            {
855
856
                int P = GetBasisNumModes(1)-1, Q = GetBasisNumModes(2)-1;
                return Q+1 + (P*(1 + 2*Q - P))/2;
857
858
            }
        }
mike turner's avatar
mike turner committed
859

860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
        int StdPyrExp::v_GetFaceIntNcoeffs(const int i) const
        {
            ASSERTL2(i >= 0 && i <= 4, "face id is out of range");

            int P = m_base[0]->GetNumModes()-1;
            int Q = m_base[1]->GetNumModes()-1;
            int R = m_base[2]->GetNumModes()-1;

            if (i == 0)
            {
                return (P-1)*(Q-1);
            }
            else if (i == 1 || i == 3)
            {
                return (P-1) * (2*(R-1) - (P-1) - 1) / 2;
            }
            else
            {
                return (Q-1) * (2*(R-1) - (Q-1) - 1) / 2;
            }
        }

882
883
884
        int StdPyrExp::v_GetFaceNumPoints(const int i) const
        {
            ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
mike turner's avatar
mike turner committed
885

886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
            if (i == 0)
            {
                return m_base[0]->GetNumPoints()*
                       m_base[1]->GetNumPoints();
            }
            else if (i == 1 || i == 3)
            {
                return m_base[0]->GetNumPoints()*
                       m_base[2]->GetNumPoints();
            }
            else
            {
                return m_base[1]->GetNumPoints()*
                       m_base[2]->GetNumPoints();
            }
        }


904
905
906
907
908
909
910
911
912
913
914
915
916
917
        const LibUtilities::BasisKey StdPyrExp::v_DetFaceBasisKey(
            const int i, const int k) const
        {
            ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
            ASSERTL2(k >= 0 && k <= 1, "basis key id is out of range");

            switch(i)
            {
                case 0:
                {
                    return EvaluateQuadFaceBasisKey(k,
                                                    m_base[k]->GetBasisType(),
                                                    m_base[k]->GetNumPoints(),
                                                    m_base[k]->GetNumModes());
mike turner's avatar
mike turner committed
918

919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
                }
                case 1:
                case 3:
                {
                    return EvaluateTriFaceBasisKey(k,
                                                   m_base[2*k]->GetBasisType(),
                                                   m_base[2*k]->GetNumPoints(),
                                                   m_base[2*k]->GetNumModes());
                }
                case 2:
                case 4:
                {
                    return EvaluateTriFaceBasisKey(k,
                                                   m_base[k+1]->GetBasisType(),
                                                   m_base[k+1]->GetNumPoints(),
                                                   m_base[k+1]->GetNumModes());
                }
            }

            // Should never get here.
            return LibUtilities::NullBasisKey;
        }

942
        int StdPyrExp::v_CalcNumberOfCoefficients(
mike turner's avatar
mike turner committed
943
            const std::vector<unsigned int> &nummodes,
944
945
            int &modes_offset)
        {
946
            int nmodes = LibUtilities::StdPyrData::getNumberOfCoefficients(
947
948
949
                nummodes[modes_offset],
                nummodes[modes_offset+1],
                nummodes[modes_offset+2]);
mike turner's avatar
mike turner committed
950

951
952
953
            modes_offset += 3;
            return nmodes;
        }
mike turner's avatar
mike turner committed
954

955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
        LibUtilities::BasisType StdPyrExp::v_GetEdgeBasisType(const int i) const
        {
            ASSERTL2(i >= 0 && i <= 7, "edge id is out of range");
            if (i == 0 || i == 2)
            {
                return GetBasisType(0);
            }
            else if (i == 1 || i == 3)
            {
                return GetBasisType(1);
            }
            else
            {
                return GetBasisType(2);
            }
        }


        //---------------------------------------
        // Mappings
        //---------------------------------------
mike turner's avatar
mike turner committed
976

977
        void StdPyrExp::v_GetFaceToElementMap(
mike turner's avatar
mike turner committed
978
            const int                  fid,
Dave Moxey's avatar
Dave Moxey committed
979
            const Orientation          faceOrient,
980
981
            Array<OneD, unsigned int> &maparray,
            Array<OneD,          int> &signarray,
mike turner's avatar
mike turner committed
982
            int                        P,
983
            int                        Q)
984
        {
985
986
987
            ASSERTL1(GetEdgeBasisType(0) == GetEdgeBasisType(1),
                     "Method only implemented if BasisType is identical"
                     "in x and y directions");
mike turner's avatar
mike turner committed
988
            ASSERTL1(GetEdgeBasisType(0) == LibUtilities::eModified_A &&
989
                     GetEdgeBasisType(4) == LibUtilities::eModifiedPyr_C,
990
                     "Method only implemented for Modified_A BasisType"
991
                     "(x and y direction) and ModifiedPyr_C BasisType (z "
992
993
                     "direction)");

994
            int i, j, k, p, q, r, nFaceCoeffs, idx = 0;
mike turner's avatar
mike turner committed
995
996
            int nummodesA=0, nummodesB=0;

997
998
999
1000
            int order0 = m_base[0]->GetNumModes();
            int order1 = m_base[1]->GetNumModes();
            int order2 = m_base[2]->GetNumModes();

1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
            switch (fid)
            {
            case 0:
                nummodesA = order0;
                nummodesB = order1;
                break;
            case 1:
            case 3:
                nummodesA = order0;
                nummodesB = order2;
                break;
            case 2:
            case 4:
                nummodesA = order1;
                nummodesB = order2;
                break;
            }
mike turner's avatar
mike turner committed
1018

1019
1020
1021
            bool CheckForZeroedModes = false;

            if (P == -1)
1022
            {
1023
1024
                P = nummodesA;
                Q = nummodesB;
1025
1026
1027
1028
                nFaceCoeffs = GetFaceNcoeffs(fid);
            }
            else if (fid > 0)
            {
1029
1030
                nFaceCoeffs = P*(2*Q-P+1)/2;
                CheckForZeroedModes = true;
1031
1032
1033
            }
            else
            {
1034
1035
                nFaceCoeffs = P*Q;
                CheckForZeroedModes = true;
1036
1037
1038
1039
1040
1041
1042
            }

            // Allocate the map array and sign array; set sign array to ones (+)
            if (maparray.num_elements() != nFaceCoeffs)
            {
                maparray = Array<OneD, unsigned int>(nFaceCoeffs);
            }
mike turner's avatar
mike turner committed
1043

1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
            if (signarray.num_elements() != nFaceCoeffs)
            {
                signarray = Array<OneD, int>(nFaceCoeffs,1);
            }
            else
            {
                fill(signarray.get(), signarray.get() + nFaceCoeffs, 1);
            }

            // Set up an array indexing for quads, since the ordering may need
            // to be transposed.
            Array<OneD, int> arrayindx(nFaceCoeffs,-1);
1056

1057
1058
            if (fid == 0)
            {
1059
                for (i = 0; i < Q; i++)
1060
                {
1061
                    for (j = 0; j < P; j++)
1062
1063
1064
                    {
                        if (faceOrient < 9)
                        {
1065
                            arrayindx[i*P+j] = i*P+j;
1066
1067
1068
                        }
                        else
                        {
1069
                            arrayindx[i*P+j] = j*Q+i;
1070
1071
1072
1073
1074
1075
1076
                        }
                    }
                }
            }

            // Set up ordering inside each 2D face. Also for triangular faces,
            // populate signarray.
mike turner's avatar
mike turner committed
1077
            switch (fid)
1078
            {
1079
            case 0: // Bottom quad
mike turner's avatar
mike turner committed
1080

1081
1082
1083
                for (q = 0; q < Q; ++q)
                {
                    for (p = 0; p < P; ++p)
1084
                    {
1085
                        maparray[arrayindx[q*P+p]] = GetMode(p,q,0);
1086
                    }
1087
1088
1089
1090
1091
1092
1093
                }
                break;
                
            case 1: // Front triangle
                for (p = 0; p < P; ++p)
                {
                    for (r = 0; r < Q-p; ++r)
1094
                    {
1095
                        if ((int)faceOrient == 7 && p > 1)
1096
                        {
1097
                            signarray[idx] = p % 2 ? -1 : 1;
1098
                        }
1099
                        maparray[idx++] = GetMode(p,0,r);
1100
                    }
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
                }
                break;
                
            case 2: // Right triangle
                maparray[idx++] = GetMode(1,0,0);
                maparray[idx++] = GetMode(0,0,1);
                for (r = 1; r < Q-1; ++r)
                {
                    maparray[idx++] = GetMode(1,0,r);
                }
1111

1112
1113
1114
                for (q = 1; q < P; ++q)
                {
                    for (r = 0; r < Q-q; ++r)
1115
                    {
1116
                        if ((int)faceOrient == 7 && q > 1)
1117
                        {
1118
                            signarray[idx] = q % 2 ? -1 : 1;
1119
                        }
1120
                        maparray[idx++] = GetMode(1,q,r);
1121
                    }
1122
1123
                }
                break;
1124

1125
1126
1127
1128
1129
1130
1131
            case 3: // Rear triangle
                maparray[idx++] = GetMode(0,1,0);
                maparray[idx++] = GetMode(0,0,1);
                for (r = 1; r < Q-1; ++r)
                {
                    maparray[idx++] = GetMode(0,1,r);
                }
mike turner's avatar
mike turner committed
1132

1133
1134
1135
                for (p = 1; p < P; ++p)
                {
                    for (r = 0; r < Q-p; ++r)
1136
                    {
1137
                        if ((int)faceOrient == 7 && p > 1)
1138
                        {
1139
                            signarray[idx] = p % 2 ? -1 : 1;
1140
                        }
1141
                        maparray[idx++] = GetMode(p, 1, r);
1142
                    }
1143
1144
1145
1146
1147
1148
1149
                }
                break;
                
            case 4: // Left triangle
                for (q = 0; q < P; ++q)
                {
                    for (r = 0; r < Q-q; ++r)
1150
                    {
1151
                        if ((int)faceOrient == 7 && q > 1)
1152
                        {
1153
                            signarray[idx] = q % 2 ? -1 : 1;
1154
                        }
1155
                        maparray[idx++] = GetMode(0,q,r);
1156
                    }
1157
1158
                }
                break;
mike turner's avatar
mike turner committed
1159

1160
1161
            default:
                ASSERTL0(false, "Face to element map unavailable.");
1162
            }
1163

1164
1165
            if (fid > 0)
            {
1166
1167
1168
1169
               if(CheckForZeroedModes)
                {
                    // zero signmap and set maparray to zero if elemental
                    // modes are not as large as face modesl
mike turner's avatar
mike turner committed
1170
                    int idx = 0;
1171
                    for (j = 0; j < P; ++j)
1172
                    {
1173
1174
                        idx += Q-j;
                        for (k = Q-j; k < Q-j; ++k)
1175
1176
1177
1178
1179
                        {
                            signarray[idx]  = 0.0;
                            maparray[idx++] = maparray[0];
                        }
                    }
mike turner's avatar
mike turner committed
1180

1181
                    for (j = P; j < P; ++j)
1182
1183
1184
1185
1186
1187
1188
1189
                    {
                        for (k = 0; k < Q-j; ++k)
                        {
                            signarray[idx]  = 0.0;
                            maparray[idx++] = maparray[0];
                        }
                    }
                }
mike turner's avatar
mike turner committed
1190

1191
1192
1193
1194
                // Triangles only have one possible orientation (base
                // direction reversed); swap edge modes.
                if ((int)faceOrient == 7)
                {
1195
1196
                    swap(maparray[0], maparray[Q]);
                    for (i = 1; i < Q-1; ++i)
1197
                    {
1198
                        swap(maparray[i+1], maparray[Q+i]);
1199
1200
1201
1202
1203
                    }
                }
            }
            else
            {
1204
1205
1206
1207
                if(CheckForZeroedModes)
                {
                    // zero signmap and set maparray to zero if elemental
                    // modes are not as large as face modesl
1208
                    for (j = 0; j < P; ++j)
1209
                    {
1210
                        for (k = Q; k < Q; ++k)
1211
1212
1213
1214
1215
1216
                        {
                            signarray[arrayindx[j+k*P]] = 0.0;
                            maparray[arrayindx[j+k*P]]  = maparray[0];
                        }
                    }

1217
                    for (j = P; j < P; ++j)
1218
1219
1220
1221
1222
1223
                    {
                        for (k = 0; k < Q; ++k)
                        {
                            signarray[arrayindx[j+k*P]] = 0.0;
                            maparray[arrayindx[j+k*P]]  = maparray[0];
                        }
mike turner's avatar
mike turner committed
1224
                    }
1225
1226
                }

1227
1228
1229
1230
1231
1232
1233
1234
1235
                // The code below is exactly the same as that taken from
                // StdHexExp and reverses the 'b' and 'a' directions as
                // appropriate (1st and 2nd if statements respectively) in
                // quadrilateral faces.
                if (faceOrient == 6 || faceOrient == 8 ||
                    faceOrient == 11 || faceOrient == 12)
                {
                    if (faceOrient < 9)
                    {
1236
                        for (i = 3; i < Q; i += 2)
1237
                        {