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.size() > 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.size() > 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.size() > 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  /** \brief Forward transform from physical quadrature space
381  stored in \a inarray and evaluate the expansion coefficients and
382  store in \a outarray
383 
384  Inputs:\n
385 
386  - \a inarray: array of physical quadrature points to be transformed
387 
388  Outputs:\n
389 
390  - \a outarray: updated array of expansion coefficients.
391 
392  */
393  void StdPyrExp::v_FwdTrans(const Array<OneD, const NekDouble> &inarray,
394  Array<OneD, NekDouble> &outarray)
395  {
396  v_IProductWRTBase(inarray,outarray);
397 
398  // get Mass matrix inverse
399  StdMatrixKey imasskey(eInvMass,DetShapeType(),*this);
400  DNekMatSharedPtr imatsys = GetStdMatrix(imasskey);
401 
402 
403  // copy inarray in case inarray == outarray
404  DNekVec in (m_ncoeffs, outarray);
405  DNekVec out(m_ncoeffs, outarray, eWrapper);
406 
407  out = (*imatsys)*in;
408  }
409 
410 
411  //---------------------------------------
412  // Inner product functions
413  //---------------------------------------
414 
415  /** \brief Inner product of \a inarray over region with respect to the
416  expansion basis m_base[0]->GetBdata(),m_base[1]->GetBdata(), m_base[2]->GetBdata() and return in \a outarray
417 
418  Wrapper call to StdPyrExp::IProductWRTBase
419 
420  Input:\n
421 
422  - \a inarray: array of function evaluated at the physical collocation points
423 
424  Output:\n
425 
426  - \a outarray: array of inner product with respect to each basis over region
427 
428  */
429  void StdPyrExp::v_IProductWRTBase(
430  const Array<OneD, const NekDouble> &inarray,
431  Array<OneD, NekDouble> &outarray)
432  {
433  if (m_base[0]->Collocation() &&
434  m_base[1]->Collocation() &&
435  m_base[2]->Collocation())
436  {
437  v_MultiplyByStdQuadratureMetric(inarray, outarray);
438  }
439  else
440  {
441  StdPyrExp::v_IProductWRTBase_SumFac(inarray,outarray);
442  }
443  }
444 
445  void StdPyrExp::v_IProductWRTBase_SumFac(
446  const Array<OneD, const NekDouble>& inarray,
447  Array<OneD, NekDouble>& outarray,
448  bool multiplybyweights)
449  {
450 
451  int nquad1 = m_base[1]->GetNumPoints();
452  int nquad2 = m_base[2]->GetNumPoints();
453  int order0 = m_base[0]->GetNumModes();
454  int order1 = m_base[1]->GetNumModes();
455 
456  Array<OneD, NekDouble> wsp(order0*nquad2*(nquad1 + order1));
457 
458  if(multiplybyweights)
459  {
460  Array<OneD, NekDouble> tmp(inarray.size());
461 
462  v_MultiplyByStdQuadratureMetric(inarray, tmp);
463 
464  v_IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
465  m_base[1]->GetBdata(),
466  m_base[2]->GetBdata(),
467  tmp,outarray,wsp,
468  true,true,true);
469  }
470  else
471  {
472  v_IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
473  m_base[1]->GetBdata(),
474  m_base[2]->GetBdata(),
475  inarray,outarray,wsp,
476  true,true,true);
477  }
478  }
479 
480  void StdPyrExp::v_IProductWRTBase_SumFacKernel(
481  const Array<OneD, const NekDouble> &base0,
482  const Array<OneD, const NekDouble> &base1,
483  const Array<OneD, const NekDouble> &base2,
484  const Array<OneD, const NekDouble> &inarray,
485  Array<OneD, NekDouble> &outarray,
486  Array<OneD, NekDouble> &wsp,
487  bool doCheckCollDir0,
488  bool doCheckCollDir1,
489  bool doCheckCollDir2)
490  {
491  boost::ignore_unused(doCheckCollDir0, doCheckCollDir1,
492  doCheckCollDir2);
493 
494  int nquad0 = m_base[0]->GetNumPoints();
495  int nquad1 = m_base[1]->GetNumPoints();
496  int nquad2 = m_base[2]->GetNumPoints();
497 
498  int order0 = m_base[0]->GetNumModes();
499  int order1 = m_base[1]->GetNumModes();
500  int order2 = m_base[2]->GetNumModes();
501 
502  ASSERTL1(wsp.size() >= nquad1*nquad2*order0 +
503  nquad2*order0*order1,
504  "Insufficient workspace size");
505 
506  Array<OneD, NekDouble > tmp1 = wsp;
507  Array<OneD, NekDouble > tmp2 = wsp + nquad1*nquad2*order0;
508 
509  int i,j, mode,mode1, cnt;
510 
511  // Inner product with respect to the '0' direction
512  Blas::Dgemm('T', 'N', nquad1*nquad2, order0, nquad0,
513  1.0, inarray.get(), nquad0,
514  base0.get(), nquad0,
515  0.0, tmp1.get(), nquad1*nquad2);
516 
517  // Inner product with respect to the '1' direction
518  for(mode=i=0; i < order0; ++i)
519  {
520  Blas::Dgemm('T', 'N', nquad2, order1, nquad1,
521  1.0, tmp1.get()+i*nquad1*nquad2, nquad1,
522  base1.get(), nquad1,
523  0.0, tmp2.get()+mode*nquad2, nquad2);
524  mode += order1;
525  }
526 
527 
528  // Inner product with respect to the '2' direction
529  mode = mode1 = cnt = 0;
530  for(i = 0; i < order0; ++i)
531  {
532  for(j = 0; j < order1; ++j, ++cnt)
533  {
534  int ijmax = max(i,j);
535 
536  Blas::Dgemv('T', nquad2, order2-ijmax,
537  1.0, base2.get()+mode*nquad2, nquad2,
538  tmp2.get()+cnt*nquad2, 1,
539  0.0, outarray.get()+mode1, 1);
540  mode += order2-ijmax;
541  mode1 += order2-ijmax;
542  }
543 
544  //increment mode in case order1!=order2
545  for(j = order1; j < order2; ++j)
546  {
547  int ijmax = max(i,j);
548  mode += order2-ijmax;
549  }
550  }
551 
552  // fix for modified basis for top singular vertex component
553  // Already have evaluated (1+c)/2 (1-b)/2 (1-a)/2
554  if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
555  {
556  // add in (1+c)/2 (1+b)/2 (1-a)/2 component
557  outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
558  &tmp2[nquad2],1);
559 
560  // add in (1+c)/2 (1-b)/2 (1+a)/2 component
561  outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
562  &tmp2[nquad2*order1],1);
563 
564  // add in (1+c)/2 (1+b)/2 (1+a)/2 component
565  outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
566  &tmp2[nquad2*order1+nquad2],1);
567  }
568  }
569 
570  void StdPyrExp::v_IProductWRTDerivBase(
571  const int dir,
572  const Array<OneD, const NekDouble>& inarray,
573  Array<OneD, NekDouble>& outarray)
574  {
575  StdPyrExp::v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
576  }
577 
578  /**
579  * @param inarray Function evaluated at physical collocation
580  * points.
581  * @param outarray Inner product with respect to each basis
582  * function over the element.
583  */
584  void StdPyrExp::v_IProductWRTDerivBase_SumFac(
585  const int dir,
586  const Array<OneD, const NekDouble>& inarray,
587  Array<OneD, NekDouble>& outarray)
588  {
589  int i;
590  int nquad0 = m_base[0]->GetNumPoints();
591  int nquad1 = m_base[1]->GetNumPoints();
592  int nquad2 = m_base[2]->GetNumPoints();
593  int nqtot = nquad0*nquad1*nquad2;
594 
595  Array<OneD, NekDouble> gfac0(nquad0);
596  Array<OneD, NekDouble> gfac1(nquad1);
597  Array<OneD, NekDouble> gfac2(nquad2);
598  Array<OneD, NekDouble> tmp0 (nqtot);
599 
600 
601  int order0 = m_base[0]->GetNumModes();
602  int order1 = m_base[1]->GetNumModes();
603 
604  Array<OneD, NekDouble> wsp(nquad1*nquad2*order0 +
605  nquad2*order0*order1);
606 
607  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
608  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
609  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
610 
611  // set up geometric factor: (1+z0)/2
612  for(i = 0; i < nquad0; ++i)
613  {
614  gfac0[i] = 0.5*(1+z0[i]);
615  }
616 
617  // set up geometric factor: (1+z1)/2
618  for(i = 0; i < nquad1; ++i)
619  {
620  gfac1[i] = 0.5*(1+z1[i]);
621  }
622 
623  // Set up geometric factor: 2/(1-z2)
624  for(i = 0; i < nquad2; ++i)
625  {
626  gfac2[i] = 2.0/(1-z2[i]);
627  }
628 
629  // Derivative in first/second direction is always scaled as follows
630  const int nq01 = nquad0*nquad1;
631  for(i = 0; i < nquad2; ++i)
632  {
633  Vmath::Smul(nq01, gfac2[i],
634  &inarray[0] + i*nq01, 1,
635  &tmp0 [0] + i*nq01, 1);
636  }
637 
638  v_MultiplyByStdQuadratureMetric(tmp0, tmp0);
639 
640  switch(dir)
641  {
642  case 0:
643  {
644  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
645  m_base[1]->GetBdata (),
646  m_base[2]->GetBdata (),
647  tmp0, outarray, wsp,
648  false, true, true);
649  break;
650  }
651  case 1:
652  {
653  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata (),
654  m_base[1]->GetDbdata(),
655  m_base[2]->GetBdata (),
656  tmp0, outarray, wsp,
657  true, false, true);
658  break;
659  }
660  case 2:
661  {
662  Array<OneD, NekDouble> tmp3(m_ncoeffs);
663  Array<OneD, NekDouble> tmp4(m_ncoeffs);
664 
665  // Scale eta_1 derivative by gfac0
666  for(i = 0; i < nquad1*nquad2; ++i)
667  {
668  Vmath::Vmul(nquad0, tmp0 .get() + i*nquad0, 1,
669  gfac0.get(), 1,
670  tmp0 .get() + i*nquad0, 1);
671  }
672  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
673  m_base[1]->GetBdata(),
674  m_base[2]->GetBdata(),
675  tmp0, tmp3, wsp,
676  false, true, true);
677 
678  // Scale eta_2 derivative by gfac1*gfac2
679  for(i = 0; i < nquad2; ++i)
680  {
681  Vmath::Smul(nq01, gfac2[i],
682  &inarray[0] + i*nq01, 1,
683  &tmp0 [0] + i*nq01, 1);
684  }
685  for(i = 0; i < nquad1*nquad2; ++i)
686  {
687  Vmath::Smul(nquad0, gfac1[i%nquad1],
688  &tmp0[0] + i*nquad0, 1,
689  &tmp0[0] + i*nquad0, 1);
690  }
691 
692  v_MultiplyByStdQuadratureMetric(tmp0, tmp0);
693  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
694  m_base[1]->GetDbdata(),
695  m_base[2]->GetBdata(),
696  tmp0, tmp4, wsp,
697  true, false, true);
698 
699  v_MultiplyByStdQuadratureMetric(inarray,tmp0);
700  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
701  m_base[1]->GetBdata(),
702  m_base[2]->GetDbdata(),
703  tmp0,outarray,wsp,
704  true, true, false);
705 
706  Vmath::Vadd(m_ncoeffs,&tmp3[0],1,&outarray[0],1,&outarray[0],1);
707  Vmath::Vadd(m_ncoeffs,&tmp4[0],1,&outarray[0],1,&outarray[0],1);
708  break;
709  }
710  default:
711  {
712  ASSERTL1(false, "input dir is out of range");
713  break;
714  }
715  }
716  }
717 
718  //---------------------------------------
719  // Evaluation functions
720  //---------------------------------------
721 
722  void StdPyrExp::v_LocCoordToLocCollapsed(
723  const Array<OneD, const NekDouble>& xi,
724  Array<OneD, NekDouble>& eta)
725  {
726  NekDouble d2 = 1.0 - xi[2];
727  if(fabs(d2) < NekConstants::kNekZeroTol)
728  {
729  if(d2>=0.)
730  {
731  d2 = NekConstants::kNekZeroTol;
732  }
733  else
734  {
735  d2 = -NekConstants::kNekZeroTol;
736  }
737  }
738  eta[2] = xi[2]; // eta_z = xi_z
739  eta[1] = 2.0*(1.0 + xi[1])/d2 - 1.0;
740  eta[0] = 2.0*(1.0 + xi[0])/d2 - 1.0;
741  }
742 
743  void StdPyrExp::v_LocCollapsedToLocCoord(
744  const Array<OneD, const NekDouble>& eta,
745  Array<OneD, NekDouble>& xi)
746  {
747  xi[0] = (1.0 + eta[0]) * (1.0 - eta[2]) * 0.5 - 1.0;
748  xi[1] = (1.0 + eta[1]) * (1.0 - eta[2]) * 0.5 - 1.0;
749  xi[2] = eta[2];
750  }
751 
752  void StdPyrExp::v_GetCoords(Array<OneD, NekDouble> &xi_x,
753  Array<OneD, NekDouble> &xi_y,
754  Array<OneD, NekDouble> &xi_z)
755  {
756  Array<OneD, const NekDouble> etaBar_x = m_base[0]->GetZ();
757  Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
758  Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
759  int Qx = GetNumPoints(0);
760  int Qy = GetNumPoints(1);
761  int Qz = GetNumPoints(2);
762 
763  // Convert collapsed coordinates into cartesian coordinates: eta --> xi
764  for (int k = 0; k < Qz; ++k )
765  {
766  for (int j = 0; j < Qy; ++j)
767  {
768  for (int i = 0; i < Qx; ++i)
769  {
770  int s = i + Qx*(j + Qy*k);
771 
772  xi_z[s] = eta_z[k];
773  xi_y[s] = (1.0 + eta_y[j]) * (1.0 - eta_z[k]) / 2.0 - 1.0;
774  xi_x[s] = (1.0 + etaBar_x[i]) * (1.0 - eta_z[k]) / 2.0 - 1.0;
775  }
776  }
777  }
778  }
779 
780  void StdPyrExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
781  {
782  Array<OneD, NekDouble> tmp(m_ncoeffs, 0.0);
783  tmp[mode] = 1.0;
784  v_BwdTrans(tmp, outarray);
785  }
786 
787  void StdPyrExp::v_GetTraceNumModes(
788  const int fid,
789  int &numModes0,
790  int &numModes1,
791  Orientation faceOrient)
792  {
793  int nummodes [3] = {m_base[0]->GetNumModes(),
794  m_base[1]->GetNumModes(),
795  m_base[2]->GetNumModes()};
796  switch(fid)
797  {
798  // quad
799  case 0:
800  {
801  numModes0 = nummodes[0];
802  numModes1 = nummodes[1];
803  }
804  break;
805  case 1:
806  case 3:
807  {
808  numModes0 = nummodes[0];
809  numModes1 = nummodes[2];
810  }
811  break;
812  case 2:
813  case 4:
814  {
815  numModes0 = nummodes[1];
816  numModes1 = nummodes[2];
817  }
818  break;
819  }
820 
821  if ( faceOrient >= 9 )
822  {
823  std::swap(numModes0, numModes1);
824  }
825  }
826 
827  NekDouble StdPyrExp::v_PhysEvaluateBasis(
828  const Array<OneD, const NekDouble>& coords,
829  int mode)
830  {
831  Array<OneD, NekDouble> coll(3);
832  LocCoordToLocCollapsed(coords, coll);
833 
834  const int nm0 = m_base[0]->GetNumModes();
835  const int nm1 = m_base[1]->GetNumModes();
836  const int nm2 = m_base[2]->GetNumModes();
837 
838  int mode0 = 0, mode1 = 0, mode2 = 0, cnt = 0;
839 
840  bool found = false;
841  for (mode0 = 0; mode0 < nm0; ++mode0)
842  {
843  for (mode1 = 0; mode1 < nm1; ++mode1)
844  {
845  int maxpq = max(mode0, mode1);
846  for (mode2 = 0; mode2 < nm2 - maxpq; ++mode2, ++cnt)
847  {
848  if (cnt == mode)
849  {
850  found = true;
851  break;
852  }
853  }
854 
855  if (found)
856  {
857  break;
858  }
859  }
860 
861  if (found)
862  {
863  break;
864  }
865 
866  for (int j = nm1; j < nm2; ++j)
867  {
868  int ijmax = max(mode0, j);
869  mode2 += nm2 - ijmax;
870  }
871  }
872 
873  if (mode == 1 &&
874  m_base[0]->GetBasisType() == LibUtilities::eModified_A)
875  {
876  return StdExpansion::BaryEvaluateBasis<2>(coll[2], 1);
877  }
878  else
879  {
880  return
881  StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
882  StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
883  StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
884  }
885  }
886 
887  //---------------------------------------
888  // Helper functions
889  //---------------------------------------
890 
891  int StdPyrExp::v_GetNverts() const
892  {
893  return 5;
894  }
895 
896  int StdPyrExp::v_GetNedges() const
897  {
898  return 8;
899  }
900 
901  int StdPyrExp::v_GetNtraces() const
902  {
903  return 5;
904  }
905 
906  LibUtilities::ShapeType StdPyrExp::v_DetShapeType() const
907  {
908  return LibUtilities::ePyramid;
909  }
910 
911  int StdPyrExp::v_NumBndryCoeffs() const
912  {
913  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
914  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
915  "BasisType is not a boundary interior form");
916  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
917  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
918  "BasisType is not a boundary interior form");
919  ASSERTL1(GetBasisType(2) == LibUtilities::eModifiedPyr_C ||
920  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
921  "BasisType is not a boundary interior form");
922 
923  int P = m_base[0]->GetNumModes();
924  int Q = m_base[1]->GetNumModes();
925  int R = m_base[2]->GetNumModes();
926 
927  return LibUtilities::StdPyrData::
928  getNumberOfBndCoefficients(P, Q, R);
929  }
930 
931  int StdPyrExp::v_NumDGBndryCoeffs() const
932  {
933  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
934  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
935  "BasisType is not a boundary interior form");
936  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
937  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
938  "BasisType is not a boundary interior form");
939  ASSERTL1(GetBasisType(2) == LibUtilities::eModifiedPyr_C ||
940  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
941  "BasisType is not a boundary interior form");
942 
943  int P = m_base[0]->GetNumModes()-1;
944  int Q = m_base[1]->GetNumModes()-1;
945  int R = m_base[2]->GetNumModes()-1;
946 
947  return (P+1)*(Q+1) // 1 rect. face on base
948  + 2*(R+1) + P*(1 + 2*R - P) // 2 tri. (P,R) faces
949  + 2*(R+1) + Q*(1 + 2*R - Q); // 2 tri. (Q,R) faces
950  }
951 
952 
953  int StdPyrExp::v_GetTraceNcoeffs(const int i) const
954  {
955  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
956 
957  if (i == 0)
958  {
959  return GetBasisNumModes(0)*GetBasisNumModes(1);
960  }
961  else if (i == 1 || i == 3)
962  {
963  int P = GetBasisNumModes(0)-1, Q = GetBasisNumModes(2)-1;
964  return Q+1 + (P*(1 + 2*Q - P))/2;
965  }
966  else
967  {
968  int P = GetBasisNumModes(1)-1, Q = GetBasisNumModes(2)-1;
969  return Q+1 + (P*(1 + 2*Q - P))/2;
970  }
971  }
972 
973  int StdPyrExp::v_GetTraceIntNcoeffs(const int i) const
974  {
975  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
976 
977  int P = m_base[0]->GetNumModes()-1;
978  int Q = m_base[1]->GetNumModes()-1;
979  int R = m_base[2]->GetNumModes()-1;
980 
981  if (i == 0)
982  {
983  return (P-1)*(Q-1);
984  }
985  else if (i == 1 || i == 3)
986  {
987  return (P-1) * (2*(R-1) - (P-1) - 1) / 2;
988  }
989  else
990  {
991  return (Q-1) * (2*(R-1) - (Q-1) - 1) / 2;
992  }
993  }
994 
995  int StdPyrExp::v_GetTraceNumPoints(const int i) const
996  {
997  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
998 
999  if (i == 0)
1000  {
1001  return m_base[0]->GetNumPoints()*
1002  m_base[1]->GetNumPoints();
1003  }
1004  else if (i == 1 || i == 3)
1005  {
1006  return m_base[0]->GetNumPoints()*
1007  m_base[2]->GetNumPoints();
1008  }
1009  else
1010  {
1011  return m_base[1]->GetNumPoints()*
1012  m_base[2]->GetNumPoints();
1013  }
1014  }
1015 
1016  int StdPyrExp::v_GetEdgeNcoeffs(const int i) const
1017  {
1018  ASSERTL2(i >= 0 && i <= 7, "edge id is out of range");
1019 
1020  if (i == 0 || i == 2)
1021  {
1022  return GetBasisNumModes(0);
1023  }
1024  else if (i == 1 || i == 3)
1025  {
1026  return GetBasisNumModes(1);
1027  }
1028  else
1029  {
1030  return GetBasisNumModes(2);
1031  }
1032  }
1033 
1034  const LibUtilities::BasisKey StdPyrExp::v_GetTraceBasisKey(
1035  const int i, const int k) const
1036  {
1037  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
1038  ASSERTL2(k >= 0 && k <= 1, "basis key id is out of range");
1039 
1040  switch(i)
1041  {
1042  case 0:
1043  {
1044  return EvaluateQuadFaceBasisKey(k,
1045  m_base[k]->GetBasisType(),
1046  m_base[k]->GetNumPoints(),
1047  m_base[k]->GetNumModes());
1048 
1049  }
1050  case 1:
1051  case 3:
1052  {
1053  return EvaluateTriFaceBasisKey(k,
1054  m_base[2*k]->GetBasisType(),
1055  m_base[2*k]->GetNumPoints(),
1056  m_base[2*k]->GetNumModes());
1057  }
1058  case 2:
1059  case 4:
1060  {
1061  return EvaluateTriFaceBasisKey(k,
1062  m_base[k+1]->GetBasisType(),
1063  m_base[k+1]->GetNumPoints(),
1064  m_base[k+1]->GetNumModes());
1065  }
1066  }
1067 
1068  // Should never get here.
1069  return LibUtilities::NullBasisKey;
1070  }
1071 
1072  int StdPyrExp::v_CalcNumberOfCoefficients(
1073  const std::vector<unsigned int> &nummodes,
1074  int &modes_offset)
1075  {
1076  int nmodes = LibUtilities::StdPyrData::getNumberOfCoefficients(
1077  nummodes[modes_offset],
1078  nummodes[modes_offset+1],
1079  nummodes[modes_offset+2]);
1080 
1081  modes_offset += 3;
1082  return nmodes;
1083  }
1084 
1085  int StdPyrExp::v_GetVertexMap(int vId, bool useCoeffPacking)
1086  {
1087  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1088  GetBasisType(1) == LibUtilities::eModified_A ||
1089  GetBasisType(2) == LibUtilities::eModifiedPyr_C,
1090  "Mapping not defined for this type of basis");
1091 
1092  int l = 0;
1093 
1094  if(useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1095  {
1096  switch (vId)
1097  {
1098  case 0:
1099  l = GetMode(0,0,0);
1100  break;
1101  case 1:
1102  l = GetMode(0,0,1);
1103  break;
1104  case 2:
1105  l = GetMode(0,1,0);
1106  break;
1107  case 3:
1108  l = GetMode(1,0,0);
1109  break;
1110  case 4:
1111  l = GetMode(1,1,0);
1112  break;
1113  default:
1114  ASSERTL0(false, "local vertex id must be between 0 and 4");
1115  }
1116  }
1117  else
1118  {
1119  switch (vId)
1120  {
1121  case 0:
1122  l = GetMode(0,0,0);
1123  break;
1124  case 1:
1125  l = GetMode(1,0,0);
1126  break;
1127  case 2:
1128  l = GetMode(1,1,0);
1129  break;
1130  case 3:
1131  l = GetMode(0,1,0);
1132  break;
1133  case 4:
1134  l = GetMode(0,0,1);
1135  break;
1136  default:
1137  ASSERTL0(false, "local vertex id must be between 0 and 4");
1138  }
1139  }
1140 
1141  return l;
1142  }
1143 
1144  void StdPyrExp::v_GetInteriorMap(Array<OneD, unsigned int> &outarray)
1145  {
1146  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1147  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1148  "BasisType is not a boundary interior form");
1149  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1150  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1151  "BasisType is not a boundary interior form");
1152  ASSERTL1(GetBasisType(2) == LibUtilities::eModifiedPyr_C ||
1153  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1154  "BasisType is not a boundary interior form");
1155 
1156 
1157  int P = m_base[0]->GetNumModes() - 1, p;
1158  int Q = m_base[1]->GetNumModes() - 1, q;
1159  int R = m_base[2]->GetNumModes() - 1, r;
1160 
1161  int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1162 
1163  if(outarray.size()!=nIntCoeffs)
1164  {
1165  outarray = Array<OneD, unsigned int>(nIntCoeffs);
1166  }
1167 
1168  int idx = 0;
1169 
1170  // Loop over all interior modes.
1171  for (p = 2; p <= P; ++p)
1172  {
1173  for (q = 2; q <= Q; ++q)
1174  {
1175  int maxpq = max(p,q);
1176  for (r = 1; r <= R-maxpq; ++r)
1177  {
1178  outarray[idx++] = GetMode(p,q,r);
1179  }
1180  }
1181  }
1182  }
1183 
1184  void StdPyrExp::v_GetBoundaryMap(Array<OneD, unsigned int> &maparray)
1185  {
1186  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1187  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1188  "BasisType is not a boundary interior form");
1189  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1190  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1191  "BasisType is not a boundary interior form");
1192  ASSERTL1(GetBasisType(2) == LibUtilities::eModifiedPyr_C ||
1193  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1194  "BasisType is not a boundary interior form");
1195 
1196  int P = m_base[0]->GetNumModes() - 1, p;
1197  int Q = m_base[1]->GetNumModes() - 1, q;
1198  int R = m_base[2]->GetNumModes() - 1, r;
1199  int idx = 0;
1200 
1201  int nBnd = NumBndryCoeffs();
1202 
1203  if (maparray.size() != nBnd)
1204  {
1205  maparray = Array<OneD, unsigned int>(nBnd);
1206  }
1207 
1208  // Loop over all boundary modes (in ascending order).
1209  for (p = 0; p <= P; ++p)
1210  {
1211  // First two q-r planes are entirely boundary modes.
1212  if (p <= 1)
1213  {
1214  for (q = 0; q <= Q; ++q)
1215  {
1216  int maxpq = max(p,q);
1217  for (r = 0; r <= R-maxpq; ++r)
1218  {
1219  maparray[idx++] = GetMode(p,q,r);
1220  }
1221  }
1222  }
1223  else
1224  {
1225  // Remaining q-r planes contain boundary modes on the two
1226  // front and back sides and edges 0 2.
1227  for (q = 0; q <= Q; ++q)
1228  {
1229  if (q <= 1)
1230  {
1231  for (r = 0; r <= R-p; ++r)
1232  {
1233  maparray[idx++] = GetMode(p,q,r);
1234  }
1235  }
1236  else
1237  {
1238  maparray[idx++] = GetMode(p,q,0);
1239  }
1240  }
1241  }
1242  }
1243  }
1244 
1245  void StdPyrExp::v_GetTraceToElementMap(
1246  const int fid,
1247  Array<OneD, unsigned int> &maparray,
1248  Array<OneD, int> &signarray,
1249  const Orientation faceOrient,
1250  int P,
1251  int Q)
1252  {
1253  ASSERTL1(GetBasisType(0) == GetBasisType(1),
1254  "Method only implemented if BasisType is identical"
1255  "in x and y directions");
1256  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A &&
1257  GetBasisType(2) == LibUtilities::eModifiedPyr_C,
1258  "Method only implemented for Modified_A BasisType"
1259  "(x and y direction) and ModifiedPyr_C BasisType (z "
1260  "direction)");
1261 
1262  int i, j, k, p, q, r, nFaceCoeffs, idx = 0;
1263  int nummodesA=0, nummodesB=0;
1264 
1265  int order0 = m_base[0]->GetNumModes();
1266  int order1 = m_base[1]->GetNumModes();
1267  int order2 = m_base[2]->GetNumModes();
1268 
1269  switch (fid)
1270  {
1271  case 0:
1272  nummodesA = order0;
1273  nummodesB = order1;
1274  break;
1275  case 1:
1276  case 3:
1277  nummodesA = order0;
1278  nummodesB = order2;
1279  break;
1280  case 2:
1281  case 4:
1282  nummodesA = order1;
1283  nummodesB = order2;
1284  break;
1285  default:
1286  ASSERTL0(false,"fid must be between 0 and 4");
1287  }
1288 
1289  bool CheckForZeroedModes = false;
1290 
1291  if (P == -1)
1292  {
1293  P = nummodesA;
1294  Q = nummodesB;
1295  nFaceCoeffs = GetTraceNcoeffs(fid);
1296  }
1297  else if (fid > 0)
1298  {
1299  nFaceCoeffs = P*(2*Q-P+1)/2;
1300  CheckForZeroedModes = true;
1301  }
1302  else
1303  {
1304  nFaceCoeffs = P*Q;
1305  CheckForZeroedModes = true;
1306  }
1307 
1308  // Allocate the map array and sign array; set sign array to ones (+)
1309  if (maparray.size() != nFaceCoeffs)
1310  {
1311  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1312  }
1313 
1314  if (signarray.size() != nFaceCoeffs)
1315  {
1316  signarray = Array<OneD, int>(nFaceCoeffs,1);
1317  }
1318  else
1319  {
1320  fill(signarray.get(), signarray.get() + nFaceCoeffs, 1);
1321  }
1322 
1323  // Set up an array indexing for quads, since the ordering may need
1324  // to be transposed.
1325  Array<OneD, int> arrayindx(nFaceCoeffs,-1);
1326 
1327  if (fid == 0)
1328  {
1329  for (i = 0; i < Q; i++)
1330  {
1331  for (j = 0; j < P; j++)
1332  {
1333  if (faceOrient < 9)
1334  {
1335  arrayindx[i*P+j] = i*P+j;
1336  }
1337  else
1338  {
1339  arrayindx[i*P+j] = j*Q+i;
1340  }
1341  }
1342  }
1343  }
1344 
1345  // Set up ordering inside each 2D face. Also for triangular faces,
1346  // populate signarray.
1347  switch (fid)
1348  {
1349  case 0: // Bottom quad
1350 
1351  for (q = 0; q < Q; ++q)
1352  {
1353  for (p = 0; p < P; ++p)
1354  {
1355  maparray[arrayindx[q*P+p]] = GetMode(p,q,0);
1356  }
1357  }
1358  break;
1359 
1360  case 1: // Front triangle
1361  for (p = 0; p < P; ++p)
1362  {
1363  for (r = 0; r < Q-p; ++r)
1364  {
1365  if ((int)faceOrient == 7 && p > 1)
1366  {
1367  signarray[idx] = p % 2 ? -1 : 1;
1368  }
1369  maparray[idx++] = GetMode(p,0,r);
1370  }
1371  }
1372  break;
1373 
1374  case 2: // Right triangle
1375  maparray[idx++] = GetMode(1,0,0);
1376  maparray[idx++] = GetMode(0,0,1);
1377  for (r = 1; r < Q-1; ++r)
1378  {
1379  maparray[idx++] = GetMode(1,0,r);
1380  }
1381 
1382  for (q = 1; q < P; ++q)
1383  {
1384  for (r = 0; r < Q-q; ++r)
1385  {
1386  if ((int)faceOrient == 7 && q > 1)
1387  {
1388  signarray[idx] = q % 2 ? -1 : 1;
1389  }
1390  maparray[idx++] = GetMode(1,q,r);
1391  }
1392  }
1393  break;
1394 
1395  case 3: // Rear triangle
1396  maparray[idx++] = GetMode(0,1,0);
1397  maparray[idx++] = GetMode(0,0,1);
1398  for (r = 1; r < Q-1; ++r)
1399  {
1400  maparray[idx++] = GetMode(0,1,r);
1401  }
1402 
1403  for (p = 1; p < P; ++p)
1404  {
1405  for (r = 0; r < Q-p; ++r)
1406  {
1407  if ((int)faceOrient == 7 && p > 1)
1408  {
1409  signarray[idx] = p % 2 ? -1 : 1;
1410  }
1411  maparray[idx++] = GetMode(p, 1, r);
1412  }
1413  }
1414  break;
1415 
1416  case 4: // Left triangle
1417  for (q = 0; q < P; ++q)
1418  {
1419  for (r = 0; r < Q-q; ++r)
1420  {
1421  if ((int)faceOrient == 7 && q > 1)
1422  {
1423  signarray[idx] = q % 2 ? -1 : 1;
1424  }
1425  maparray[idx++] = GetMode(0,q,r);
1426  }
1427  }
1428  break;
1429 
1430  default:
1431  ASSERTL0(false, "Face to element map unavailable.");
1432  }
1433 
1434  if (fid > 0)
1435  {
1436  if(CheckForZeroedModes)
1437  {
1438  // zero signmap and set maparray to zero if elemental
1439  // modes are not as large as face modesl
1440  int idx = 0;
1441  for (j = 0; j < P; ++j)
1442  {
1443  idx += Q-j;
1444  for (k = Q-j; k < Q-j; ++k)
1445  {
1446  signarray[idx] = 0.0;
1447  maparray[idx++] = maparray[0];
1448  }
1449  }
1450 
1451  for (j = P; j < P; ++j)
1452  {
1453  for (k = 0; k < Q-j; ++k)
1454  {
1455  signarray[idx] = 0.0;
1456  maparray[idx++] = maparray[0];
1457  }
1458  }
1459  }
1460 
1461  // Triangles only have one possible orientation (base
1462  // direction reversed); swap edge modes.
1463  if ((int)faceOrient == 7)
1464  {
1465  swap(maparray[0], maparray[Q]);
1466  for (i = 1; i < Q-1; ++i)
1467  {
1468  swap(maparray[i+1], maparray[Q+i]);
1469  }
1470  }
1471  }
1472  else
1473  {
1474  if(CheckForZeroedModes)
1475  {
1476  // zero signmap and set maparray to zero if elemental
1477  // modes are not as large as face modesl
1478  for (j = 0; j < P; ++j)
1479  {
1480  for (k = Q; k < Q; ++k)
1481  {
1482  signarray[arrayindx[j+k*P]] = 0.0;
1483  maparray[arrayindx[j+k*P]] = maparray[0];
1484  }
1485  }
1486 
1487  for (j = P; j < P; ++j)
1488  {
1489  for (k = 0; k < Q; ++k)
1490  {
1491  signarray[arrayindx[j+k*P]] = 0.0;
1492  maparray[arrayindx[j+k*P]] = maparray[0];
1493  }
1494  }
1495  }
1496 
1497  // The code below is exactly the same as that taken from
1498  // StdHexExp and reverses the 'b' and 'a' directions as
1499  // appropriate (1st and 2nd if statements respectively) in
1500  // quadrilateral faces.
1501  if (faceOrient == 6 || faceOrient == 8 ||
1502  faceOrient == 11 || faceOrient == 12)
1503  {
1504  if (faceOrient < 9)
1505  {
1506  for (i = 3; i < Q; i += 2)
1507  {
1508  for (j = 0; j < P; j++)
1509  {
1510  signarray[arrayindx[i*P+j]] *= -1;
1511  }
1512  }
1513 
1514  for (i = 0; i < P; i++)
1515  {
1516  swap(maparray [i], maparray [i+P]);
1517  swap(signarray[i], signarray[i+P]);
1518  }
1519  }
1520  else
1521  {
1522  for (i = 0; i < Q; i++)
1523  {
1524  for (j = 3; j < P; j += 2)
1525  {
1526  signarray[arrayindx[i*P+j]] *= -1;
1527  }
1528  }
1529 
1530  for (i = 0; i < Q; i++)
1531  {
1532  swap (maparray [i], maparray [i+Q]);
1533  swap (signarray[i], signarray[i+Q]);
1534  }
1535  }
1536  }
1537 
1538  if (faceOrient == 7 || faceOrient == 8 ||
1539  faceOrient == 10 || faceOrient == 12)
1540  {
1541  if (faceOrient < 9)
1542  {
1543  for (i = 0; i < Q; i++)
1544  {
1545  for (j = 3; j < P; j += 2)
1546  {
1547  signarray[arrayindx[i*P+j]] *= -1;
1548  }
1549  }
1550 
1551  for(i = 0; i < Q; i++)
1552  {
1553  swap(maparray [i*P], maparray [i*P+1]);
1554  swap(signarray[i*P], signarray[i*P+1]);
1555  }
1556  }
1557  else
1558  {
1559  for (i = 3; i < Q; i += 2)
1560  {
1561  for (j = 0; j < P; j++)
1562  {
1563  signarray[arrayindx[i*P+j]] *= -1;
1564  }
1565  }
1566 
1567  for (i = 0; i < P; i++)
1568  {
1569  swap(maparray [i*Q], maparray [i*Q+1]);
1570  swap(signarray[i*Q], signarray[i*Q+1]);
1571  }
1572  }
1573  }
1574  }
1575  }
1576 
1577  void StdPyrExp::v_GetEdgeInteriorToElementMap(
1578  const int eid,
1579  Array<OneD, unsigned int> &maparray,
1580  Array<OneD, int> &signarray,
1581  const Orientation edgeOrient)
1582  {
1583  int i;
1584  bool signChange;
1585  const int P = m_base[0]->GetNumModes() - 1;
1586  const int Q = m_base[1]->GetNumModes() - 1;
1587  const int R = m_base[2]->GetNumModes() - 1;
1588  const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid) - 2;
1589 
1590  if (maparray.size() != nEdgeIntCoeffs)
1591  {
1592  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1593  }
1594 
1595  if(signarray.size() != nEdgeIntCoeffs)
1596  {
1597  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1598  }
1599  else
1600  {
1601  fill(signarray.get(), signarray.get()+nEdgeIntCoeffs, 1);
1602  }
1603 
1604  // If edge is oriented backwards, change sign of modes which have
1605  // degree 2n+1, n >= 1.
1606  signChange = edgeOrient == eBackwards;
1607 
1608  switch (eid)
1609  {
1610  case 0:
1611  for (i = 2; i <= P; ++i)
1612  {
1613  maparray[i-2] = GetMode(i,0,0);
1614  }
1615  break;
1616 
1617  case 1:
1618  for (i = 2; i <= Q; ++i)
1619  {
1620  maparray[i-2] = GetMode(1,i,0);
1621  }
1622  break;
1623  case 2:
1624  for (i = 2; i <= P; ++i)
1625  {
1626  maparray[i-2] = GetMode(i,1,0);
1627  }
1628  break;
1629 
1630  case 3:
1631  for (i = 2; i <= Q; ++i)
1632  {
1633  maparray[i-2] = GetMode(0,i,0);
1634  }
1635  break;
1636  case 4:
1637  for (i = 2; i <= R; ++i)
1638  {
1639  maparray[i-2] = GetMode(0,0,i);
1640  }
1641  break;
1642 
1643  case 5:
1644  for (i = 1; i <= R-1; ++i)
1645  {
1646  maparray[i-1] = GetMode(1,0,i);
1647  }
1648  break;
1649  case 6:
1650  for (i = 1; i <= R-1; ++i)
1651  {
1652  maparray[i-1] = GetMode(1,1,i);
1653  }
1654  break;
1655 
1656  case 7:
1657  for (i = 1; i <= R-1; ++i)
1658  {
1659  maparray[i-1] = GetMode(0,1,i);
1660  }
1661  break;
1662  default:
1663  ASSERTL0(false, "Edge not defined.");
1664  break;
1665  }
1666 
1667  if (signChange)
1668  {
1669  for (i = 1; i < nEdgeIntCoeffs; i += 2)
1670  {
1671  signarray[i] = -1;
1672  }
1673  }
1674  }
1675 
1676  void StdPyrExp::v_GetTraceInteriorToElementMap(
1677  const int fid,
1678  Array<OneD, unsigned int> &maparray,
1679  Array<OneD, int> &signarray,
1680  const Orientation faceOrient)
1681  {
1682  const int P = m_base[0]->GetNumModes() - 1;
1683  const int Q = m_base[1]->GetNumModes() - 1;
1684  const int R = m_base[2]->GetNumModes() - 1;
1685  const int nFaceIntCoeffs = v_GetTraceIntNcoeffs(fid);
1686  int p, q, r, idx = 0;
1687  int nummodesA = 0;
1688  int nummodesB = 0;
1689  int i, j;
1690 
1691  if (maparray.size() != nFaceIntCoeffs)
1692  {
1693  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1694  }
1695 
1696  if (signarray.size() != nFaceIntCoeffs)
1697  {
1698  signarray = Array<OneD, int>(nFaceIntCoeffs, 1);
1699  }
1700  else
1701  {
1702  fill(signarray.get(), signarray.get()+nFaceIntCoeffs, 1);
1703  }
1704 
1705  // Set up an array indexing for quad faces, since the ordering may
1706  // need to be transposed depending on orientation.
1707  Array<OneD, int> arrayindx(nFaceIntCoeffs);
1708  if (fid == 0)
1709  {
1710  nummodesA = P-1;
1711  nummodesB = Q-1;
1712 
1713  for (i = 0; i < nummodesB; i++)
1714  {
1715  for (j = 0; j < nummodesA; j++)
1716  {
1717  if (faceOrient < 9)
1718  {
1719  arrayindx[i*nummodesA+j] = i*nummodesA+j;
1720  }
1721  else
1722  {
1723  arrayindx[i*nummodesA+j] = j*nummodesB+i;
1724  }
1725  }
1726  }
1727  }
1728 
1729  switch (fid)
1730  {
1731  case 0: // Bottom quad
1732  for (q = 2; q <= Q; ++q)
1733  {
1734  for (p = 2; p <= P; ++p)
1735  {
1736  maparray[arrayindx[(q-2)*nummodesA+(p-2)]] = GetMode(p,q,0);
1737  }
1738  }
1739  break;
1740  case 1: // Front triangle
1741  for (p = 2; p <= P; ++p)
1742  {
1743  for (r = 1; r <= R-p; ++r)
1744  {
1745  if ((int)faceOrient == 7)
1746  {
1747  signarray[idx] = p % 2 ? -1 : 1;
1748  }
1749  maparray[idx++] = GetMode(p,0,r);
1750  }
1751  }
1752  break;
1753  case 2: // Right triangle
1754  for (q = 2; q <= Q; ++q)
1755  {
1756  for (r = 1; r <= R-q; ++r)
1757  {
1758  if ((int)faceOrient == 7)
1759  {
1760  signarray[idx] = q % 2 ? -1 : 1;
1761  }
1762  maparray[idx++] = GetMode(1, q, r);
1763  }
1764  }
1765  break;
1766 
1767  case 3: // Rear triangle
1768  for (p = 2; p <= P; ++p)
1769  {
1770  for (r = 1; r <= R-p; ++r)
1771  {
1772  if ((int)faceOrient == 7)
1773  {
1774  signarray[idx] = p % 2 ? -1 : 1;
1775  }
1776  maparray[idx++] = GetMode(p, 1, r);
1777  }
1778  }
1779  break;
1780 
1781  case 4: // Left triangle
1782  for (q = 2; q <= Q; ++q)
1783  {
1784  for (r = 1; r <= R-q; ++r)
1785  {
1786  if ((int)faceOrient == 7)
1787  {
1788  signarray[idx] = q % 2 ? -1 : 1;
1789  }
1790  maparray[idx++] = GetMode(0, q, r);
1791  }
1792  }
1793  break;
1794  default:
1795  ASSERTL0(false, "Face interior map not available.");
1796  }
1797 
1798  // Triangular faces are processed in the above switch loop; for
1799  // remaining quad faces, set up orientation if necessary.
1800  if (fid > 0)
1801  {
1802  return;
1803  }
1804 
1805  if (faceOrient == 6 || faceOrient == 8 ||
1806  faceOrient == 11 || faceOrient == 12)
1807  {
1808  if (faceOrient < 9)
1809  {
1810  for (i = 1; i < nummodesB; i += 2)
1811  {
1812  for (j = 0; j < nummodesA; j++)
1813  {
1814  signarray[arrayindx[i*nummodesA+j]] *= -1;
1815  }
1816  }
1817  }
1818  else
1819  {
1820  for (i = 0; i < nummodesB; i++)
1821  {
1822  for (j = 1; j < nummodesA; j += 2)
1823  {
1824  signarray[arrayindx[i*nummodesA+j]] *= -1;
1825  }
1826  }
1827  }
1828  }
1829 
1830  if (faceOrient == 7 || faceOrient == 8 ||
1831  faceOrient == 10 || faceOrient == 12)
1832  {
1833  if (faceOrient < 9)
1834  {
1835  for (i = 0; i < nummodesB; i++)
1836  {
1837  for (j = 1; j < nummodesA; j += 2)
1838  {
1839  signarray[arrayindx[i*nummodesA+j]] *= -1;
1840  }
1841  }
1842  }
1843  else
1844  {
1845  for (i = 1; i < nummodesB; i += 2)
1846  {
1847  for (j = 0; j < nummodesA; j++)
1848  {
1849  signarray[arrayindx[i*nummodesA+j]] *= -1;
1850  }
1851  }
1852  }
1853  }
1854  }
1855 
1856  //---------------------------------------
1857  // Wrapper functions
1858  //---------------------------------------
1859 
1860  DNekMatSharedPtr StdPyrExp::v_GenMatrix(const StdMatrixKey &mkey)
1861  {
1862  return CreateGeneralMatrix(mkey);
1863  }
1864 
1865  DNekMatSharedPtr StdPyrExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
1866  {
1867  return v_GenMatrix(mkey);
1868  }
1869 
1870 
1871  /**
1872  * @brief Compute the mode number in the expansion for a
1873  * particular tensorial combination.
1874  *
1875  * Modes are numbered with the r index travelling fastest,
1876  * followed by q and then p, and each q-r plane is of size
1877  *
1878  * (R+1-p)*(Q+1) - l(l+1)/2 where l = max(0,Q-p)
1879  *
1880  * For example, when P=2, Q=3 and R=4 the indexing inside each
1881  * q-r plane (with r increasing upwards and q to the right)
1882  * is:
1883  *
1884  * p = 0: p = 1: p = 2:
1885  * ----------------------------------
1886  * 4
1887  * 3 8 17 21
1888  * 2 7 11 16 20 24 29 32 35
1889  * 1 6 10 13 15 19 23 26 28 31 34 37
1890  * 0 5 9 12 14 18 22 25 27 30 33 36
1891  *
1892  * Note that in this element, we must have that \f$ P,Q \leq
1893  * R\f$.
1894  */
1895  int StdPyrExp::GetMode(const int I, const int J, const int K)
1896  {
1897  const int Q = m_base[1]->GetNumModes()-1;
1898  const int R = m_base[2]->GetNumModes()-1;
1899 
1900  int i,l;
1901  int cnt = 0;
1902 
1903  // Traverse to q-r plane number I
1904  for (i = 0; i < I; ++i)
1905  {
1906  // Size of triangle part
1907  l = max(0,Q-i);
1908 
1909  // Size of rectangle part
1910  cnt += (R+1-i)*(Q+1) - l*(l+1)/2;
1911  }
1912 
1913  // Traverse to q column J (Pretend this is a face of width J)
1914  l = max(0,J-1-I);
1915  cnt += (R+1-I)*J - l*(l+1)/2;
1916 
1917  // Traverse up stacks to K
1918  cnt += K;
1919 
1920  return cnt;
1921  }
1922 
1923 
1924  void StdPyrExp::v_MultiplyByStdQuadratureMetric(
1925  const Array<OneD, const NekDouble>& inarray,
1926  Array<OneD, NekDouble>& outarray)
1927  {
1928  int i, j;
1929 
1930  int nquad0 = m_base[0]->GetNumPoints();
1931  int nquad1 = m_base[1]->GetNumPoints();
1932  int nquad2 = m_base[2]->GetNumPoints();
1933 
1934  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1935  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1936  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
1937 
1938  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1939 
1940  // Multiply by integration constants in x-direction
1941  for(i = 0; i < nquad1*nquad2; ++i)
1942  {
1943  Vmath::Vmul(nquad0, inarray.get()+i*nquad0, 1,
1944  w0.get(), 1, outarray.get()+i*nquad0,1);
1945  }
1946 
1947  // Multiply by integration constants in y-direction
1948  for(j = 0; j < nquad2; ++j)
1949  {
1950  for(i = 0; i < nquad1; ++i)
1951  {
1952  Blas::Dscal(nquad0,w1[i], &outarray[0]+i*nquad0 +
1953  j*nquad0*nquad1,1);
1954  }
1955  }
1956 
1957  // Multiply by integration constants in z-direction; need to
1958  // incorporate factor [(1-eta_3)/2]^2 into weights, but only if
1959  // using GLL quadrature points.
1960  switch(m_base[2]->GetPointsType())
1961  {
1962  // (2,0) Jacobi inner product.
1963  case LibUtilities::eGaussRadauMAlpha2Beta0:
1964  for(i = 0; i < nquad2; ++i)
1965  {
1966  Blas::Dscal(nquad0*nquad1, 0.25*w2[i],
1967  &outarray[0]+i*nquad0*nquad1, 1);
1968  }
1969  break;
1970 
1971  default:
1972  for(i = 0; i < nquad2; ++i)
1973  {
1974  Blas::Dscal(nquad0*nquad1,0.25*(1-z2[i])*(1-z2[i])*w2[i],
1975  &outarray[0]+i*nquad0*nquad1,1);
1976  }
1977  break;
1978  }
1979  }
1980 
1981 
1982  void StdPyrExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
1983  const StdMatrixKey &mkey)
1984  {
1985  // Generate an orthonogal expansion
1986  int qa = m_base[0]->GetNumPoints();
1987  int qb = m_base[1]->GetNumPoints();
1988  int qc = m_base[2]->GetNumPoints();
1989  int nmodes_a = m_base[0]->GetNumModes();
1990  int nmodes_b = m_base[1]->GetNumModes();
1991  int nmodes_c = m_base[2]->GetNumModes();
1992  // Declare orthogonal basis.
1993  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
1994  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
1995  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
1996 
1997  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
1998  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A,nmodes_b,pb);
1999  LibUtilities::BasisKey Bc(LibUtilities::eOrthoPyr_C,nmodes_c,pc);
2000  StdPyrExp OrthoExp(Ba,Bb,Bc);
2001 
2002  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2003  int i,j,k,cnt = 0;
2004 
2005  // project onto modal space.
2006  OrthoExp.FwdTrans(array,orthocoeffs);
2007 
2008  if(mkey.ConstFactorExists(eFactorSVVPowerKerDiffCoeff))
2009  {
2010  // Rodrigo's power kernel
2011  NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
2012  NekDouble SvvDiffCoeff =
2013  mkey.GetConstFactor(eFactorSVVPowerKerDiffCoeff)*
2014  mkey.GetConstFactor(eFactorSVVDiffCoeff);
2015 
2016  for(int i = 0; i < nmodes_a; ++i)
2017  {
2018  for(int j = 0; j < nmodes_b; ++j)
2019  {
2020  int maxij = max(i,j);
2021  NekDouble fac1 = std::max(
2022  pow((1.0*i)/(nmodes_a-1),cutoff*nmodes_a),
2023  pow((1.0*j)/(nmodes_b-1),cutoff*nmodes_b));
2024 
2025  for(int k = 0; k < nmodes_c-maxij; ++k)
2026  {
2027  NekDouble fac = std::max(fac1,
2028  pow((1.0*k)/(nmodes_c-1),cutoff*nmodes_c));
2029 
2030  orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2031  cnt++;
2032  }
2033  }
2034  }
2035  }
2036  else if(mkey.ConstFactorExists(eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
2037  {
2038  NekDouble SvvDiffCoeff =
2039  mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff)*
2040  mkey.GetConstFactor(eFactorSVVDiffCoeff);
2041 
2042  int max_abc = max(nmodes_a-kSVVDGFiltermodesmin,
2043  nmodes_b-kSVVDGFiltermodesmin);
2044  max_abc = max(max_abc, nmodes_c-kSVVDGFiltermodesmin);
2045  // clamp max_abc
2046  max_abc = max(max_abc,0);
2047  max_abc = min(max_abc,kSVVDGFiltermodesmax-kSVVDGFiltermodesmin);
2048 
2049  for(int i = 0; i < nmodes_a; ++i)
2050  {
2051  for(int j = 0; j < nmodes_b; ++j)
2052  {
2053  int maxij = max(i,j);
2054 
2055  for(int k = 0; k < nmodes_c-maxij; ++k)
2056  {
2057  int maxijk = max(maxij,k);
2058  maxijk = min(maxijk,kSVVDGFiltermodesmax-1);
2059 
2060  orthocoeffs[cnt] *= SvvDiffCoeff *
2061  kSVVDGFilter[max_abc][maxijk];
2062  cnt++;
2063  }
2064  }
2065  }
2066  }
2067  else
2068  {
2069  //SVV filter paramaters (how much added diffusion relative
2070  // to physical one and fraction of modes from which you
2071  // start applying this added diffusion)
2072  //
2073  NekDouble SvvDiffCoeff = mkey.GetConstFactor(StdRegions::eFactorSVVDiffCoeff);
2074  NekDouble SVVCutOff = mkey.GetConstFactor(StdRegions::eFactorSVVCutoffRatio);
2075 
2076  //Defining the cut of mode
2077  int cutoff_a = (int) (SVVCutOff*nmodes_a);
2078  int cutoff_b = (int) (SVVCutOff*nmodes_b);
2079  int cutoff_c = (int) (SVVCutOff*nmodes_c);
2080  //To avoid the fac[j] from blowing up
2081  NekDouble epsilon = 1;
2082 
2083  int nmodes = min(min(nmodes_a,nmodes_b),nmodes_c);
2084  NekDouble cutoff = min(min(cutoff_a,cutoff_b),cutoff_c);
2085 
2086  for(i = 0; i < nmodes_a; ++i)//P
2087  {
2088  for(j = 0; j < nmodes_b; ++j) //Q
2089  {
2090  int maxij = max(i,j);
2091  for(k = 0; k < nmodes_c-maxij; ++k) //R
2092  {
2093  if(j + k >= cutoff || i + k >= cutoff)
2094  {
2095  orthocoeffs[cnt] *=
2096  (SvvDiffCoeff*exp(-(i+k-nmodes)*(i+k-nmodes)/
2097  ((NekDouble)((i+k-cutoff+epsilon)*
2098  (i+k-cutoff+epsilon))))*
2099  exp(-(j-nmodes)*(j-nmodes)/
2100  ((NekDouble)((j-cutoff+epsilon)*
2101  (j-cutoff+epsilon)))));
2102  }
2103  else
2104  {
2105  orthocoeffs[cnt] *= 0.0;
2106  }
2107  cnt++;
2108  }
2109  }
2110  }
2111  }
2112 
2113  // backward transform to physical space
2114  OrthoExp.BwdTrans(orthocoeffs,array);
2115  }
2116 
2117  void StdPyrExp::v_ReduceOrderCoeffs(
2118  int numMin,
2119  const Array<OneD, const NekDouble> &inarray,
2120  Array<OneD, NekDouble> &outarray)
2121  {
2122  int nquad0 = m_base[0]->GetNumPoints();
2123  int nquad1 = m_base[1]->GetNumPoints();
2124  int nquad2 = m_base[2]->GetNumPoints();
2125  int nqtot = nquad0*nquad1*nquad2;
2126  int nmodes0 = m_base[0]->GetNumModes();
2127  int nmodes1 = m_base[1]->GetNumModes();
2128  int nmodes2 = m_base[2]->GetNumModes();
2129  int numMax = nmodes0;
2130 
2131  Array<OneD, NekDouble> coeff (m_ncoeffs);
2132  Array<OneD, NekDouble> coeff_tmp1(m_ncoeffs, 0.0);
2133  Array<OneD, NekDouble> phys_tmp (nqtot, 0.0);
2134  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
2135 
2136 
2137  const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
2138  const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
2139  const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
2140 
2141  LibUtilities::BasisKey bortho0(
2142  LibUtilities::eOrtho_A, nmodes0, Pkey0);
2143  LibUtilities::BasisKey bortho1(
2144  LibUtilities::eOrtho_A, nmodes1, Pkey1);
2145  LibUtilities::BasisKey bortho2(
2146  LibUtilities::eOrthoPyr_C, nmodes2, Pkey2);
2147 
2148  int cnt = 0;
2149  int u = 0;
2150  int i = 0;
2151  StdRegions::StdPyrExpSharedPtr OrthoPyrExp;
2152 
2153  OrthoPyrExp = MemoryManager<StdRegions::StdPyrExp>
2154  ::AllocateSharedPtr(bortho0, bortho1, bortho2);
2155 
2156  BwdTrans(inarray,phys_tmp);
2157  OrthoPyrExp->FwdTrans(phys_tmp, coeff);
2158 
2159  // filtering
2160  for (u = 0; u < numMin; ++u)
2161  {
2162  for (i = 0; i < numMin; ++i)
2163  {
2164 
2165  int maxui = max(u,i);
2166  Vmath::Vcopy(numMin - maxui, tmp = coeff + cnt, 1,
2167  tmp2 = coeff_tmp1 + cnt, 1);
2168  cnt += nmodes2 - maxui;
2169  }
2170 
2171  for (i = numMin; i < nmodes1; ++i)
2172  {
2173  int maxui = max(u,i);
2174  cnt += numMax - maxui;
2175  }
2176  }
2177 
2178  OrthoPyrExp->BwdTrans(coeff_tmp1, phys_tmp);
2179  StdPyrExp::FwdTrans(phys_tmp, outarray);
2180  }
2181 
2182  }
2183 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:274
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::LibUtilities::BasisKey::GetBasisType
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:144
Nektar::LibUtilities::BasisKey::GetNumModes
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:86
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:60
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:171
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::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:63
Nektar::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:124
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:158
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:617
Nektar::StdRegions::StdExpansion::LocCoordToLocCollapsed
void LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Convert local cartesian coordinate xi into local collapsed coordinates eta.
Definition: StdExpansion.h:982
Nektar::StdRegions::StdExpansion::CreateGeneralMatrix
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
Definition: StdExpansion.cpp:258
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:376
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:433
Nektar::StdRegions::StdExpansion::GetTraceNcoeffs
int GetTraceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th trace.
Definition: StdExpansion.h:265
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:208
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:1745
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:221
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:171
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1162
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:121
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:130
Nektar::StdRegions::StdPyrExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdPyrExp.cpp:1144
Nektar::StdRegions::StdPyrExp::v_MultiplyByStdQuadratureMetric
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdPyrExp.cpp:1924
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:752
Nektar::StdRegions::StdPyrExp::v_GetTraceNumModes
virtual void v_GetTraceNumModes(const int fid, int &numModes0, int &numModes1, Orientation faceOrient=eDir1FwdDir1_Dir2FwdDir2)
Definition: StdPyrExp.cpp:787
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::StdPyrExp::v_GetTraceNumPoints
virtual int v_GetTraceNumPoints(const int i) const
Definition: StdPyrExp.cpp:995
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:584
Nektar::StdRegions::StdPyrExp::v_PhysEvaluateBasis
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) final
Definition: StdPyrExp.cpp:827
Nektar::StdRegions::StdPyrExp::v_GetTraceBasisKey
virtual const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int k) const
Definition: StdPyrExp.cpp:1034
Nektar::StdRegions::StdPyrExp::v_LocCoordToLocCollapsed
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Definition: StdPyrExp.cpp:722
Nektar::StdRegions::StdPyrExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Definition: StdPyrExp.cpp:1860
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:570
Nektar::StdRegions::StdPyrExp::v_GetEdgeNcoeffs
virtual int v_GetEdgeNcoeffs(const int i) const
Definition: StdPyrExp.cpp:1016
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::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:445
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::StdRegions::StdPyrExp::v_CalcNumberOfCoefficients
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdPyrExp.cpp:1072
Nektar::StdRegions::StdPyrExp::v_GetNverts
virtual int v_GetNverts() const
Definition: StdPyrExp.cpp:891
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_GetTraceToElementMap
virtual void v_GetTraceToElementMap(const int fid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation faceOrient, int P, int Q)
Definition: StdPyrExp.cpp:1245
Nektar::StdRegions::StdPyrExp::v_GetTraceNcoeffs
virtual int v_GetTraceNcoeffs(const int i) const
Definition: StdPyrExp.cpp:953
Nektar::StdRegions::StdPyrExp::v_GetBoundaryMap
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdPyrExp.cpp:1184
Nektar::StdRegions::StdPyrExp::v_NumDGBndryCoeffs
virtual int v_NumDGBndryCoeffs() const
Definition: StdPyrExp.cpp:931
Nektar::StdRegions::StdPyrExp::v_DetShapeType
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdPyrExp.cpp:906
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:1895
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:480
Nektar::StdRegions::StdPyrExp::v_NumBndryCoeffs
virtual int v_NumBndryCoeffs() const
Definition: StdPyrExp.cpp:911
Nektar::StdRegions::StdPyrExp::v_FillMode
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdPyrExp.cpp:780
Nektar::StdRegions::StdPyrExp::v_GetTraceIntNcoeffs
virtual int v_GetTraceIntNcoeffs(const int i) const
Definition: StdPyrExp.cpp:973
Nektar::StdRegions::StdPyrExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdPyrExp.cpp:1865
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(),...
Definition: StdPyrExp.cpp:429
Nektar::StdRegions::StdPyrExp::v_ReduceOrderCoeffs
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdPyrExp.cpp:2117
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:393
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::StdRegions::StdPyrExp::v_GetNedges
virtual int v_GetNedges() const
Definition: StdPyrExp.cpp:896
Nektar::StdRegions::StdPyrExp::v_GetTraceInteriorToElementMap
virtual void v_GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
Definition: StdPyrExp.cpp:1676
Nektar::StdRegions::StdPyrExp::v_GetVertexMap
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdPyrExp.cpp:1085
Nektar::StdRegions::StdPyrExp::v_GetEdgeInteriorToElementMap
virtual void v_GetEdgeInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
Definition: StdPyrExp.cpp:1577
Nektar::StdRegions::StdPyrExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdPyrExp.cpp:1982
Nektar::StdRegions::StdPyrExp::v_GetNtraces
virtual int v_GetNtraces() const
Definition: StdPyrExp.cpp:901
Nektar::StdRegions::StdPyrExp::v_LocCollapsedToLocCoord
virtual void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi)
Definition: StdPyrExp.cpp:743
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:265
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:182
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:197
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 op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition: Blas.hpp:394
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:167
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::LibUtilities::StdPyrData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:240
Nektar::LibUtilities::StdPyrData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:267
Nektar::LibUtilities::NullBasisKey
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
Nektar::LibUtilities::eGaussRadauMAlpha2Beta0
@ eGaussRadauMAlpha2Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:59
Nektar::LibUtilities::eOrtho_A
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:45
Nektar::LibUtilities::eGLL_Lagrange
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:54
Nektar::LibUtilities::eModifiedPyr_C
@ eModifiedPyr_C
Principle Modified Functions .
Definition: BasisType.h:52
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
Nektar::LibUtilities::eOrthoPyr_C
@ eOrthoPyr_C
Principle Orthogonal Functions .
Definition: BasisType.h:51
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:51
Nektar::StdRegions::EvaluateTriFaceBasisKey
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Definition: StdExpansion3D.cpp:499
Nektar::StdRegions::kSVVDGFiltermodesmin
const int kSVVDGFiltermodesmin
Definition: StdRegions.hpp:388
Nektar::StdRegions::StdPyrExpSharedPtr
std::shared_ptr< StdPyrExp > StdPyrExpSharedPtr
Definition: StdPyrExp.h:278
Nektar::StdRegions::kSVVDGFiltermodesmax
const int kSVVDGFiltermodesmax
Definition: StdRegions.hpp:389
Nektar::StdRegions::EvaluateQuadFaceBasisKey
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Definition: StdExpansion3D.cpp:441
Nektar::StdRegions::kSVVDGFilter
const NekDouble kSVVDGFilter[9][11]
Definition: StdRegions.hpp:391
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:820
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
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:192
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:322
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
main.P
P
Definition: main.py:133
Generated on Wed Nov 24 2021 10:33:29 for Nektar++ by doxygen 1.9.1