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