Nektar++
StdPyrExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdPyrExp.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: pyramidic routines built upon StdExpansion3D
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
37 #include <StdRegions/StdPyrExp.h>
38 #include <LibUtilities/Foundations/ManagerAccess.h>
39 #include <iomanip>
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45  namespace StdRegions
46  {
47  StdPyrExp::StdPyrExp() // Deafult construct of standard expansion directly called.
48  {
49  }
50 
51  StdPyrExp::StdPyrExp(const LibUtilities::BasisKey &Ba,
52  const LibUtilities::BasisKey &Bb,
53  const LibUtilities::BasisKey &Bc)
54  : StdExpansion (LibUtilities::StdPyrData::getNumberOfCoefficients(
55  Ba.GetNumModes(),
56  Bb.GetNumModes(),
57  Bc.GetNumModes()),
58  3, Ba, Bb, Bc),
59  StdExpansion3D(LibUtilities::StdPyrData::getNumberOfCoefficients(
60  Ba.GetNumModes(),
61  Bb.GetNumModes(),
62  Bc.GetNumModes()),
63  Ba, Bb, Bc)
64  {
65 
66  ASSERTL0(Ba.GetNumModes() <= Bc.GetNumModes(), "order in 'a' direction is higher "
67  "than order in 'c' direction");
68  ASSERTL0(Bb.GetNumModes() <= Bc.GetNumModes(), "order in 'b' direction is higher "
69  "than order in 'c' direction");
70  ASSERTL1(Bc.GetBasisType() == LibUtilities::eModifiedPyr_C ||
71  Bc.GetBasisType() == LibUtilities::eOrthoPyr_C,
72  "Expected basis type in 'c' direction to be ModifiedPyr_C or OrthoPyr_C");
73 
74  }
75 
76  StdPyrExp::StdPyrExp(const StdPyrExp &T)
77  : StdExpansion (T),
78  StdExpansion3D(T)
79  {
80  }
81 
82 
83  // Destructor
84  StdPyrExp::~StdPyrExp()
85  {
86  }
87 
88 
89  //---------------------------------------
90  // Differentiation/integration Methods
91  //---------------------------------------
92 
93  /**
94  * \brief Calculate the derivative of the physical points
95  *
96  * The derivative is evaluated at the nodal physical points.
97  * Derivatives with respect to the local Cartesian coordinates.
98  *
99  * \f$\begin{Bmatrix} \frac {\partial} {\partial \xi_1} \\ \frac
100  * {\partial} {\partial \xi_2} \\ \frac {\partial} {\partial \xi_3}
101  * \end{Bmatrix} = \begin{Bmatrix} \frac 2 {(1-\eta_3)} \frac \partial
102  * {\partial \bar \eta_1} \\ \frac {\partial} {\partial \xi_2} \ \
103  * \frac {(1 + \bar \eta_1)} {(1 - \eta_3)} \frac \partial {\partial
104  * \bar \eta_1} + \frac {\partial} {\partial \eta_3} \end{Bmatrix}\f$
105  */
106  void StdPyrExp::v_PhysDeriv(
107  const Array<OneD, const NekDouble> &u_physical,
108  Array<OneD, NekDouble> &out_dxi1,
109  Array<OneD, NekDouble> &out_dxi2,
110  Array<OneD, NekDouble> &out_dxi3)
111  {
112  // PhysDerivative implementation based on Spen's book page 152.
113  int Qx = m_base[0]->GetNumPoints();
114  int Qy = m_base[1]->GetNumPoints();
115  int Qz = m_base[2]->GetNumPoints();
116 
117  Array<OneD, NekDouble> dEta_bar1(Qx*Qy*Qz,0.0);
118  Array<OneD, NekDouble> dXi2 (Qx*Qy*Qz,0.0);
119  Array<OneD, NekDouble> dEta3 (Qx*Qy*Qz,0.0);
120  PhysTensorDeriv(u_physical, dEta_bar1, dXi2, dEta3);
121 
122  Array<OneD, const NekDouble> eta_x, eta_y, eta_z;
123  eta_x = m_base[0]->GetZ();
124  eta_y = m_base[1]->GetZ();
125  eta_z = m_base[2]->GetZ();
126 
127  int i, j, k, n;
128 
129  if (out_dxi1.num_elements() > 0)
130  {
131  for (k = 0, n = 0; k < Qz; ++k)
132  {
133  NekDouble fac = 2.0/(1.0 - eta_z[k]);
134  for (j = 0; j < Qy; ++j)
135  {
136  for (i = 0; i < Qx; ++i, ++n)
137  {
138  out_dxi1[n] = fac * dEta_bar1[n];
139  }
140  }
141  }
142  }
143 
144  if (out_dxi2.num_elements() > 0)
145  {
146  for (k = 0, n = 0; k < Qz; ++k)
147  {
148  NekDouble fac = 2.0/(1.0 - eta_z[k]);
149  for (j = 0; j < Qy; ++j)
150  {
151  for (i = 0; i < Qx; ++i, ++n)
152  {
153  out_dxi2[n] = fac * dXi2[n];
154  }
155  }
156  }
157  }
158 
159  if (out_dxi3.num_elements() > 0)
160  {
161  for (k = 0, n = 0; k < Qz; ++k)
162  {
163  NekDouble fac = 1.0/(1.0 - eta_z[k]);
164  for (j = 0; j < Qy; ++j)
165  {
166  NekDouble fac1 = (1.0+eta_y[j]);
167  for (i = 0; i < Qx; ++i, ++n)
168  {
169  out_dxi3[n] = (1.0+eta_x[i])*fac*dEta_bar1[n] +
170  fac1*fac*dXi2[n] + dEta3[n];
171  }
172  }
173  }
174  }
175  }
176 
177  void StdPyrExp::v_PhysDeriv(const int dir,
178  const Array<OneD, const NekDouble>& inarray,
179  Array<OneD, NekDouble>& outarray)
180  {
181  switch(dir)
182  {
183  case 0:
184  {
185  v_PhysDeriv(inarray, outarray, NullNekDouble1DArray,
186  NullNekDouble1DArray);
187  break;
188  }
189 
190  case 1:
191  {
192  v_PhysDeriv(inarray, NullNekDouble1DArray, outarray,
193  NullNekDouble1DArray);
194  break;
195  }
196 
197  case 2:
198  {
199  v_PhysDeriv(inarray, NullNekDouble1DArray,
200  NullNekDouble1DArray, outarray);
201  break;
202  }
203 
204  default:
205  {
206  ASSERTL1(false,"input dir is out of range");
207  }
208  break;
209  }
210  }
211 
212  void StdPyrExp::v_StdPhysDeriv(const Array<OneD, const NekDouble> &inarray,
213  Array<OneD, NekDouble> &out_d0,
214  Array<OneD, NekDouble> &out_d1,
215  Array<OneD, NekDouble> &out_d2)
216  {
217  StdPyrExp::v_PhysDeriv(inarray, out_d0, out_d1, out_d2);
218  }
219 
220  void StdPyrExp::v_StdPhysDeriv(const int dir,
221  const Array<OneD, const NekDouble>& inarray,
222  Array<OneD, NekDouble>& outarray)
223  {
224  StdPyrExp::v_PhysDeriv(dir, inarray, outarray);
225  }
226 
227  //---------------------------------------
228  // Transforms
229  //---------------------------------------
230 
231  /**
232  * \brief Backward transformation is evaluated at the quadrature
233  * points.
234  *
235  * \f$ u^{\delta} (\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{m(pqr)} \hat
236  * u_{pqr} \phi_{pqr} (\xi_{1i}, \xi_{2j}, \xi_{3k})\f$
237  *
238  * Backward transformation is three dimensional tensorial expansion
239  *
240  * \f$ u (\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{p=0}^{Q_x} \psi_p^a
241  * (\xi_{1i}) \lbrace { \sum_{q=0}^{Q_y} \psi_{q}^a (\xi_{2j})
242  * \lbrace { \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{pqr}^c (\xi_{3k})
243  * \rbrace} \rbrace}. \f$
244  *
245  * And sumfactorizing step of the form is as:\ \ \f$ f_{pqr}
246  * (\xi_{3k}) = \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{pqr}^c
247  * (\xi_{3k}),\\ g_{p} (\xi_{2j}, \xi_{3k}) = \sum_{r=0}^{Q_y}
248  * \psi_{p}^a (\xi_{2j}) f_{pqr} (\xi_{3k}),\\ u(\xi_{1i}, \xi_{2j},
249  * \xi_{3k}) = \sum_{p=0}^{Q_x} \psi_{p}^a (\xi_{1i}) g_{p}
250  * (\xi_{2j}, \xi_{3k}). \f$
251  **/
252  void StdPyrExp::v_BwdTrans(const Array<OneD, const NekDouble> &inarray,
253  Array<OneD, NekDouble> &outarray)
254  {
255  if (m_base[0]->Collocation() &&
256  m_base[1]->Collocation() &&
257  m_base[2]->Collocation())
258  {
259  Vmath::Vcopy(m_base[0]->GetNumPoints()*
260  m_base[1]->GetNumPoints()*
261  m_base[2]->GetNumPoints(),
262  inarray, 1, outarray, 1);
263  }
264  else
265  {
266  StdPyrExp::v_BwdTrans_SumFac(inarray,outarray);
267  }
268  }
269 
270  /**
271  * Sum-factorisation implementation of the BwdTrans operation.
272  */
273  void StdPyrExp::v_BwdTrans_SumFac(
274  const Array<OneD, const NekDouble>& inarray,
275  Array<OneD, NekDouble>& outarray)
276  {
277  int nquad0 = m_base[0]->GetNumPoints();
278  int nquad1 = m_base[1]->GetNumPoints();
279  int nquad2 = m_base[2]->GetNumPoints();
280  int order0 = m_base[0]->GetNumModes();
281  int order1 = m_base[1]->GetNumModes();
282 
283  Array<OneD, NekDouble> wsp(nquad2*order0*order1+
284  nquad2*nquad1*nquad0);
285 
286  v_BwdTrans_SumFacKernel(m_base[0]->GetBdata(),
287  m_base[1]->GetBdata(),
288  m_base[2]->GetBdata(),
289  inarray,outarray,wsp,
290  true,true,true);
291  }
292 
293  void StdPyrExp::v_BwdTrans_SumFacKernel(
294  const Array<OneD, const NekDouble> &base0,
295  const Array<OneD, const NekDouble> &base1,
296  const Array<OneD, const NekDouble> &base2,
297  const Array<OneD, const NekDouble> &inarray,
298  Array<OneD, NekDouble> &outarray,
299  Array<OneD, NekDouble> &wsp,
300  bool doCheckCollDir0,
301  bool doCheckCollDir1,
302  bool doCheckCollDir2)
303  {
304  boost::ignore_unused(doCheckCollDir0, doCheckCollDir1,
305  doCheckCollDir2);
306 
307  int nquad0 = m_base[0]->GetNumPoints();
308  int nquad1 = m_base[1]->GetNumPoints();
309  int nquad2 = m_base[2]->GetNumPoints();
310 
311  int order0 = m_base[0]->GetNumModes();
312  int order1 = m_base[1]->GetNumModes();
313  int order2 = m_base[2]->GetNumModes();
314 
315  Array<OneD, NekDouble > tmp = wsp;
316  Array<OneD, NekDouble > tmp1 = tmp + nquad2*order0*order1;
317 
318  int i, j, mode, mode1, cnt;
319 
320  // Perform summation over '2' direction
321  mode = mode1 = cnt = 0;
322  for(i = 0; i < order0; ++i)
323  {
324  for(j = 0; j < order1; ++j, ++cnt)
325  {
326  int ijmax = max(i,j);
327  Blas::Dgemv('N', nquad2, order2-ijmax,
328  1.0, base2.get()+mode*nquad2, nquad2,
329  inarray.get()+mode1, 1,
330  0.0, tmp.get()+cnt*nquad2, 1);
331  mode += order2-ijmax;
332  mode1 += order2-ijmax;
333  }
334  //increment mode in case order1!=order2
335  for(j = order1; j < order2-i; ++j)
336  {
337  int ijmax = max(i,j);
338  mode += order2-ijmax;
339  }
340  }
341 
342  // fix for modified basis by adding split of top singular
343  // vertex mode - currently (1+c)/2 x (1-b)/2 x (1-a)/2
344  // component is evaluated
345  if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
346  {
347 
348  // Not sure why we could not use basis as 1.0
349  // top singular vertex - (1+c)/2 x (1+b)/2 x (1-a)/2 component
350  Blas::Daxpy(nquad2,inarray[1],base2.get()+nquad2,1,
351  &tmp[0]+nquad2,1);
352 
353  // top singular vertex - (1+c)/2 x (1-b)/2 x (1+a)/2 component
354  Blas::Daxpy(nquad2,inarray[1],base2.get()+nquad2,1,
355  &tmp[0]+order1*nquad2,1);
356 
357  // top singular vertex - (1+c)/2 x (1+b)/2 x (1+a)/2 component
358  Blas::Daxpy(nquad2,inarray[1],base2.get()+nquad2,1,
359  &tmp[0]+order1*nquad2+nquad2,1);
360 }
361 
362  // Perform summation over '1' direction
363  mode = 0;
364  for(i = 0; i < order0; ++i)
365  {
366  Blas::Dgemm('N', 'T', nquad1, nquad2, order1,
367  1.0, base1.get(), nquad1,
368  tmp.get()+mode*nquad2, nquad2,
369  0.0, tmp1.get()+i*nquad1*nquad2, nquad1);
370  mode += order1;
371  }
372 
373  // Perform summation over '0' direction
374  Blas::Dgemm('N', 'T', nquad0, nquad1*nquad2, order0,
375  1.0, base0.get(), nquad0,
376  tmp1.get(), nquad1*nquad2,
377  0.0, outarray.get(), nquad0);
378 
379  }
380 
381  /** \brief Forward transform from physical quadrature space
382  stored in \a inarray and evaluate the expansion coefficients and
383  store in \a outarray
384 
385  Inputs:\n
386 
387  - \a inarray: array of physical quadrature points to be transformed
388 
389  Outputs:\n
390 
391  - \a outarray: updated array of expansion coefficients.
392 
393  */
394  void StdPyrExp::v_FwdTrans(const Array<OneD, const NekDouble> &inarray,
395  Array<OneD, NekDouble> &outarray)
396  {
397  v_IProductWRTBase(inarray,outarray);
398 
399  // get Mass matrix inverse
400  StdMatrixKey imasskey(eInvMass,DetShapeType(),*this);
401  DNekMatSharedPtr imatsys = GetStdMatrix(imasskey);
402 
403 
404  // copy inarray in case inarray == outarray
405  DNekVec in (m_ncoeffs, outarray);
406  DNekVec out(m_ncoeffs, outarray, eWrapper);
407 
408  out = (*imatsys)*in;
409  }
410 
411 
412  //---------------------------------------
413  // Inner product functions
414  //---------------------------------------
415 
416  /** \brief Inner product of \a inarray over region with respect to the
417  expansion basis m_base[0]->GetBdata(),m_base[1]->GetBdata(), m_base[2]->GetBdata() and return in \a outarray
418 
419  Wrapper call to StdPyrExp::IProductWRTBase
420 
421  Input:\n
422 
423  - \a inarray: array of function evaluated at the physical collocation points
424 
425  Output:\n
426 
427  - \a outarray: array of inner product with respect to each basis over region
428 
429  */
430  void StdPyrExp::v_IProductWRTBase(
431  const Array<OneD, const NekDouble> &inarray,
432  Array<OneD, NekDouble> &outarray)
433  {
434  if (m_base[0]->Collocation() &&
435  m_base[1]->Collocation() &&
436  m_base[2]->Collocation())
437  {
438  MultiplyByStdQuadratureMetric(inarray, outarray);
439  }
440  else
441  {
442  StdPyrExp::v_IProductWRTBase_SumFac(inarray,outarray);
443  }
444  }
445 
446  void StdPyrExp::v_IProductWRTBase_SumFac(
447  const Array<OneD, const NekDouble>& inarray,
448  Array<OneD, NekDouble>& outarray,
449  bool multiplybyweights)
450  {
451 
452  int nquad1 = m_base[1]->GetNumPoints();
453  int nquad2 = m_base[2]->GetNumPoints();
454  int order0 = m_base[0]->GetNumModes();
455  int order1 = m_base[1]->GetNumModes();
456 
457  Array<OneD, NekDouble> wsp(order0*nquad2*(nquad1 + order1));
458 
459  if(multiplybyweights)
460  {
461  Array<OneD, NekDouble> tmp(inarray.num_elements());
462 
463  v_MultiplyByStdQuadratureMetric(inarray, tmp);
464 
465  v_IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
466  m_base[1]->GetBdata(),
467  m_base[2]->GetBdata(),
468  tmp,outarray,wsp,
469  true,true,true);
470  }
471  else
472  {
473  v_IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
474  m_base[1]->GetBdata(),
475  m_base[2]->GetBdata(),
476  inarray,outarray,wsp,
477  true,true,true);
478  }
479  }
480 
481  void StdPyrExp::v_IProductWRTBase_SumFacKernel(
482  const Array<OneD, const NekDouble> &base0,
483  const Array<OneD, const NekDouble> &base1,
484  const Array<OneD, const NekDouble> &base2,
485  const Array<OneD, const NekDouble> &inarray,
486  Array<OneD, NekDouble> &outarray,
487  Array<OneD, NekDouble> &wsp,
488  bool doCheckCollDir0,
489  bool doCheckCollDir1,
490  bool doCheckCollDir2)
491  {
492  boost::ignore_unused(doCheckCollDir0, doCheckCollDir1,
493  doCheckCollDir2);
494 
495  int nquad0 = m_base[0]->GetNumPoints();
496  int nquad1 = m_base[1]->GetNumPoints();
497  int nquad2 = m_base[2]->GetNumPoints();
498 
499  int order0 = m_base[0]->GetNumModes();
500  int order1 = m_base[1]->GetNumModes();
501  int order2 = m_base[2]->GetNumModes();
502 
503  ASSERTL1(wsp.num_elements() >= nquad1*nquad2*order0 +
504  nquad2*order0*order1,
505  "Insufficient workspace size");
506 
507  Array<OneD, NekDouble > tmp1 = wsp;
508  Array<OneD, NekDouble > tmp2 = wsp + nquad1*nquad2*order0;
509 
510  int i,j, mode,mode1, cnt;
511 
512  // Inner product with respect to the '0' direction
513  Blas::Dgemm('T', 'N', nquad1*nquad2, order0, nquad0,
514  1.0, inarray.get(), nquad0,
515  base0.get(), nquad0,
516  0.0, tmp1.get(), nquad1*nquad2);
517 
518  // Inner product with respect to the '1' direction
519  for(mode=i=0; i < order0; ++i)
520  {
521  Blas::Dgemm('T', 'N', nquad2, order1, nquad1,
522  1.0, tmp1.get()+i*nquad1*nquad2, nquad1,
523  base1.get(), nquad1,
524  0.0, tmp2.get()+mode*nquad2, nquad2);
525  mode += order1;
526  }
527 
528 
529  // Inner product with respect to the '2' direction
530  mode = mode1 = cnt = 0;
531  for(i = 0; i < order0; ++i)
532  {
533  for(j = 0; j < order1; ++j, ++cnt)
534  {
535  int ijmax = max(i,j);
536 
537  Blas::Dgemv('T', nquad2, order2-ijmax,
538  1.0, base2.get()+mode*nquad2, nquad2,
539  tmp2.get()+cnt*nquad2, 1,
540  0.0, outarray.get()+mode1, 1);
541  mode += order2-ijmax;
542  mode1 += order2-ijmax;
543  }
544 
545  //increment mode in case order1!=order2
546  for(j = order1; j < order2; ++j)
547  {
548  int ijmax = max(i,j);
549  mode += order2-ijmax;
550  }
551  }
552 
553  // fix for modified basis for top singular vertex component
554  // Already have evaluated (1+c)/2 (1-b)/2 (1-a)/2
555  if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
556  {
557  // add in (1+c)/2 (1+b)/2 (1-a)/2 component
558  outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
559  &tmp2[nquad2],1);
560 
561  // add in (1+c)/2 (1-b)/2 (1+a)/2 component
562  outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
563  &tmp2[nquad2*order1],1);
564 
565  // add in (1+c)/2 (1+b)/2 (1+a)/2 component
566  outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
567  &tmp2[nquad2*order1+nquad2],1);
568  }
569  }
570 
571  void StdPyrExp::v_IProductWRTDerivBase(
572  const int dir,
573  const Array<OneD, const NekDouble>& inarray,
574  Array<OneD, NekDouble>& outarray)
575  {
576  StdPyrExp::v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
577  }
578 
579  /**
580  * @param inarray Function evaluated at physical collocation
581  * points.
582  * @param outarray Inner product with respect to each basis
583  * function over the element.
584  */
585  void StdPyrExp::v_IProductWRTDerivBase_SumFac(
586  const int dir,
587  const Array<OneD, const NekDouble>& inarray,
588  Array<OneD, NekDouble>& outarray)
589  {
590  int i;
591  int nquad0 = m_base[0]->GetNumPoints();
592  int nquad1 = m_base[1]->GetNumPoints();
593  int nquad2 = m_base[2]->GetNumPoints();
594  int nqtot = nquad0*nquad1*nquad2;
595 
596  Array<OneD, NekDouble> gfac0(nquad0);
597  Array<OneD, NekDouble> gfac1(nquad1);
598  Array<OneD, NekDouble> gfac2(nquad2);
599  Array<OneD, NekDouble> tmp0 (nqtot);
600 
601 
602  int order0 = m_base[0]->GetNumModes();
603  int order1 = m_base[1]->GetNumModes();
604 
605  Array<OneD, NekDouble> wsp(nquad1*nquad2*order0 +
606  nquad2*order0*order1);
607 
608  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
609  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
610  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
611 
612  // set up geometric factor: (1+z0)/2
613  for(i = 0; i < nquad0; ++i)
614  {
615  gfac0[i] = 0.5*(1+z0[i]);
616  }
617 
618  // set up geometric factor: (1+z1)/2
619  for(i = 0; i < nquad1; ++i)
620  {
621  gfac1[i] = 0.5*(1+z1[i]);
622  }
623 
624  // Set up geometric factor: 2/(1-z2)
625  for(i = 0; i < nquad2; ++i)
626  {
627  gfac2[i] = 2.0/(1-z2[i]);
628  }
629 
630  // Derivative in first/second direction is always scaled as follows
631  const int nq01 = nquad0*nquad1;
632  for(i = 0; i < nquad2; ++i)
633  {
634  Vmath::Smul(nq01, gfac2[i],
635  &inarray[0] + i*nq01, 1,
636  &tmp0 [0] + i*nq01, 1);
637  }
638 
639  MultiplyByStdQuadratureMetric(tmp0, tmp0);
640 
641  switch(dir)
642  {
643  case 0:
644  {
645  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
646  m_base[1]->GetBdata (),
647  m_base[2]->GetBdata (),
648  tmp0, outarray, wsp,
649  false, true, true);
650  break;
651  }
652  case 1:
653  {
654  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata (),
655  m_base[1]->GetDbdata(),
656  m_base[2]->GetBdata (),
657  tmp0, outarray, wsp,
658  true, false, true);
659  break;
660  }
661  case 2:
662  {
663  Array<OneD, NekDouble> tmp3(m_ncoeffs);
664  Array<OneD, NekDouble> tmp4(m_ncoeffs);
665 
666  // Scale eta_1 derivative by gfac0
667  for(i = 0; i < nquad1*nquad2; ++i)
668  {
669  Vmath::Vmul(nquad0, tmp0 .get() + i*nquad0, 1,
670  gfac0.get(), 1,
671  tmp0 .get() + i*nquad0, 1);
672  }
673  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
674  m_base[1]->GetBdata(),
675  m_base[2]->GetBdata(),
676  tmp0, tmp3, wsp,
677  false, true, true);
678 
679  // Scale eta_2 derivative by gfac1*gfac2
680  for(i = 0; i < nquad2; ++i)
681  {
682  Vmath::Smul(nq01, gfac2[i],
683  &inarray[0] + i*nq01, 1,
684  &tmp0 [0] + i*nq01, 1);
685  }
686  for(i = 0; i < nquad1*nquad2; ++i)
687  {
688  Vmath::Smul(nquad0, gfac1[i%nquad1],
689  &tmp0[0] + i*nquad0, 1,
690  &tmp0[0] + i*nquad0, 1);
691  }
692 
693  MultiplyByStdQuadratureMetric(tmp0, tmp0);
694  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
695  m_base[1]->GetDbdata(),
696  m_base[2]->GetBdata(),
697  tmp0, tmp4, wsp,
698  true, false, true);
699 
700  MultiplyByStdQuadratureMetric(inarray,tmp0);
701  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
702  m_base[1]->GetBdata(),
703  m_base[2]->GetDbdata(),
704  tmp0,outarray,wsp,
705  true, true, false);
706 
707  Vmath::Vadd(m_ncoeffs,&tmp3[0],1,&outarray[0],1,&outarray[0],1);
708  Vmath::Vadd(m_ncoeffs,&tmp4[0],1,&outarray[0],1,&outarray[0],1);
709  break;
710  }
711  default:
712  {
713  ASSERTL1(false, "input dir is out of range");
714  break;
715  }
716  }
717  }
718 
719  //---------------------------------------
720  // Evaluation functions
721  //---------------------------------------
722 
723  void StdPyrExp::v_LocCoordToLocCollapsed(
724  const Array<OneD, const NekDouble>& xi,
725  Array<OneD, NekDouble>& eta)
726  {
727  if (fabs(xi[2]-1.0) < NekConstants::kNekZeroTol)
728  {
729  // Very top point of the pyramid
730  eta[0] = -1.0;
731  eta[1] = -1.0;
732  eta[2] = xi[2];
733  }
734  else
735  {
736  // Below the line-singularity -- Common case
737  eta[2] = xi[2]; // eta_z = xi_z
738  eta[1] = 2.0*(1.0 + xi[1])/(1.0 - xi[2]) - 1.0;
739  eta[0] = 2.0*(1.0 + xi[0])/(1.0 - xi[2]) - 1.0;
740  }
741  }
742 
743  void StdPyrExp::v_GetCoords(Array<OneD, NekDouble> &xi_x,
744  Array<OneD, NekDouble> &xi_y,
745  Array<OneD, NekDouble> &xi_z)
746  {
747  Array<OneD, const NekDouble> etaBar_x = m_base[0]->GetZ();
748  Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
749  Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
750  int Qx = GetNumPoints(0);
751  int Qy = GetNumPoints(1);
752  int Qz = GetNumPoints(2);
753 
754  // Convert collapsed coordinates into cartesian coordinates: eta --> xi
755  for (int k = 0; k < Qz; ++k )
756  {
757  for (int j = 0; j < Qy; ++j)
758  {
759  for (int i = 0; i < Qx; ++i)
760  {
761  int s = i + Qx*(j + Qy*k);
762 
763  xi_z[s] = eta_z[k];
764  xi_y[s] = (1.0 + eta_y[j]) * (1.0 - eta_z[k]) / 2.0 - 1.0;
765  xi_x[s] = (1.0 + etaBar_x[i]) * (1.0 - eta_z[k]) / 2.0 - 1.0;
766  }
767  }
768  }
769  }
770 
771  void StdPyrExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
772  {
773  Array<OneD, NekDouble> tmp(m_ncoeffs, 0.0);
774  tmp[mode] = 1.0;
775  v_BwdTrans(tmp, outarray);
776  }
777 
778  void StdPyrExp::v_GetFaceNumModes(
779  const int fid,
780  const Orientation faceOrient,
781  int &numModes0,
782  int &numModes1)
783  {
784  int nummodes [3] = {m_base[0]->GetNumModes(),
785  m_base[1]->GetNumModes(),
786  m_base[2]->GetNumModes()};
787  switch(fid)
788  {
789  // quad
790  case 0:
791  {
792  numModes0 = nummodes[0];
793  numModes1 = nummodes[1];
794  }
795  break;
796  case 1:
797  case 3:
798  {
799  numModes0 = nummodes[0];
800  numModes1 = nummodes[2];
801  }
802  break;
803  case 2:
804  case 4:
805  {
806  numModes0 = nummodes[1];
807  numModes1 = nummodes[2];
808  }
809  break;
810  }
811 
812  if ( faceOrient >= 9 )
813  {
814  std::swap(numModes0, numModes1);
815  }
816  }
817 
818  //---------------------------------------
819  // Helper functions
820  //---------------------------------------
821 
822  int StdPyrExp::v_GetNverts() const
823  {
824  return 5;
825  }
826 
827  int StdPyrExp::v_GetNedges() const
828  {
829  return 8;
830  }
831 
832  int StdPyrExp::v_GetNfaces() const
833  {
834  return 5;
835  }
836 
837  LibUtilities::ShapeType StdPyrExp::v_DetShapeType() const
838  {
839  return LibUtilities::ePyramid;
840  }
841 
842  int StdPyrExp::v_NumBndryCoeffs() const
843  {
844  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
845  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
846  "BasisType is not a boundary interior form");
847  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
848  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
849  "BasisType is not a boundary interior form");
850  ASSERTL1(GetBasisType(2) == LibUtilities::eModifiedPyr_C ||
851  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
852  "BasisType is not a boundary interior form");
853 
854  int P = m_base[0]->GetNumModes();
855  int Q = m_base[1]->GetNumModes();
856  int R = m_base[2]->GetNumModes();
857 
858  return LibUtilities::StdPyrData::
859  getNumberOfBndCoefficients(P, Q, R);
860  }
861 
862  int StdPyrExp::v_NumDGBndryCoeffs() const
863  {
864  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
865  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
866  "BasisType is not a boundary interior form");
867  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
868  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
869  "BasisType is not a boundary interior form");
870  ASSERTL1(GetBasisType(2) == LibUtilities::eModifiedPyr_C ||
871  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
872  "BasisType is not a boundary interior form");
873 
874  int P = m_base[0]->GetNumModes()-1;
875  int Q = m_base[1]->GetNumModes()-1;
876  int R = m_base[2]->GetNumModes()-1;
877 
878  return (P+1)*(Q+1) // 1 rect. face on base
879  + 2*(R+1) + P*(1 + 2*R - P) // 2 tri. (P,R) faces
880  + 2*(R+1) + Q*(1 + 2*R - Q); // 2 tri. (Q,R) faces
881  }
882 
883  int StdPyrExp::v_GetEdgeNcoeffs(const int i) const
884  {
885  ASSERTL2(i >= 0 && i <= 7, "edge id is out of range");
886 
887  if (i == 0 || i == 2)
888  {
889  return GetBasisNumModes(0);
890  }
891  else if (i == 1 || i == 3)
892  {
893  return GetBasisNumModes(1);
894  }
895  else
896  {
897  return GetBasisNumModes(2);
898  }
899  }
900 
901  int StdPyrExp::v_GetFaceNcoeffs(const int i) const
902  {
903  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
904 
905  if (i == 0)
906  {
907  return GetBasisNumModes(0)*GetBasisNumModes(1);
908  }
909  else if (i == 1 || i == 3)
910  {
911  int P = GetBasisNumModes(0)-1, Q = GetBasisNumModes(2)-1;
912  return Q+1 + (P*(1 + 2*Q - P))/2;
913  }
914  else
915  {
916  int P = GetBasisNumModes(1)-1, Q = GetBasisNumModes(2)-1;
917  return Q+1 + (P*(1 + 2*Q - P))/2;
918  }
919  }
920 
921  int StdPyrExp::v_GetFaceIntNcoeffs(const int i) const
922  {
923  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
924 
925  int P = m_base[0]->GetNumModes()-1;
926  int Q = m_base[1]->GetNumModes()-1;
927  int R = m_base[2]->GetNumModes()-1;
928 
929  if (i == 0)
930  {
931  return (P-1)*(Q-1);
932  }
933  else if (i == 1 || i == 3)
934  {
935  return (P-1) * (2*(R-1) - (P-1) - 1) / 2;
936  }
937  else
938  {
939  return (Q-1) * (2*(R-1) - (Q-1) - 1) / 2;
940  }
941  }
942 
943  int StdPyrExp::v_GetFaceNumPoints(const int i) const
944  {
945  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
946 
947  if (i == 0)
948  {
949  return m_base[0]->GetNumPoints()*
950  m_base[1]->GetNumPoints();
951  }
952  else if (i == 1 || i == 3)
953  {
954  return m_base[0]->GetNumPoints()*
955  m_base[2]->GetNumPoints();
956  }
957  else
958  {
959  return m_base[1]->GetNumPoints()*
960  m_base[2]->GetNumPoints();
961  }
962  }
963 
964 
965  const LibUtilities::BasisKey StdPyrExp::v_DetFaceBasisKey(
966  const int i, const int k) const
967  {
968  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
969  ASSERTL2(k >= 0 && k <= 1, "basis key id is out of range");
970 
971  switch(i)
972  {
973  case 0:
974  {
975  return EvaluateQuadFaceBasisKey(k,
976  m_base[k]->GetBasisType(),
977  m_base[k]->GetNumPoints(),
978  m_base[k]->GetNumModes());
979 
980  }
981  case 1:
982  case 3:
983  {
984  return EvaluateTriFaceBasisKey(k,
985  m_base[2*k]->GetBasisType(),
986  m_base[2*k]->GetNumPoints(),
987  m_base[2*k]->GetNumModes());
988  }
989  case 2:
990  case 4:
991  {
992  return EvaluateTriFaceBasisKey(k,
993  m_base[k+1]->GetBasisType(),
994  m_base[k+1]->GetNumPoints(),
995  m_base[k+1]->GetNumModes());
996  }
997  }
998 
999  // Should never get here.
1000  return LibUtilities::NullBasisKey;
1001  }
1002 
1003  int StdPyrExp::v_CalcNumberOfCoefficients(
1004  const std::vector<unsigned int> &nummodes,
1005  int &modes_offset)
1006  {
1007  int nmodes = LibUtilities::StdPyrData::getNumberOfCoefficients(
1008  nummodes[modes_offset],
1009  nummodes[modes_offset+1],
1010  nummodes[modes_offset+2]);
1011 
1012  modes_offset += 3;
1013  return nmodes;
1014  }
1015 
1016  LibUtilities::BasisType StdPyrExp::v_GetEdgeBasisType(const int i) const
1017  {
1018  ASSERTL2(i >= 0 && i <= 7, "edge id is out of range");
1019  if (i == 0 || i == 2)
1020  {
1021  return GetBasisType(0);
1022  }
1023  else if (i == 1 || i == 3)
1024  {
1025  return GetBasisType(1);
1026  }
1027  else
1028  {
1029  return GetBasisType(2);
1030  }
1031  }
1032 
1033 
1034  //---------------------------------------
1035  // Mappings
1036  //---------------------------------------
1037 
1038  void StdPyrExp::v_GetFaceToElementMap(
1039  const int fid,
1040  const Orientation faceOrient,
1041  Array<OneD, unsigned int> &maparray,
1042  Array<OneD, int> &signarray,
1043  int P,
1044  int Q)
1045  {
1046  ASSERTL1(GetEdgeBasisType(0) == GetEdgeBasisType(1),
1047  "Method only implemented if BasisType is identical"
1048  "in x and y directions");
1049  ASSERTL1(GetEdgeBasisType(0) == LibUtilities::eModified_A &&
1050  GetEdgeBasisType(4) == LibUtilities::eModifiedPyr_C,
1051  "Method only implemented for Modified_A BasisType"
1052  "(x and y direction) and ModifiedPyr_C BasisType (z "
1053  "direction)");
1054 
1055  int i, j, k, p, q, r, nFaceCoeffs, idx = 0;
1056  int nummodesA=0, nummodesB=0;
1057 
1058  int order0 = m_base[0]->GetNumModes();
1059  int order1 = m_base[1]->GetNumModes();
1060  int order2 = m_base[2]->GetNumModes();
1061 
1062  switch (fid)
1063  {
1064  case 0:
1065  nummodesA = order0;
1066  nummodesB = order1;
1067  break;
1068  case 1:
1069  case 3:
1070  nummodesA = order0;
1071  nummodesB = order2;
1072  break;
1073  case 2:
1074  case 4:
1075  nummodesA = order1;
1076  nummodesB = order2;
1077  break;
1078  default:
1079  ASSERTL0(false,"fid must be between 0 and 4");
1080  }
1081 
1082  bool CheckForZeroedModes = false;
1083 
1084  if (P == -1)
1085  {
1086  P = nummodesA;
1087  Q = nummodesB;
1088  nFaceCoeffs = GetFaceNcoeffs(fid);
1089  }
1090  else if (fid > 0)
1091  {
1092  nFaceCoeffs = P*(2*Q-P+1)/2;
1093  CheckForZeroedModes = true;
1094  }
1095  else
1096  {
1097  nFaceCoeffs = P*Q;
1098  CheckForZeroedModes = true;
1099  }
1100 
1101  // Allocate the map array and sign array; set sign array to ones (+)
1102  if (maparray.num_elements() != nFaceCoeffs)
1103  {
1104  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1105  }
1106 
1107  if (signarray.num_elements() != nFaceCoeffs)
1108  {
1109  signarray = Array<OneD, int>(nFaceCoeffs,1);
1110  }
1111  else
1112  {
1113  fill(signarray.get(), signarray.get() + nFaceCoeffs, 1);
1114  }
1115 
1116  // Set up an array indexing for quads, since the ordering may need
1117  // to be transposed.
1118  Array<OneD, int> arrayindx(nFaceCoeffs,-1);
1119 
1120  if (fid == 0)
1121  {
1122  for (i = 0; i < Q; i++)
1123  {
1124  for (j = 0; j < P; j++)
1125  {
1126  if (faceOrient < 9)
1127  {
1128  arrayindx[i*P+j] = i*P+j;
1129  }
1130  else
1131  {
1132  arrayindx[i*P+j] = j*Q+i;
1133  }
1134  }
1135  }
1136  }
1137 
1138  // Set up ordering inside each 2D face. Also for triangular faces,
1139  // populate signarray.
1140  switch (fid)
1141  {
1142  case 0: // Bottom quad
1143 
1144  for (q = 0; q < Q; ++q)
1145  {
1146  for (p = 0; p < P; ++p)
1147  {
1148  maparray[arrayindx[q*P+p]] = GetMode(p,q,0);
1149  }
1150  }
1151  break;
1152 
1153  case 1: // Front triangle
1154  for (p = 0; p < P; ++p)
1155  {
1156  for (r = 0; r < Q-p; ++r)
1157  {
1158  if ((int)faceOrient == 7 && p > 1)
1159  {
1160  signarray[idx] = p % 2 ? -1 : 1;
1161  }
1162  maparray[idx++] = GetMode(p,0,r);
1163  }
1164  }
1165  break;
1166 
1167  case 2: // Right triangle
1168  maparray[idx++] = GetMode(1,0,0);
1169  maparray[idx++] = GetMode(0,0,1);
1170  for (r = 1; r < Q-1; ++r)
1171  {
1172  maparray[idx++] = GetMode(1,0,r);
1173  }
1174 
1175  for (q = 1; q < P; ++q)
1176  {
1177  for (r = 0; r < Q-q; ++r)
1178  {
1179  if ((int)faceOrient == 7 && q > 1)
1180  {
1181  signarray[idx] = q % 2 ? -1 : 1;
1182  }
1183  maparray[idx++] = GetMode(1,q,r);
1184  }
1185  }
1186  break;
1187 
1188  case 3: // Rear triangle
1189  maparray[idx++] = GetMode(0,1,0);
1190  maparray[idx++] = GetMode(0,0,1);
1191  for (r = 1; r < Q-1; ++r)
1192  {
1193  maparray[idx++] = GetMode(0,1,r);
1194  }
1195 
1196  for (p = 1; p < P; ++p)
1197  {
1198  for (r = 0; r < Q-p; ++r)
1199  {
1200  if ((int)faceOrient == 7 && p > 1)
1201  {
1202  signarray[idx] = p % 2 ? -1 : 1;
1203  }
1204  maparray[idx++] = GetMode(p, 1, r);
1205  }
1206  }
1207  break;
1208 
1209  case 4: // Left triangle
1210  for (q = 0; q < P; ++q)
1211  {
1212  for (r = 0; r < Q-q; ++r)
1213  {
1214  if ((int)faceOrient == 7 && q > 1)
1215  {
1216  signarray[idx] = q % 2 ? -1 : 1;
1217  }
1218  maparray[idx++] = GetMode(0,q,r);
1219  }
1220  }
1221  break;
1222 
1223  default:
1224  ASSERTL0(false, "Face to element map unavailable.");
1225  }
1226 
1227  if (fid > 0)
1228  {
1229  if(CheckForZeroedModes)
1230  {
1231  // zero signmap and set maparray to zero if elemental
1232  // modes are not as large as face modesl
1233  int idx = 0;
1234  for (j = 0; j < P; ++j)
1235  {
1236  idx += Q-j;
1237  for (k = Q-j; k < Q-j; ++k)
1238  {
1239  signarray[idx] = 0.0;
1240  maparray[idx++] = maparray[0];
1241  }
1242  }
1243 
1244  for (j = P; j < P; ++j)
1245  {
1246  for (k = 0; k < Q-j; ++k)
1247  {
1248  signarray[idx] = 0.0;
1249  maparray[idx++] = maparray[0];
1250  }
1251  }
1252  }
1253 
1254  // Triangles only have one possible orientation (base
1255  // direction reversed); swap edge modes.
1256  if ((int)faceOrient == 7)
1257  {
1258  swap(maparray[0], maparray[Q]);
1259  for (i = 1; i < Q-1; ++i)
1260  {
1261  swap(maparray[i+1], maparray[Q+i]);
1262  }
1263  }
1264  }
1265  else
1266  {
1267  if(CheckForZeroedModes)
1268  {
1269  // zero signmap and set maparray to zero if elemental
1270  // modes are not as large as face modesl
1271  for (j = 0; j < P; ++j)
1272  {
1273  for (k = Q; k < Q; ++k)
1274  {
1275  signarray[arrayindx[j+k*P]] = 0.0;
1276  maparray[arrayindx[j+k*P]] = maparray[0];
1277  }
1278  }
1279 
1280  for (j = P; j < P; ++j)
1281  {
1282  for (k = 0; k < Q; ++k)
1283  {
1284  signarray[arrayindx[j+k*P]] = 0.0;
1285  maparray[arrayindx[j+k*P]] = maparray[0];
1286  }
1287  }
1288  }
1289 
1290  // The code below is exactly the same as that taken from
1291  // StdHexExp and reverses the 'b' and 'a' directions as
1292  // appropriate (1st and 2nd if statements respectively) in
1293  // quadrilateral faces.
1294  if (faceOrient == 6 || faceOrient == 8 ||
1295  faceOrient == 11 || faceOrient == 12)
1296  {
1297  if (faceOrient < 9)
1298  {
1299  for (i = 3; i < Q; i += 2)
1300  {
1301  for (j = 0; j < P; j++)
1302  {
1303  signarray[arrayindx[i*P+j]] *= -1;
1304  }
1305  }
1306 
1307  for (i = 0; i < P; i++)
1308  {
1309  swap(maparray [i], maparray [i+P]);
1310  swap(signarray[i], signarray[i+P]);
1311  }
1312  }
1313  else
1314  {
1315  for (i = 0; i < Q; i++)
1316  {
1317  for (j = 3; j < P; j += 2)
1318  {
1319  signarray[arrayindx[i*P+j]] *= -1;
1320  }
1321  }
1322 
1323  for (i = 0; i < Q; i++)
1324  {
1325  swap (maparray [i], maparray [i+Q]);
1326  swap (signarray[i], signarray[i+Q]);
1327  }
1328  }
1329  }
1330 
1331  if (faceOrient == 7 || faceOrient == 8 ||
1332  faceOrient == 10 || faceOrient == 12)
1333  {
1334  if (faceOrient < 9)
1335  {
1336  for (i = 0; i < Q; i++)
1337  {
1338  for (j = 3; j < P; j += 2)
1339  {
1340  signarray[arrayindx[i*P+j]] *= -1;
1341  }
1342  }
1343 
1344  for(i = 0; i < Q; i++)
1345  {
1346  swap(maparray [i*P], maparray [i*P+1]);
1347  swap(signarray[i*P], signarray[i*P+1]);
1348  }
1349  }
1350  else
1351  {
1352  for (i = 3; i < Q; i += 2)
1353  {
1354  for (j = 0; j < P; j++)
1355  {
1356  signarray[arrayindx[i*P+j]] *= -1;
1357  }
1358  }
1359 
1360  for (i = 0; i < P; i++)
1361  {
1362  swap(maparray [i*Q], maparray [i*Q+1]);
1363  swap(signarray[i*Q], signarray[i*Q+1]);
1364  }
1365  }
1366  }
1367  }
1368  }
1369 
1370  int StdPyrExp::v_GetVertexMap(int vId, bool useCoeffPacking)
1371  {
1372  ASSERTL0(GetEdgeBasisType(vId) == LibUtilities::eModified_A ||
1373  GetEdgeBasisType(vId) == LibUtilities::eModified_A ||
1374  GetEdgeBasisType(vId) == LibUtilities::eModifiedPyr_C,
1375  "Mapping not defined for this type of basis");
1376 
1377 
1378 
1379  int l = 0;
1380 
1381  if(useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1382  {
1383  switch (vId)
1384  {
1385  case 0:
1386  l = GetMode(0,0,0);
1387  break;
1388  case 1:
1389  l = GetMode(0,0,1);
1390  break;
1391  case 2:
1392  l = GetMode(0,1,0);
1393  break;
1394  case 3:
1395  l = GetMode(1,0,0);
1396  break;
1397  case 4:
1398  l = GetMode(1,1,0);
1399  break;
1400  default:
1401  ASSERTL0(false, "local vertex id must be between 0 and 4");
1402  }
1403  }
1404  else
1405  {
1406  switch (vId)
1407  {
1408  case 0:
1409  l = GetMode(0,0,0);
1410  break;
1411  case 1:
1412  l = GetMode(1,0,0);
1413  break;
1414  case 2:
1415  l = GetMode(1,1,0);
1416  break;
1417  case 3:
1418  l = GetMode(0,1,0);
1419  break;
1420  case 4:
1421  l = GetMode(0,0,1);
1422  break;
1423  default:
1424  ASSERTL0(false, "local vertex id must be between 0 and 4");
1425  }
1426  }
1427 
1428  return l;
1429  }
1430 
1431  void StdPyrExp::v_GetEdgeInteriorMap(
1432  const int eid,
1433  const Orientation edgeOrient,
1434  Array<OneD, unsigned int> &maparray,
1435  Array<OneD, int> &signarray)
1436  {
1437  int i;
1438  bool signChange;
1439  const int P = m_base[0]->GetNumModes() - 1;
1440  const int Q = m_base[1]->GetNumModes() - 1;
1441  const int R = m_base[2]->GetNumModes() - 1;
1442  const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid) - 2;
1443 
1444  if (maparray.num_elements() != nEdgeIntCoeffs)
1445  {
1446  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1447  }
1448 
1449  if(signarray.num_elements() != nEdgeIntCoeffs)
1450  {
1451  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1452  }
1453  else
1454  {
1455  fill(signarray.get(), signarray.get()+nEdgeIntCoeffs, 1);
1456  }
1457 
1458  // If edge is oriented backwards, change sign of modes which have
1459  // degree 2n+1, n >= 1.
1460  signChange = edgeOrient == eBackwards;
1461 
1462  switch (eid)
1463  {
1464  case 0:
1465  for (i = 2; i <= P; ++i)
1466  {
1467  maparray[i-2] = GetMode(i,0,0);
1468  }
1469  break;
1470 
1471  case 1:
1472  for (i = 2; i <= Q; ++i)
1473  {
1474  maparray[i-2] = GetMode(1,i,0);
1475  }
1476  break;
1477  case 2:
1478  for (i = 2; i <= P; ++i)
1479  {
1480  maparray[i-2] = GetMode(i,1,0);
1481  }
1482  break;
1483 
1484  case 3:
1485  for (i = 2; i <= Q; ++i)
1486  {
1487  maparray[i-2] = GetMode(0,i,0);
1488  }
1489  break;
1490  case 4:
1491  for (i = 2; i <= R; ++i)
1492  {
1493  maparray[i-2] = GetMode(0,0,i);
1494  }
1495  break;
1496 
1497  case 5:
1498  for (i = 1; i <= R-1; ++i)
1499  {
1500  maparray[i-1] = GetMode(1,0,i);
1501  }
1502  break;
1503  case 6:
1504  for (i = 1; i <= R-1; ++i)
1505  {
1506  maparray[i-1] = GetMode(1,1,i);
1507  }
1508  break;
1509 
1510  case 7:
1511  for (i = 1; i <= R-1; ++i)
1512  {
1513  maparray[i-1] = GetMode(0,1,i);
1514  }
1515  break;
1516  default:
1517  ASSERTL0(false, "Edge not defined.");
1518  break;
1519  }
1520 
1521  if (signChange)
1522  {
1523  for (i = 1; i < nEdgeIntCoeffs; i += 2)
1524  {
1525  signarray[i] = -1;
1526  }
1527  }
1528  }
1529 
1530  void StdPyrExp::v_GetFaceInteriorMap(
1531  const int fid,
1532  const Orientation faceOrient,
1533  Array<OneD, unsigned int> &maparray,
1534  Array<OneD, int> &signarray)
1535  {
1536  const int P = m_base[0]->GetNumModes() - 1;
1537  const int Q = m_base[1]->GetNumModes() - 1;
1538  const int R = m_base[2]->GetNumModes() - 1;
1539  const int nFaceIntCoeffs = v_GetFaceIntNcoeffs(fid);
1540  int p, q, r, idx = 0;
1541  int nummodesA = 0;
1542  int nummodesB = 0;
1543  int i, j;
1544 
1545  if (maparray.num_elements() != nFaceIntCoeffs)
1546  {
1547  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1548  }
1549 
1550  if (signarray.num_elements() != nFaceIntCoeffs)
1551  {
1552  signarray = Array<OneD, int>(nFaceIntCoeffs, 1);
1553  }
1554  else
1555  {
1556  fill(signarray.get(), signarray.get()+nFaceIntCoeffs, 1);
1557  }
1558 
1559  // Set up an array indexing for quad faces, since the ordering may
1560  // need to be transposed depending on orientation.
1561  Array<OneD, int> arrayindx(nFaceIntCoeffs);
1562  if (fid == 0)
1563  {
1564  nummodesA = P-1;
1565  nummodesB = Q-1;
1566 
1567  for (i = 0; i < nummodesB; i++)
1568  {
1569  for (j = 0; j < nummodesA; j++)
1570  {
1571  if (faceOrient < 9)
1572  {
1573  arrayindx[i*nummodesA+j] = i*nummodesA+j;
1574  }
1575  else
1576  {
1577  arrayindx[i*nummodesA+j] = j*nummodesB+i;
1578  }
1579  }
1580  }
1581  }
1582 
1583  switch (fid)
1584  {
1585  case 0: // Bottom quad
1586  for (q = 2; q <= Q; ++q)
1587  {
1588  for (p = 2; p <= P; ++p)
1589  {
1590  maparray[arrayindx[(q-2)*nummodesA+(p-2)]] = GetMode(p,q,0);
1591  }
1592  }
1593  break;
1594  case 1: // Front triangle
1595  for (p = 2; p <= P; ++p)
1596  {
1597  for (r = 1; r <= R-p; ++r)
1598  {
1599  if ((int)faceOrient == 7)
1600  {
1601  signarray[idx] = p % 2 ? -1 : 1;
1602  }
1603  maparray[idx++] = GetMode(p,0,r);
1604  }
1605  }
1606  break;
1607  case 2: // Right triangle
1608  for (q = 2; q <= Q; ++q)
1609  {
1610  for (r = 1; r <= R-q; ++r)
1611  {
1612  if ((int)faceOrient == 7)
1613  {
1614  signarray[idx] = q % 2 ? -1 : 1;
1615  }
1616  maparray[idx++] = GetMode(1, q, r);
1617  }
1618  }
1619  break;
1620 
1621  case 3: // Rear triangle
1622  for (p = 2; p <= P; ++p)
1623  {
1624  for (r = 1; r <= R-p; ++r)
1625  {
1626  if ((int)faceOrient == 7)
1627  {
1628  signarray[idx] = p % 2 ? -1 : 1;
1629  }
1630  maparray[idx++] = GetMode(p, 1, r);
1631  }
1632  }
1633  break;
1634 
1635  case 4: // Left triangle
1636  for (q = 2; q <= Q; ++q)
1637  {
1638  for (r = 1; r <= R-q; ++r)
1639  {
1640  if ((int)faceOrient == 7)
1641  {
1642  signarray[idx] = q % 2 ? -1 : 1;
1643  }
1644  maparray[idx++] = GetMode(0, q, r);
1645  }
1646  }
1647  break;
1648  default:
1649  ASSERTL0(false, "Face interior map not available.");
1650  }
1651 
1652  // Triangular faces are processed in the above switch loop; for
1653  // remaining quad faces, set up orientation if necessary.
1654  if (fid > 0)
1655  {
1656  return;
1657  }
1658 
1659  if (faceOrient == 6 || faceOrient == 8 ||
1660  faceOrient == 11 || faceOrient == 12)
1661  {
1662  if (faceOrient < 9)
1663  {
1664  for (i = 1; i < nummodesB; i += 2)
1665  {
1666  for (j = 0; j < nummodesA; j++)
1667  {
1668  signarray[arrayindx[i*nummodesA+j]] *= -1;
1669  }
1670  }
1671  }
1672  else
1673  {
1674  for (i = 0; i < nummodesB; i++)
1675  {
1676  for (j = 1; j < nummodesA; j += 2)
1677  {
1678  signarray[arrayindx[i*nummodesA+j]] *= -1;
1679  }
1680  }
1681  }
1682  }
1683 
1684  if (faceOrient == 7 || faceOrient == 8 ||
1685  faceOrient == 10 || faceOrient == 12)
1686  {
1687  if (faceOrient < 9)
1688  {
1689  for (i = 0; i < nummodesB; i++)
1690  {
1691  for (j = 1; j < nummodesA; j += 2)
1692  {
1693  signarray[arrayindx[i*nummodesA+j]] *= -1;
1694  }
1695  }
1696  }
1697  else
1698  {
1699  for (i = 1; i < nummodesB; i += 2)
1700  {
1701  for (j = 0; j < nummodesA; j++)
1702  {
1703  signarray[arrayindx[i*nummodesA+j]] *= -1;
1704  }
1705  }
1706  }
1707  }
1708  }
1709 
1710  void StdPyrExp::v_GetInteriorMap(Array<OneD, unsigned int> &outarray)
1711  {
1712  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1713  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1714  "BasisType is not a boundary interior form");
1715  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1716  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1717  "BasisType is not a boundary interior form");
1718  ASSERTL1(GetBasisType(2) == LibUtilities::eModifiedPyr_C ||
1719  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1720  "BasisType is not a boundary interior form");
1721 
1722 
1723  int P = m_base[0]->GetNumModes() - 1, p;
1724  int Q = m_base[1]->GetNumModes() - 1, q;
1725  int R = m_base[2]->GetNumModes() - 1, r;
1726 
1727  int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1728 
1729  if(outarray.num_elements()!=nIntCoeffs)
1730  {
1731  outarray = Array<OneD, unsigned int>(nIntCoeffs);
1732  }
1733 
1734  int idx = 0;
1735 
1736  // Loop over all interior modes.
1737  for (p = 2; p <= P; ++p)
1738  {
1739  for (q = 2; q <= Q; ++q)
1740  {
1741  int maxpq = max(p,q);
1742  for (r = 1; r <= R-maxpq; ++r)
1743  {
1744  outarray[idx++] = GetMode(p,q,r);
1745  }
1746  }
1747  }
1748  }
1749 
1750  void StdPyrExp::v_GetBoundaryMap(Array<OneD, unsigned int> &maparray)
1751  {
1752  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1753  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1754  "BasisType is not a boundary interior form");
1755  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1756  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1757  "BasisType is not a boundary interior form");
1758  ASSERTL1(GetBasisType(2) == LibUtilities::eModifiedPyr_C ||
1759  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1760  "BasisType is not a boundary interior form");
1761 
1762  int P = m_base[0]->GetNumModes() - 1, p;
1763  int Q = m_base[1]->GetNumModes() - 1, q;
1764  int R = m_base[2]->GetNumModes() - 1, r;
1765  int idx = 0;
1766 
1767  int nBnd = NumBndryCoeffs();
1768 
1769  if (maparray.num_elements() != nBnd)
1770  {
1771  maparray = Array<OneD, unsigned int>(nBnd);
1772  }
1773 
1774  // Loop over all boundary modes (in ascending order).
1775  for (p = 0; p <= P; ++p)
1776  {
1777  // First two q-r planes are entirely boundary modes.
1778  if (p <= 1)
1779  {
1780  for (q = 0; q <= Q; ++q)
1781  {
1782  int maxpq = max(p,q);
1783  for (r = 0; r <= R-maxpq; ++r)
1784  {
1785  maparray[idx++] = GetMode(p,q,r);
1786  }
1787  }
1788  }
1789  else
1790  {
1791  // Remaining q-r planes contain boundary modes on the two
1792  // front and back sides and edges 0 2.
1793  for (q = 0; q <= Q; ++q)
1794  {
1795  if (q <= 1)
1796  {
1797  for (r = 0; r <= R-p; ++r)
1798  {
1799  maparray[idx++] = GetMode(p,q,r);
1800  }
1801  }
1802  else
1803  {
1804  maparray[idx++] = GetMode(p,q,0);
1805  }
1806  }
1807  }
1808  }
1809  }
1810 
1811  //---------------------------------------
1812  // Wrapper functions
1813  //---------------------------------------
1814 
1815  DNekMatSharedPtr StdPyrExp::v_GenMatrix(const StdMatrixKey &mkey)
1816  {
1817  return CreateGeneralMatrix(mkey);
1818  }
1819 
1820  DNekMatSharedPtr StdPyrExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
1821  {
1822  return v_GenMatrix(mkey);
1823  }
1824 
1825 
1826  /**
1827  * @brief Compute the mode number in the expansion for a
1828  * particular tensorial combination.
1829  *
1830  * Modes are numbered with the r index travelling fastest,
1831  * followed by q and then p, and each q-r plane is of size
1832  *
1833  * (R+1-p)*(Q+1) - l(l+1)/2 where l = max(0,Q-p)
1834  *
1835  * For example, when P=2, Q=3 and R=4 the indexing inside each
1836  * q-r plane (with r increasing upwards and q to the right)
1837  * is:
1838  *
1839  * p = 0: p = 1: p = 2:
1840  * ----------------------------------
1841  * 4
1842  * 3 8 17 21
1843  * 2 7 11 16 20 24 29 32 35
1844  * 1 6 10 13 15 19 23 26 28 31 34 37
1845  * 0 5 9 12 14 18 22 25 27 30 33 36
1846  *
1847  * Note that in this element, we must have that \f$ P,Q \leq
1848  * R\f$.
1849  */
1850  int StdPyrExp::GetMode(const int I, const int J, const int K)
1851  {
1852  const int Q = m_base[1]->GetNumModes()-1;
1853  const int R = m_base[2]->GetNumModes()-1;
1854 
1855  int i,l;
1856  int cnt = 0;
1857 
1858  // Traverse to q-r plane number I
1859  for (i = 0; i < I; ++i)
1860  {
1861  // Size of triangle part
1862  l = max(0,Q-i);
1863 
1864  // Size of rectangle part
1865  cnt += (R+1-i)*(Q+1) - l*(l+1)/2;
1866  }
1867 
1868  // Traverse to q column J (Pretend this is a face of width J)
1869  l = max(0,J-1-I);
1870  cnt += (R+1-I)*J - l*(l+1)/2;
1871 
1872  // Traverse up stacks to K
1873  cnt += K;
1874 
1875  return cnt;
1876  }
1877 
1878 
1879  void StdPyrExp::v_MultiplyByStdQuadratureMetric(
1880  const Array<OneD, const NekDouble>& inarray,
1881  Array<OneD, NekDouble>& outarray)
1882  {
1883  int i, j;
1884 
1885  int nquad0 = m_base[0]->GetNumPoints();
1886  int nquad1 = m_base[1]->GetNumPoints();
1887  int nquad2 = m_base[2]->GetNumPoints();
1888 
1889  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1890  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1891  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
1892 
1893  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1894 
1895  // Multiply by integration constants in x-direction
1896  for(i = 0; i < nquad1*nquad2; ++i)
1897  {
1898  Vmath::Vmul(nquad0, inarray.get()+i*nquad0, 1,
1899  w0.get(), 1, outarray.get()+i*nquad0,1);
1900  }
1901 
1902  // Multiply by integration constants in y-direction
1903  for(j = 0; j < nquad2; ++j)
1904  {
1905  for(i = 0; i < nquad1; ++i)
1906  {
1907  Blas::Dscal(nquad0,w1[i], &outarray[0]+i*nquad0 +
1908  j*nquad0*nquad1,1);
1909  }
1910  }
1911 
1912  // Multiply by integration constants in z-direction; need to
1913  // incorporate factor [(1-eta_3)/2]^2 into weights, but only if
1914  // using GLL quadrature points.
1915  switch(m_base[2]->GetPointsType())
1916  {
1917  // (2,0) Jacobi inner product.
1918  case LibUtilities::eGaussRadauMAlpha2Beta0:
1919  for(i = 0; i < nquad2; ++i)
1920  {
1921  Blas::Dscal(nquad0*nquad1, 0.25*w2[i],
1922  &outarray[0]+i*nquad0*nquad1, 1);
1923  }
1924  break;
1925 
1926  default:
1927  for(i = 0; i < nquad2; ++i)
1928  {
1929  Blas::Dscal(nquad0*nquad1,0.25*(1-z2[i])*(1-z2[i])*w2[i],
1930  &outarray[0]+i*nquad0*nquad1,1);
1931  }
1932  break;
1933  }
1934  }
1935 
1936 
1937  void StdPyrExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
1938  const StdMatrixKey &mkey)
1939  {
1940  // Generate an orthonogal expansion
1941  int qa = m_base[0]->GetNumPoints();
1942  int qb = m_base[1]->GetNumPoints();
1943  int qc = m_base[2]->GetNumPoints();
1944  int nmodes_a = m_base[0]->GetNumModes();
1945  int nmodes_b = m_base[1]->GetNumModes();
1946  int nmodes_c = m_base[2]->GetNumModes();
1947  // Declare orthogonal basis.
1948  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
1949  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
1950  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
1951 
1952  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
1953  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A,nmodes_b,pb);
1954  LibUtilities::BasisKey Bc(LibUtilities::eOrthoPyr_C,nmodes_c,pc);
1955  StdPyrExp OrthoExp(Ba,Bb,Bc);
1956 
1957  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1958  int i,j,k,cnt = 0;
1959 
1960  // project onto modal space.
1961  OrthoExp.FwdTrans(array,orthocoeffs);
1962 
1963  if(mkey.ConstFactorExists(eFactorSVVPowerKerDiffCoeff))
1964  {
1965  // Rodrigo's power kernel
1966  NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
1967  NekDouble SvvDiffCoeff =
1968  mkey.GetConstFactor(eFactorSVVPowerKerDiffCoeff)*
1969  mkey.GetConstFactor(eFactorSVVDiffCoeff);
1970 
1971  for(int i = 0; i < nmodes_a; ++i)
1972  {
1973  for(int j = 0; j < nmodes_b; ++j)
1974  {
1975  int maxij = max(i,j);
1976  NekDouble fac1 = std::max(
1977  pow((1.0*i)/(nmodes_a-1),cutoff*nmodes_a),
1978  pow((1.0*j)/(nmodes_b-1),cutoff*nmodes_b));
1979 
1980  for(int k = 0; k < nmodes_c-maxij; ++k)
1981  {
1982  NekDouble fac = std::max(fac1,
1983  pow((1.0*k)/(nmodes_c-1),cutoff*nmodes_c));
1984 
1985  orthocoeffs[cnt] *= SvvDiffCoeff * fac;
1986  cnt++;
1987  }
1988  }
1989  }
1990  }
1991  else if(mkey.ConstFactorExists(eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
1992  {
1993  NekDouble SvvDiffCoeff =
1994  mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff)*
1995  mkey.GetConstFactor(eFactorSVVDiffCoeff);
1996 
1997  int max_abc = max(nmodes_a-kSVVDGFiltermodesmin,
1998  nmodes_b-kSVVDGFiltermodesmin);
1999  max_abc = max(max_abc, nmodes_c-kSVVDGFiltermodesmin);
2000  // clamp max_abc
2001  max_abc = max(max_abc,0);
2002  max_abc = min(max_abc,kSVVDGFiltermodesmax-kSVVDGFiltermodesmin);
2003 
2004  for(int i = 0; i < nmodes_a; ++i)
2005  {
2006  for(int j = 0; j < nmodes_b; ++j)
2007  {
2008  int maxij = max(i,j);
2009 
2010  for(int k = 0; k < nmodes_c-maxij; ++k)
2011  {
2012  int maxijk = max(maxij,k);
2013  maxijk = min(maxijk,kSVVDGFiltermodesmax-1);
2014 
2015  orthocoeffs[cnt] *= SvvDiffCoeff *
2016  kSVVDGFilter[max_abc][maxijk];
2017  cnt++;
2018  }
2019  }
2020  }
2021  }
2022  else
2023  {
2024  //SVV filter paramaters (how much added diffusion relative
2025  // to physical one and fraction of modes from which you
2026  // start applying this added diffusion)
2027  //
2028  NekDouble SvvDiffCoeff = mkey.GetConstFactor(StdRegions::eFactorSVVDiffCoeff);
2029  NekDouble SVVCutOff = mkey.GetConstFactor(StdRegions::eFactorSVVCutoffRatio);
2030 
2031  //Defining the cut of mode
2032  int cutoff_a = (int) (SVVCutOff*nmodes_a);
2033  int cutoff_b = (int) (SVVCutOff*nmodes_b);
2034  int cutoff_c = (int) (SVVCutOff*nmodes_c);
2035  //To avoid the fac[j] from blowing up
2036  NekDouble epsilon = 1;
2037 
2038  int nmodes = min(min(nmodes_a,nmodes_b),nmodes_c);
2039  NekDouble cutoff = min(min(cutoff_a,cutoff_b),cutoff_c);
2040 
2041  for(i = 0; i < nmodes_a; ++i)//P
2042  {
2043  for(j = 0; j < nmodes_b; ++j) //Q
2044  {
2045  int maxij = max(i,j);
2046  for(k = 0; k < nmodes_c-maxij; ++k) //R
2047  {
2048  if(j + k >= cutoff || i + k >= cutoff)
2049  {
2050  orthocoeffs[cnt] *=
2051  (SvvDiffCoeff*exp(-(i+k-nmodes)*(i+k-nmodes)/
2052  ((NekDouble)((i+k-cutoff+epsilon)*
2053  (i+k-cutoff+epsilon))))*
2054  exp(-(j-nmodes)*(j-nmodes)/
2055  ((NekDouble)((j-cutoff+epsilon)*
2056  (j-cutoff+epsilon)))));
2057  }
2058  else
2059  {
2060  orthocoeffs[cnt] *= 0.0;
2061  }
2062  cnt++;
2063  }
2064  }
2065  }
2066  }
2067 
2068  // backward transform to physical space
2069  OrthoExp.BwdTrans(orthocoeffs,array);
2070  }
2071 
2072  void StdPyrExp::v_ReduceOrderCoeffs(
2073  int numMin,
2074  const Array<OneD, const NekDouble> &inarray,
2075  Array<OneD, NekDouble> &outarray)
2076  {
2077  int nquad0 = m_base[0]->GetNumPoints();
2078  int nquad1 = m_base[1]->GetNumPoints();
2079  int nquad2 = m_base[2]->GetNumPoints();
2080  int nqtot = nquad0*nquad1*nquad2;
2081  int nmodes0 = m_base[0]->GetNumModes();
2082  int nmodes1 = m_base[1]->GetNumModes();
2083  int nmodes2 = m_base[2]->GetNumModes();
2084  int numMax = nmodes0;
2085 
2086  Array<OneD, NekDouble> coeff (m_ncoeffs);
2087  Array<OneD, NekDouble> coeff_tmp1(m_ncoeffs, 0.0);
2088  Array<OneD, NekDouble> phys_tmp (nqtot, 0.0);
2089  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
2090 
2091 
2092  const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
2093  const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
2094  const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
2095 
2096  LibUtilities::BasisKey bortho0(
2097  LibUtilities::eOrtho_A, nmodes0, Pkey0);
2098  LibUtilities::BasisKey bortho1(
2099  LibUtilities::eOrtho_A, nmodes1, Pkey1);
2100  LibUtilities::BasisKey bortho2(
2101  LibUtilities::eOrthoPyr_C, nmodes2, Pkey2);
2102 
2103  int cnt = 0;
2104  int u = 0;
2105  int i = 0;
2106  StdRegions::StdPyrExpSharedPtr OrthoPyrExp;
2107 
2108  OrthoPyrExp = MemoryManager<StdRegions::StdPyrExp>
2109  ::AllocateSharedPtr(bortho0, bortho1, bortho2);
2110 
2111  BwdTrans(inarray,phys_tmp);
2112  OrthoPyrExp->FwdTrans(phys_tmp, coeff);
2113 
2114  // filtering
2115  for (u = 0; u < numMin; ++u)
2116  {
2117  for (i = 0; i < numMin; ++i)
2118  {
2119 
2120  int maxui = max(u,i);
2121  Vmath::Vcopy(numMin - maxui, tmp = coeff + cnt, 1,
2122  tmp2 = coeff_tmp1 + cnt, 1);
2123  cnt += nmodes2 - maxui;
2124  }
2125 
2126  for (i = numMin; i < nmodes1; ++i)
2127  {
2128  int maxui = max(u,i);
2129  cnt += numMax - maxui;
2130  }
2131  }
2132 
2133  OrthoPyrExp->BwdTrans(coeff_tmp1, phys_tmp);
2134  StdPyrExp::FwdTrans(phys_tmp, outarray);
2135  }
2136 
2137  }
2138 }
Nektar::StdRegions::StdPyrExp::v_GetFaceNcoeffs
virtual int v_GetFaceNcoeffs(const int i) const
Definition: StdPyrExp.cpp:901
Nektar::StdRegions::StdPyrExp::v_GetFaceInteriorMap
virtual void v_GetFaceInteriorMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
Definition: StdPyrExp.cpp:1530
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::StdRegions::kSVVDGFiltermodesmax
const int kSVVDGFiltermodesmax
Definition: StdRegions.hpp:385
Nektar::StdRegions::StdPyrExp::v_GetCoords
virtual void v_GetCoords(Array< OneD, NekDouble > &xi_x, Array< OneD, NekDouble > &xi_y, Array< OneD, NekDouble > &xi_z)
Definition: StdPyrExp.cpp:743
Nektar::StdRegions::StdPyrExp::v_MultiplyByStdQuadratureMetric
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdPyrExp.cpp:1879
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:787
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:469
Nektar::LibUtilities::eModifiedPyr_C
Principle Modified Functions .
Definition: BasisType.h:52
Nektar::StdRegions::StdExpansion::GetNumPoints
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:228
Nektar::StdRegions::StdExpansion::GetBasisNumModes
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
Definition: StdExpansion.h:177
Nektar::StdRegions::StdExpansion::MultiplyByStdQuadratureMetric
void MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:952
Nektar::StdRegions::StdPyrExp::v_FillMode
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdPyrExp.cpp:771
Nektar::StdRegions::StdPyrExp::GetMode
int GetMode(int I, int J, int K)
Compute the mode number in the expansion for a particular tensorial combination.
Definition: StdPyrExp.cpp:1850
Nektar::StdRegions::StdPyrExp::v_GetEdgeBasisType
virtual LibUtilities::BasisType v_GetEdgeBasisType(const int i) const
Definition: StdPyrExp.cpp:1016
Nektar::StdRegions::EvaluateTriFaceBasisKey
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Definition: StdExpansion3D.cpp:409
Nektar::LibUtilities::BasisKey::GetBasisType
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:144
Nektar::StdRegions::StdPyrExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Definition: StdPyrExp.cpp:1815
Nektar::LibUtilities::eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
Nektar::StdRegions::StdPyrExp::v_GetEdgeInteriorMap
virtual void v_GetEdgeInteriorMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
Definition: StdPyrExp.cpp:1431
std
STL namespace.
Nektar::StdRegions::StdPyrExp::v_GetFaceIntNcoeffs
virtual int v_GetFaceIntNcoeffs(const int i) const
Definition: StdPyrExp.cpp:921
Nektar::LibUtilities::StdPyrData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:267
Nektar::StdRegions::StdPyrExp::v_GetFaceNumPoints
virtual int v_GetFaceNumPoints(const int i) const
Definition: StdPyrExp.cpp:943
Nektar::StdRegions::StdPyrExp::v_GetEdgeNcoeffs
virtual int v_GetEdgeNcoeffs(const int i) const
Definition: StdPyrExp.cpp:883
Nektar::StdRegions::kSVVDGFiltermodesmin
const int kSVVDGFiltermodesmin
Definition: StdRegions.hpp:384
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
Nektar::StdRegions::StdPyrExp::v_FwdTrans
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Forward transform from physical quadrature space stored in inarray and evaluate the expansion coeffic...
Definition: StdPyrExp.cpp:394
Nektar::StdRegions::StdPyrExp::v_GetNfaces
virtual int v_GetNfaces() const
Definition: StdPyrExp.cpp:832
Blas::Dgemm
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where A[m x n], B[n x k], C[m x k].
Definition: Blas.hpp:213
Nektar::StdRegions::StdPyrExp::v_IProductWRTDerivBase
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdPyrExp.cpp:571
class_topology.P
P
Definition: class_topology.py:290
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:714
Nektar::StdRegions::StdPyrExp::v_PhysDeriv
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Calculate the derivative of the physical points.
Definition: StdPyrExp.cpp:106
Nektar::StdRegions::StdPyrExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdPyrExp.cpp:1937
Nektar::StdRegions::StdExpansion::GetFaceNcoeffs
int GetFaceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th face.
Definition: StdExpansion.h:353
Nektar::StdRegions::StdPyrExp::v_CalcNumberOfCoefficients
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdPyrExp.cpp:1003
Nektar::StdRegions::StdPyrExp::v_IProductWRTBase_SumFac
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: StdPyrExp.cpp:446
Nektar::StdRegions::StdExpansion3D::PhysTensorDeriv
void PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d1, Array< OneD, NekDouble > &outarray_d2, Array< OneD, NekDouble > &outarray_d3)
Calculate the 3D derivative in the local tensor/collapsed coordinate at the physical points...
Definition: StdExpansion3D.cpp:69
Nektar::StdRegions::StdExpansion
The base class for all shapes.
Definition: StdExpansion.h:68
Nektar::StdRegions::StdExpansion3D::IProductWRTBase_SumFacKernel
void 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)
Definition: StdExpansion3D.cpp:128
Nektar::StdRegions::StdPyrExp::v_GetBoundaryMap
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdPyrExp.cpp:1750
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:50
Nektar::LibUtilities::BasisKey::GetNumModes
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:86
Nektar::StdRegions::StdPyrExp::v_GetNedges
virtual int v_GetNedges() const
Definition: StdPyrExp.cpp:827
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
Nektar::StdRegions::EvaluateQuadFaceBasisKey
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Definition: StdExpansion3D.cpp:351
Nektar::StdRegions::StdPyrExp::v_LocCoordToLocCollapsed
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Definition: StdPyrExp.cpp:723
ManagerAccess.h
StdPyrExp.h
Nektar::StdRegions::StdPyrExp::v_BwdTrans_SumFac
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdPyrExp.cpp:273
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:161
Nektar::StdRegions::StdPyrExp::v_NumBndryCoeffs
virtual int v_NumBndryCoeffs() const
Definition: StdPyrExp.cpp:842
Nektar::StdRegions::StdPyrExp::v_GetVertexMap
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdPyrExp.cpp:1370
Nektar::StdRegions::StdPyrExp::v_NumDGBndryCoeffs
virtual int v_NumDGBndryCoeffs() const
Definition: StdPyrExp.cpp:862
Nektar::StdRegions::StdExpansion::CreateGeneralMatrix
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
Definition: StdExpansion.cpp:324
Nektar::LibUtilities::eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:45
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:59
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Blas::Dgemv
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = A x where A[m x n].
Definition: Blas.hpp:168
Nektar::StdRegions::StdPyrExp::v_DetFaceBasisKey
virtual const LibUtilities::BasisKey v_DetFaceBasisKey(const int i, const int k) const
Definition: StdPyrExp.cpp:965
Nektar::StdRegions::StdPyrExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdPyrExp.cpp:1710
Nektar::StdRegions::StdPyrExp::v_IProductWRTBase_SumFacKernel
virtual void 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)
Definition: StdPyrExp.cpp:481
Nektar::StdRegions::StdPyrExp::v_ReduceOrderCoeffs
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdPyrExp.cpp:2072
Nektar::StdRegions::StdPyrExp::v_GetFaceNumModes
virtual void v_GetFaceNumModes(const int fid, const Orientation faceOrient, int &numModes0, int &numModes1)
Definition: StdPyrExp.cpp:778
Nektar::StdRegions::StdPyrExp::v_BwdTrans_SumFacKernel
virtual void 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)
Definition: StdPyrExp.cpp:293
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:121
Nektar::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:130
Nektar::StdRegions::StdPyrExp::v_GetFaceToElementMap
virtual void v_GetFaceToElementMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int nummodesA=-1, int nummodesB=-1)
Definition: StdPyrExp.cpp:1038
Nektar::StdRegions::kSVVDGFilter
const NekDouble kSVVDGFilter[9][11]
Definition: StdRegions.hpp:387
Nektar::StdRegions::StdPyrExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdPyrExp.cpp:1820
Nektar::LibUtilities::eOrthoPyr_C
Principle Orthogonal Functions .
Definition: BasisType.h:51
Nektar::LibUtilities::StdPyrData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:240
Blas::Dscal
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:125
Blas::Ddot
static double Ddot(const int &n, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: output = .
Definition: Blas.hpp:140
Nektar::StdRegions::StdExpansion::GetPointsType
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:215
Nektar::StdRegions::StdPyrExp::v_DetShapeType
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdPyrExp.cpp:837
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:130
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:274
Nektar::StdRegions::StdExpansion::GetBasisType
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:164
Nektar::LibUtilities::eGaussRadauMAlpha2Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:59
Nektar::StdRegions::StdExpansion::BwdTrans
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space...
Definition: StdExpansion.h:530
Nektar::StdRegions::StdPyrExp::v_IProductWRTBase
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Inner product of inarray over region with respect to the expansion basis m_base[0]->GetBdata(),m_base[1]->GetBdata(), m_base[2]->GetBdata() and return in outarray.
Definition: StdPyrExp.cpp:430
Nektar::LibUtilities::NullBasisKey
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
Nektar::StdRegions::StdPyrExp::v_StdPhysDeriv
virtual void v_StdPhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Definition: StdPyrExp.cpp:212
Nektar::LibUtilities::eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:54
Nektar::StdRegions::StdPyrExp::v_IProductWRTDerivBase_SumFac
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdPyrExp.cpp:585
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1429
Blas::Daxpy
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
Definition: Blas.hpp:110
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Nektar::StdRegions::StdPyrExp::v_BwdTrans
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Backward transformation is evaluated at the quadrature points.
Definition: StdPyrExp.cpp:252
Nektar::StdRegions::StdExpansion::GetEdgeBasisType
LibUtilities::BasisType GetEdgeBasisType(const int i) const
This function returns the type of expansion basis on the i-th edge.
Definition: StdExpansion.h:412
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:49
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186
Nektar::StdRegions::StdPyrExp::v_GetNverts
virtual int v_GetNverts() const
Definition: StdPyrExp.cpp:822
Nektar::StdRegions::StdPyrExpSharedPtr
std::shared_ptr< StdPyrExp > StdPyrExpSharedPtr
Definition: StdPyrExp.h:270
Nektar::StdRegions::StdExpansion::FwdTrans
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space...
Definition: StdExpansion.h:1957
Generated on Wed Jan 20 2021 03:33:42 for Nektar++ by   doxygen 1.8.13