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.size() > 0)? true:false;
118  bool Do_3 = (out_dxi3.size() > 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.size());
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.size() >= 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  void StdPrismExp::v_LocCoordToLocCollapsed(
684  const Array<OneD, const NekDouble>& xi,
685  Array<OneD, NekDouble>& eta)
686  {
687  NekDouble d2 = 1.0 - xi[2];
688  if(fabs(d2) < NekConstants::kNekZeroTol)
689  {
690  if(d2>=0.)
691  {
692  d2 = NekConstants::kNekZeroTol;
693  }
694  else
695  {
696  d2 = -NekConstants::kNekZeroTol;
697  }
698  }
699  eta[2] = xi[2]; // eta_z = xi_z
700  eta[1] = xi[1]; //eta_y = xi_y
701  eta[0] = 2.0*(1.0 + xi[0])/d2 - 1.0;
702  }
703 
704  void StdPrismExp::v_LocCollapsedToLocCoord(
705  const Array<OneD, const NekDouble>& eta,
706  Array<OneD, NekDouble>& xi)
707  {
708  xi[0] = (1.0 + eta[0]) * (1.0 - eta[2]) * 0.5 - 1.0;
709  xi[1] = eta[1];
710  xi[2] = eta[2];
711  }
712 
713  void StdPrismExp::v_GetCoords(Array<OneD, NekDouble>& xi_x,
714  Array<OneD, NekDouble>& xi_y,
715  Array<OneD, NekDouble>& xi_z)
716  {
717  Array<OneD, const NekDouble> etaBar_x = m_base[0]->GetZ();
718  Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
719  Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
720  int Qx = GetNumPoints(0);
721  int Qy = GetNumPoints(1);
722  int Qz = GetNumPoints(2);
723 
724  // Convert collapsed coordinates into cartesian coordinates: eta --> xi
725  for (int k = 0; k < Qz; ++k) {
726  for (int j = 0; j < Qy; ++j) {
727  for (int i = 0; i < Qx; ++i) {
728  int s = i + Qx*(j + Qy*k);
729  xi_x[s] = (1.0 - eta_z[k])*(1.0 + etaBar_x[i]) / 2.0 - 1.0;
730  xi_y[s] = eta_y[j];
731  xi_z[s] = eta_z[k];
732  }
733  }
734  }
735  }
736 
737  void StdPrismExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
738  {
739  Array<OneD, NekDouble> tmp(m_ncoeffs,0.0);
740  tmp[mode] = 1.0;
741  StdPrismExp::v_BwdTrans(tmp, outarray);
742  }
743 
744  NekDouble StdPrismExp::v_PhysEvaluateBasis(
745  const Array<OneD, const NekDouble>& coords,
746  int mode)
747  {
748  Array<OneD, NekDouble> coll(3);
749  LocCoordToLocCollapsed(coords, coll);
750 
751  const int nm1 = m_base[1]->GetNumModes();
752  const int nm2 = m_base[2]->GetNumModes();
753  const int b = 2 * nm2 + 1;
754 
755  const int mode0 = floor(0.5 * (b - sqrt(b * b - 8.0 * mode / nm1)));
756  const int tmp =
757  mode - nm1*(mode0 * (nm2-1) + 1 - (mode0 - 2)*(mode0 - 1) / 2);
758  const int mode1 = tmp / (nm2 - mode0);
759  const int mode2 = tmp % (nm2 - mode0);
760 
761  if (mode0 == 0 && mode2 == 1 &&
762  m_base[0]->GetBasisType() == LibUtilities::eModified_A)
763  {
764  // handle collapsed top edge to remove mode0 terms
765  return
766  StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
767  StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
768  }
769  else
770  {
771  return
772  StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
773  StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
774  StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
775  }
776  }
777 
778  void StdPrismExp::v_GetTraceNumModes(
779  const int fid,
780  int &numModes0,
781  int &numModes1,
782  Orientation faceOrient)
783  {
784  int nummodes [3] = {m_base[0]->GetNumModes(),
785  m_base[1]->GetNumModes(),
786  m_base[2]->GetNumModes()};
787  switch(fid)
788  {
789  // base quad
790  case 0:
791  {
792  numModes0 = nummodes[0];
793  numModes1 = nummodes[1];
794  }
795  break;
796  // front and back quad
797  case 2:
798  case 4:
799  {
800  numModes0 = nummodes[1];
801  numModes1 = nummodes[2];
802  }
803  break;
804  // triangles
805  case 1:
806  case 3:
807  {
808  numModes0 = nummodes[0];
809  numModes1 = nummodes[2];
810  }
811  break;
812  }
813 
814  if ( faceOrient >= 9 )
815  {
816  std::swap(numModes0, numModes1);
817  }
818  }
819 
820  int StdPrismExp::v_GetEdgeNcoeffs(const int i) const
821  {
822  ASSERTL2(i >= 0 && i <= 8, "edge id is out of range");
823 
824  if (i == 0 || i == 2)
825  {
826  return GetBasisNumModes(0);
827  }
828  else if (i == 1 || i == 3 || i == 8)
829  {
830  return GetBasisNumModes(1);
831  }
832  else
833  {
834  return GetBasisNumModes(2);
835  }
836  }
837 
838  //---------------------------------------
839  // Helper functions
840  //---------------------------------------
841 
842  int StdPrismExp::v_GetNverts() const
843  {
844  return 6;
845  }
846 
847  int StdPrismExp::v_GetNedges() const
848  {
849  return 9;
850  }
851 
852  int StdPrismExp::v_GetNtraces() const
853  {
854  return 5;
855  }
856 
857  /**
858  * \brief Return Shape of region, using ShapeType enum list;
859  * i.e. prism.
860  */
861  LibUtilities::ShapeType StdPrismExp::v_DetShapeType() const
862  {
863  return LibUtilities::ePrism;
864  }
865 
866  int StdPrismExp::v_NumBndryCoeffs() const
867  {
868  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
869  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
870  "BasisType is not a boundary interior form");
871  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
872  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
873  "BasisType is not a boundary interior form");
874  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_B ||
875  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
876  "BasisType is not a boundary interior form");
877 
878  int P = m_base[0]->GetNumModes();
879  int Q = m_base[1]->GetNumModes();
880  int R = m_base[2]->GetNumModes();
881 
882  return LibUtilities::StdPrismData::
883  getNumberOfBndCoefficients(P,Q,R);
884  }
885 
886  int StdPrismExp::v_NumDGBndryCoeffs() const
887  {
888  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
889  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
890  "BasisType is not a boundary interior form");
891  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
892  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
893  "BasisType is not a boundary interior form");
894  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_B ||
895  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
896  "BasisType is not a boundary interior form");
897 
898  int P = m_base[0]->GetNumModes()-1;
899  int Q = m_base[1]->GetNumModes()-1;
900  int R = m_base[2]->GetNumModes()-1;
901 
902  return (P+1)*(Q+1) // 1 rect. face on base
903  + 2*(Q+1)*(R+1) // other 2 rect. faces
904  + 2*(R+1) + P*(1 + 2*R - P); // 2 tri. faces
905  }
906 
907 
908  int StdPrismExp::v_GetTraceNcoeffs(const int i) const
909  {
910  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
911  if (i == 0)
912  {
913  return GetBasisNumModes(0)*GetBasisNumModes(1);
914  }
915  else if (i == 1 || i == 3)
916  {
917  int P = GetBasisNumModes(0)-1, Q = GetBasisNumModes(2)-1;
918  return Q+1 + (P*(1 + 2*Q - P))/2;
919  }
920  else
921  {
922  return GetBasisNumModes(1)*GetBasisNumModes(2);
923  }
924  }
925 
926  int StdPrismExp::v_GetTraceIntNcoeffs(const int i) const
927  {
928  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
929 
930  int Pi = GetBasisNumModes(0) - 2;
931  int Qi = GetBasisNumModes(1) - 2;
932  int Ri = GetBasisNumModes(2) - 2;
933 
934  if (i == 0)
935  {
936  return Pi * Qi;
937  }
938  else if (i == 1 || i == 3)
939  {
940  return Pi * (2*Ri - Pi - 1) / 2;
941  }
942  else
943  {
944  return Qi * Ri;
945  }
946  }
947 
948  int StdPrismExp::v_GetTotalTraceIntNcoeffs() const
949  {
950  int Pi = GetBasisNumModes(0) - 2;
951  int Qi = GetBasisNumModes(1) - 2;
952  int Ri = GetBasisNumModes(2) - 2;
953 
954  return Pi * Qi +
955  Pi * (2*Ri - Pi - 1) +
956  2* Qi * Ri;
957  }
958 
959  int StdPrismExp::v_GetTraceNumPoints(const int i) const
960  {
961  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
962 
963  if (i == 0)
964  {
965  return m_base[0]->GetNumPoints()*
966  m_base[1]->GetNumPoints();
967  }
968  else if (i == 1 || i == 3)
969  {
970  return m_base[0]->GetNumPoints()*
971  m_base[2]->GetNumPoints();
972  }
973  else
974  {
975  return m_base[1]->GetNumPoints()*
976  m_base[2]->GetNumPoints();
977  }
978  }
979 
980  LibUtilities::PointsKey StdPrismExp::v_GetTracePointsKey(
981  const int i, const int j) const
982  {
983  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
984  ASSERTL2(j == 0 || j == 1, "face direction is out of range");
985 
986  if (i == 0)
987  {
988  return m_base[j]->GetPointsKey();
989  }
990  else if (i == 1 || i == 3)
991  {
992  return m_base[2*j]->GetPointsKey();
993  }
994  else
995  {
996  return m_base[j+1]->GetPointsKey();
997  }
998  }
999 
1000  const LibUtilities::BasisKey StdPrismExp::v_GetTraceBasisKey(
1001  const int i, const int k) const
1002  {
1003  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
1004  ASSERTL2(k >= 0 && k <= 1, "basis key id is out of range");
1005 
1006  switch(i)
1007  {
1008  case 0:
1009  {
1010  return EvaluateQuadFaceBasisKey(k,
1011  m_base[k]->GetBasisType(),
1012  m_base[k]->GetNumPoints(),
1013  m_base[k]->GetNumModes());
1014  }
1015  case 2:
1016  case 4:
1017  {
1018  return EvaluateQuadFaceBasisKey(k,
1019  m_base[k+1]->GetBasisType(),
1020  m_base[k+1]->GetNumPoints(),
1021  m_base[k+1]->GetNumModes());
1022  }
1023  case 1:
1024  case 3:
1025  {
1026  return EvaluateTriFaceBasisKey(k,
1027  m_base[2*k]->GetBasisType(),
1028  m_base[2*k]->GetNumPoints(),
1029  m_base[2*k]->GetNumModes());
1030 
1031  }
1032  break;
1033  }
1034 
1035  // Should never get here.
1036  return LibUtilities::NullBasisKey;
1037  }
1038 
1039  int StdPrismExp::v_CalcNumberOfCoefficients(const std::vector<unsigned int> &nummodes,
1040  int &modes_offset)
1041  {
1042  int nmodes = LibUtilities::StdPrismData::getNumberOfCoefficients(
1043  nummodes[modes_offset],
1044  nummodes[modes_offset+1],
1045  nummodes[modes_offset+2]);
1046 
1047  modes_offset += 3;
1048  return nmodes;
1049  }
1050 
1051  bool StdPrismExp::v_IsBoundaryInteriorExpansion()
1052  {
1053  return (m_base[0]->GetBasisType() == LibUtilities::eModified_A) &&
1054  (m_base[1]->GetBasisType() == LibUtilities::eModified_A) &&
1055  (m_base[2]->GetBasisType() == LibUtilities::eModified_B);
1056  }
1057 
1058  //---------------------------------------
1059  // Mappings
1060  //---------------------------------------
1061 
1062  int StdPrismExp::v_GetVertexMap(const int vId, bool useCoeffPacking)
1063  {
1064  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1065  GetBasisType(1) == LibUtilities::eModified_A ||
1066  GetBasisType(2) == LibUtilities::eModified_B,
1067  "Mapping not defined for this type of basis");
1068 
1069  int l = 0;
1070 
1071  if(useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1072  {
1073  switch (vId)
1074  {
1075  case 0:
1076  l = GetMode(0,0,0);
1077  break;
1078  case 1:
1079  l = GetMode(0,0,1);
1080  break;
1081  case 2:
1082  l = GetMode(0,1,0);
1083  break;
1084  case 3:
1085  l = GetMode(0,1,1);
1086  break;
1087  case 4:
1088  l = GetMode(1,0,0);
1089  break;
1090  case 5:
1091  l = GetMode(1,1,0);
1092  break;
1093  default:
1094  ASSERTL0(false, "local vertex id must be between 0 and 5");
1095  }
1096  }
1097  else
1098  {
1099  switch (vId)
1100  {
1101  case 0:
1102  l = GetMode(0,0,0);
1103  break;
1104  case 1:
1105  l = GetMode(1,0,0);
1106  break;
1107  case 2:
1108  l = GetMode(1,1,0);
1109  break;
1110  case 3:
1111  l = GetMode(0,1,0);
1112  break;
1113  case 4:
1114  l = GetMode(0,0,1);
1115  break;
1116  case 5:
1117  l = GetMode(0,1,1);
1118  break;
1119  default:
1120  ASSERTL0(false, "local vertex id must be between 0 and 5");
1121  }
1122  }
1123 
1124  return l;
1125  }
1126 
1127  void StdPrismExp::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
1128  {
1129  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1130  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1131  "BasisType is not a boundary interior form");
1132  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1133  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1134  "BasisType is not a boundary interior form");
1135  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_B ||
1136  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1137  "BasisType is not a boundary interior form");
1138 
1139  int P = m_base[0]->GetNumModes() - 1, p;
1140  int Q = m_base[1]->GetNumModes() - 1, q;
1141  int R = m_base[2]->GetNumModes() - 1, r;
1142 
1143  int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1144 
1145  if(outarray.size()!=nIntCoeffs)
1146  {
1147  outarray = Array<OneD, unsigned int>(nIntCoeffs);
1148  }
1149 
1150  int idx = 0;
1151 
1152  // Loop over all interior modes.
1153  for (p = 2; p <= P; ++p)
1154  {
1155  for (q = 2; q <= Q; ++q)
1156  {
1157  for (r = 1; r <= R-p; ++r)
1158  {
1159  outarray[idx++] = GetMode(p,q,r);
1160  }
1161  }
1162  }
1163  }
1164 
1165  void StdPrismExp::v_GetBoundaryMap(Array<OneD, unsigned int> &maparray)
1166  {
1167  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1168  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1169  "BasisType is not a boundary interior form");
1170  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1171  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1172  "BasisType is not a boundary interior form");
1173  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_B ||
1174  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1175  "BasisType is not a boundary interior form");
1176 
1177  int P = m_base[0]->GetNumModes() - 1, p;
1178  int Q = m_base[1]->GetNumModes() - 1, q;
1179  int R = m_base[2]->GetNumModes() - 1, r;
1180  int idx = 0;
1181 
1182  int nBnd = NumBndryCoeffs();
1183 
1184  if (maparray.size() != nBnd)
1185  {
1186  maparray = Array<OneD, unsigned int>(nBnd);
1187  }
1188 
1189  // Loop over all boundary modes (in ascending order).
1190  for (p = 0; p <= P; ++p)
1191  {
1192  // First two q-r planes are entirely boundary modes.
1193  if (p <= 1)
1194  {
1195  for (q = 0; q <= Q; ++q)
1196  {
1197  for (r = 0; r <= R-p; ++r)
1198  {
1199  maparray[idx++] = GetMode(p,q,r);
1200  }
1201  }
1202  }
1203  else
1204  {
1205  // Remaining q-r planes contain boundary modes on the two
1206  // left-hand sides and bottom edge.
1207  for (q = 0; q <= Q; ++q)
1208  {
1209  if (q <= 1)
1210  {
1211  for (r = 0; r <= R-p; ++r)
1212  {
1213  maparray[idx++] = GetMode(p,q,r);
1214  }
1215  }
1216  else
1217  {
1218  maparray[idx++] = GetMode(p,q,0);
1219  }
1220  }
1221  }
1222  }
1223  }
1224 
1225  void StdPrismExp::v_GetTraceToElementMap(
1226  const int fid,
1227  Array<OneD, unsigned int> &maparray,
1228  Array<OneD, int> &signarray,
1229  const Orientation faceOrient,
1230  int P,
1231  int Q)
1232  {
1233  ASSERTL1(GetBasisType(0) == GetBasisType(1),
1234  "Method only implemented if BasisType is identical"
1235  "in x and y directions");
1236  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A &&
1237  GetBasisType(2) == LibUtilities::eModified_B,
1238  "Method only implemented for Modified_A BasisType"
1239  "(x and y direction) and Modified_B BasisType (z "
1240  "direction)");
1241 
1242  int i, j, k, p, q, r, nFaceCoeffs, idx = 0;
1243  int nummodesA=0, nummodesB=0;
1244 
1245  switch (fid)
1246  {
1247  case 0:
1248  nummodesA = m_base[0]->GetNumModes();
1249  nummodesB = m_base[1]->GetNumModes();
1250  break;
1251  case 1:
1252  case 3:
1253  nummodesA = m_base[0]->GetNumModes();
1254  nummodesB = m_base[2]->GetNumModes();
1255  break;
1256  case 2:
1257  case 4:
1258  nummodesA = m_base[1]->GetNumModes();
1259  nummodesB = m_base[2]->GetNumModes();
1260  break;
1261  default:
1262  ASSERTL0(false,"fid must be between 0 and 4");
1263  }
1264 
1265  bool CheckForZeroedModes = false;
1266 
1267  if (P == -1)
1268  {
1269  P = nummodesA;
1270  Q = nummodesB;
1271  nFaceCoeffs = GetTraceNcoeffs(fid);
1272  }
1273  else if (fid == 1 || fid == 3)
1274  {
1275  nFaceCoeffs = P*(2*Q-P+1)/2;
1276  CheckForZeroedModes = true;
1277  }
1278  else
1279  {
1280  nFaceCoeffs = P*Q;
1281  CheckForZeroedModes = true;
1282  }
1283 
1284  // Allocate the map array and sign array; set sign array to ones (+)
1285  if (maparray.size() != nFaceCoeffs)
1286  {
1287  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1288  }
1289 
1290  if (signarray.size() != nFaceCoeffs)
1291  {
1292  signarray = Array<OneD, int>(nFaceCoeffs,1);
1293  }
1294  else
1295  {
1296  fill(signarray.get(), signarray.get() + nFaceCoeffs, 1);
1297  }
1298 
1299  // Set up an array indexing for quads, since the ordering may need
1300  // to be transposed.
1301  Array<OneD, int> arrayindx(nFaceCoeffs,-1);
1302 
1303  if (fid != 1 && fid != 3)
1304  {
1305  for (i = 0; i < Q; i++)
1306  {
1307  for (j = 0; j < P; j++)
1308  {
1309  if (faceOrient < 9)
1310  {
1311  arrayindx[i*P+j] = i*P+j;
1312  }
1313  else
1314  {
1315  arrayindx[i*P+j] = j*Q+i;
1316  }
1317  }
1318  }
1319  }
1320 
1321  // Set up ordering inside each 2D face. Also for triangular faces,
1322  // populate signarray.
1323  switch (fid)
1324  {
1325  case 0: // Bottom quad
1326  for (q = 0; q < Q; ++q)
1327  {
1328  for (p = 0; p < P; ++p)
1329  {
1330  maparray[arrayindx[q*P+p]] = GetMode(p,q,0);
1331  }
1332  }
1333  break;
1334 
1335  case 1: // Left triangle
1336  for (p = 0; p < P; ++p)
1337  {
1338  for (r = 0; r < Q-p; ++r)
1339  {
1340  if ((int)faceOrient == 7 && p > 1)
1341  {
1342  signarray[idx] = p % 2 ? -1 : 1;
1343  }
1344  maparray[idx++] = GetMode(p,0,r);
1345  }
1346  }
1347  break;
1348 
1349  case 2: // Slanted quad
1350  for (q = 0; q < P; ++q)
1351  {
1352  maparray[arrayindx[q]] = GetMode(1,q,0);
1353  }
1354  for (q = 0; q < P; ++q)
1355  {
1356  maparray[arrayindx[P+q]] = GetMode(0,q,1);
1357  }
1358  for (r = 1; r < Q-1; ++r)
1359  {
1360  for (q = 0; q < P; ++q)
1361  {
1362  maparray[arrayindx[(r+1)*P+q]] = GetMode(1,q,r);
1363  }
1364  }
1365  break;
1366 
1367  case 3: // Right triangle
1368  for (p = 0; p < P; ++p)
1369  {
1370  for (r = 0; r < Q-p; ++r)
1371  {
1372  if ((int)faceOrient == 7 && p > 1)
1373  {
1374  signarray[idx] = p % 2 ? -1 : 1;
1375  }
1376  maparray[idx++] = GetMode(p, 1, r);
1377  }
1378  }
1379  break;
1380 
1381  case 4: // Rear quad
1382  for (r = 0; r < Q; ++r)
1383  {
1384  for (q = 0; q < P; ++q)
1385  {
1386  maparray[arrayindx[r*P+q]] = GetMode(0, q, r);
1387  }
1388  }
1389  break;
1390 
1391  default:
1392  ASSERTL0(false, "Face to element map unavailable.");
1393  }
1394 
1395  if (fid == 1 || fid == 3)
1396  {
1397  if(CheckForZeroedModes)
1398  {
1399  // zero signmap and set maparray to zero if elemental
1400  // modes are not as large as face modesl
1401  idx = 0;
1402  for (j = 0; j < nummodesA; ++j)
1403  {
1404  idx += nummodesB-j;
1405  for (k = nummodesB-j; k < Q-j; ++k)
1406  {
1407  signarray[idx] = 0.0;
1408  maparray[idx++] = maparray[0];
1409  }
1410  }
1411 
1412  for (j = nummodesA; j < P; ++j)
1413  {
1414  for (k = 0; k < Q-j; ++k)
1415  {
1416  signarray[idx] = 0.0;
1417  maparray[idx++] = maparray[0];
1418  }
1419  }
1420  }
1421 
1422 
1423  // Triangles only have one possible orientation (base
1424  // direction reversed); swap edge modes.
1425  if ((int)faceOrient == 7)
1426  {
1427  swap(maparray[0], maparray[Q]);
1428  for (i = 1; i < Q-1; ++i)
1429  {
1430  swap(maparray[i+1], maparray[Q+i]);
1431  }
1432  }
1433  }
1434  else
1435  {
1436 
1437  if(CheckForZeroedModes)
1438  {
1439  // zero signmap and set maparray to zero if elemental
1440  // modes are not as large as face modesl
1441  for (j = 0; j < nummodesA; ++j)
1442  {
1443  for (k = nummodesB; k < Q; ++k)
1444  {
1445  signarray[arrayindx[j+k*P]] = 0.0;
1446  maparray[arrayindx[j+k*P]] = maparray[0];
1447  }
1448  }
1449 
1450  for (j = nummodesA; j < P; ++j)
1451  {
1452  for (k = 0; k < Q; ++k)
1453  {
1454  signarray[arrayindx[j+k*P]] = 0.0;
1455  maparray[arrayindx[j+k*P]] = maparray[0];
1456  }
1457  }
1458  }
1459 
1460  // The code below is exactly the same as that taken from
1461  // StdHexExp and reverses the 'b' and 'a' directions as
1462  // appropriate (1st and 2nd if statements respectively) in
1463  // quadrilateral faces.
1464  if (faceOrient == 6 || faceOrient == 8 ||
1465  faceOrient == 11 || faceOrient == 12)
1466  {
1467  if (faceOrient < 9)
1468  {
1469  for (i = 3; i < Q; i += 2)
1470  {
1471  for (j = 0; j < P; j++)
1472  {
1473  signarray[arrayindx[i*P+j]] *= -1;
1474  }
1475  }
1476 
1477  for (i = 0; i < P; i++)
1478  {
1479  swap(maparray [i], maparray [i+P]);
1480  swap(signarray[i], signarray[i+P]);
1481  }
1482  }
1483  else
1484  {
1485  for (i = 0; i < Q; i++)
1486  {
1487  for (j = 3; j < P; j += 2)
1488  {
1489  signarray[arrayindx[i*P+j]] *= -1;
1490  }
1491  }
1492 
1493  for (i = 0; i < Q; i++)
1494  {
1495  swap (maparray [i], maparray [i+Q]);
1496  swap (signarray[i], signarray[i+Q]);
1497  }
1498  }
1499  }
1500 
1501  if (faceOrient == 7 || faceOrient == 8 ||
1502  faceOrient == 10 || faceOrient == 12)
1503  {
1504  if (faceOrient < 9)
1505  {
1506  for (i = 0; i < Q; i++)
1507  {
1508  for (j = 3; j < P; j += 2)
1509  {
1510  signarray[arrayindx[i*P+j]] *= -1;
1511  }
1512  }
1513 
1514  for(i = 0; i < Q; i++)
1515  {
1516  swap(maparray [i*P], maparray [i*P+1]);
1517  swap(signarray[i*P], signarray[i*P+1]);
1518  }
1519  }
1520  else
1521  {
1522  for (i = 3; i < Q; i += 2)
1523  {
1524  for (j = 0; j < P; j++)
1525  {
1526  signarray[arrayindx[i*P+j]] *= -1;
1527  }
1528  }
1529 
1530  for (i = 0; i < P; i++)
1531  {
1532  swap(maparray [i*Q], maparray [i*Q+1]);
1533  swap(signarray[i*Q], signarray[i*Q+1]);
1534  }
1535  }
1536  }
1537  }
1538  }
1539 
1540 
1541  void StdPrismExp::v_GetEdgeInteriorToElementMap(
1542  const int eid,
1543  Array<OneD, unsigned int> &maparray,
1544  Array<OneD, int> &signarray,
1545  const Orientation edgeOrient)
1546  {
1547  int i;
1548  bool signChange;
1549  const int P = m_base[0]->GetNumModes() - 1;
1550  const int Q = m_base[1]->GetNumModes() - 1;
1551  const int R = m_base[2]->GetNumModes() - 1;
1552  const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid) - 2;
1553 
1554  if (maparray.size() != nEdgeIntCoeffs)
1555  {
1556  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1557  }
1558 
1559  if(signarray.size() != nEdgeIntCoeffs)
1560  {
1561  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1562  }
1563  else
1564  {
1565  fill(signarray.get(), signarray.get()+nEdgeIntCoeffs, 1);
1566  }
1567 
1568  // If edge is oriented backwards, change sign of modes which have
1569  // degree 2n+1, n >= 1.
1570  signChange = edgeOrient == eBackwards;
1571 
1572  switch (eid)
1573  {
1574  case 0:
1575  for (i = 2; i <= P; ++i)
1576  {
1577  maparray[i-2] = GetMode(i,0,0);
1578  }
1579  break;
1580 
1581  case 1:
1582  for (i = 2; i <= Q; ++i)
1583  {
1584  maparray[i-2] = GetMode(1,i,0);
1585  }
1586  break;
1587 
1588  case 2:
1589  // Base quad; reverse direction.
1590  //signChange = !signChange;
1591  for (i = 2; i <= P; ++i)
1592  {
1593  maparray[i-2] = GetMode(i,1,0);
1594  }
1595  break;
1596 
1597  case 3:
1598  // Base quad; reverse direction.
1599  //signChange = !signChange;
1600  for (i = 2; i <= Q; ++i)
1601  {
1602  maparray[i-2] = GetMode(0,i,0);
1603  }
1604  break;
1605 
1606  case 4:
1607  for (i = 2; i <= R; ++i)
1608  {
1609  maparray[i-2] = GetMode(0,0,i);
1610  }
1611  break;
1612 
1613  case 5:
1614  for (i = 1; i <= R-1; ++i)
1615  {
1616  maparray[i-1] = GetMode(1,0,i);
1617  }
1618  break;
1619 
1620  case 6:
1621  for (i = 1; i <= R-1; ++i)
1622  {
1623  maparray[i-1] = GetMode(1,1,i);
1624  }
1625  break;
1626 
1627  case 7:
1628  for (i = 2; i <= R; ++i)
1629  {
1630  maparray[i-2] = GetMode(0,1,i);
1631  }
1632  break;
1633 
1634  case 8:
1635  for (i = 2; i <= Q; ++i)
1636  {
1637  maparray[i-2] = GetMode(0,i,1);
1638  }
1639  break;
1640 
1641  default:
1642  ASSERTL0(false, "Edge not defined.");
1643  break;
1644  }
1645 
1646  if (signChange)
1647  {
1648  for (i = 1; i < nEdgeIntCoeffs; i += 2)
1649  {
1650  signarray[i] = -1;
1651  }
1652  }
1653  }
1654 
1655  void StdPrismExp::v_GetTraceInteriorToElementMap(
1656  const int fid,
1657  Array<OneD, unsigned int> &maparray,
1658  Array<OneD, int> &signarray,
1659  const Orientation faceOrient)
1660  {
1661  const int P = m_base[0]->GetNumModes() - 1;
1662  const int Q = m_base[1]->GetNumModes() - 1;
1663  const int R = m_base[2]->GetNumModes() - 1;
1664  const int nFaceIntCoeffs = v_GetTraceIntNcoeffs(fid);
1665  int p, q, r, idx = 0;
1666  int nummodesA = 0;
1667  int nummodesB = 0;
1668  int i = 0;
1669  int j = 0;
1670 
1671  if (maparray.size() != nFaceIntCoeffs)
1672  {
1673  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1674  }
1675 
1676  if (signarray.size() != nFaceIntCoeffs)
1677  {
1678  signarray = Array<OneD, int>(nFaceIntCoeffs, 1);
1679  }
1680  else
1681  {
1682  fill(signarray.get(), signarray.get()+nFaceIntCoeffs, 1);
1683  }
1684 
1685  // Set up an array indexing for quad faces, since the ordering may
1686  // need to be transposed depending on orientation.
1687  Array<OneD, int> arrayindx(nFaceIntCoeffs);
1688  if (fid != 1 && fid != 3)
1689  {
1690  if (fid == 0) // Base quad
1691  {
1692  nummodesA = P-1;
1693  nummodesB = Q-1;
1694  }
1695  else // front and back quad
1696  {
1697  nummodesA = Q-1;
1698  nummodesB = R-1;
1699  }
1700 
1701  for (i = 0; i < nummodesB; i++)
1702  {
1703  for (j = 0; j < nummodesA; j++)
1704  {
1705  if (faceOrient < 9)
1706  {
1707  arrayindx[i*nummodesA+j] = i*nummodesA+j;
1708  }
1709  else
1710  {
1711  arrayindx[i*nummodesA+j] = j*nummodesB+i;
1712  }
1713  }
1714  }
1715  }
1716 
1717  switch (fid)
1718  {
1719  case 0: // Bottom quad
1720  for (q = 2; q <= Q; ++q)
1721  {
1722  for (p = 2; p <= P; ++p)
1723  {
1724  maparray[arrayindx[(q-2)*nummodesA+(p-2)]] = GetMode(p,q,0);
1725  }
1726  }
1727  break;
1728 
1729  case 1: // Left triangle
1730  for (p = 2; p <= P; ++p)
1731  {
1732  for (r = 1; r <= R-p; ++r)
1733  {
1734  if ((int)faceOrient == 7)
1735  {
1736  signarray[idx] = p % 2 ? -1 : 1;
1737  }
1738  maparray[idx++] = GetMode(p,0,r);
1739  }
1740  }
1741  break;
1742 
1743  case 2: // Slanted quad
1744  for (r = 1; r <= R-1; ++r)
1745  {
1746  for (q = 2; q <= Q; ++q)
1747  {
1748  maparray[arrayindx[(r-1)*nummodesA+(q-2)]] = GetMode(1, q, r);
1749  }
1750  }
1751  break;
1752 
1753  case 3: // Right triangle
1754  for (p = 2; p <= P; ++p)
1755  {
1756  for (r = 1; r <= R-p; ++r)
1757  {
1758  if ((int)faceOrient == 7)
1759  {
1760  signarray[idx] = p % 2 ? -1 : 1;
1761  }
1762  maparray[idx++] = GetMode(p, 1, r);
1763  }
1764  }
1765  break;
1766 
1767  case 4: // Back quad
1768  for (r = 2; r <= R; ++r)
1769  {
1770  for (q = 2; q <= Q; ++q)
1771  {
1772  maparray[arrayindx[(r-2)*nummodesA+(q-2)]] = GetMode(0, q, r);
1773  }
1774  }
1775  break;
1776 
1777  default:
1778  ASSERTL0(false, "Face interior map not available.");
1779  }
1780 
1781  // Triangular faces are processed in the above switch loop; for
1782  // remaining quad faces, set up orientation if necessary.
1783  if (fid == 1 || fid == 3)
1784  return;
1785 
1786  if (faceOrient == 6 || faceOrient == 8 ||
1787  faceOrient == 11 || faceOrient == 12)
1788  {
1789  if (faceOrient < 9)
1790  {
1791  for (i = 1; i < nummodesB; i += 2)
1792  {
1793  for (j = 0; j < nummodesA; j++)
1794  {
1795  signarray[arrayindx[i*nummodesA+j]] *= -1;
1796  }
1797  }
1798  }
1799  else
1800  {
1801  for (i = 0; i < nummodesB; i++)
1802  {
1803  for (j = 1; j < nummodesA; j += 2)
1804  {
1805  signarray[arrayindx[i*nummodesA+j]] *= -1;
1806  }
1807  }
1808  }
1809  }
1810 
1811  if (faceOrient == 7 || faceOrient == 8 ||
1812  faceOrient == 10 || faceOrient == 12)
1813  {
1814  if (faceOrient < 9)
1815  {
1816  for (i = 0; i < nummodesB; i++)
1817  {
1818  for (j = 1; j < nummodesA; j += 2)
1819  {
1820  signarray[arrayindx[i*nummodesA+j]] *= -1;
1821  }
1822  }
1823  }
1824  else
1825  {
1826  for (i = 1; i < nummodesB; i += 2)
1827  {
1828  for (j = 0; j < nummodesA; j++)
1829  {
1830  signarray[arrayindx[i*nummodesA+j]] *= -1;
1831  }
1832  }
1833  }
1834  }
1835  }
1836 
1837  //---------------------------------------
1838  // Wrapper functions
1839  //---------------------------------------
1840 
1841  DNekMatSharedPtr StdPrismExp::v_GenMatrix(const StdMatrixKey &mkey)
1842  {
1843 
1844  MatrixType mtype = mkey.GetMatrixType();
1845 
1846  DNekMatSharedPtr Mat;
1847 
1848  switch(mtype)
1849  {
1850  case ePhysInterpToEquiSpaced:
1851  {
1852  int nq0 = m_base[0]->GetNumPoints();
1853  int nq1 = m_base[1]->GetNumPoints();
1854  int nq2 = m_base[2]->GetNumPoints();
1855  int nq;
1856 
1857  // take definition from key
1858  if(mkey.ConstFactorExists(eFactorConst))
1859  {
1860  nq = (int) mkey.GetConstFactor(eFactorConst);
1861  }
1862  else
1863  {
1864  nq = max(nq0,max(nq1,nq2));
1865  }
1866 
1867  int neq = LibUtilities::StdPrismData::
1868  getNumberOfCoefficients (nq,nq,nq);
1869  Array<OneD, Array<OneD, NekDouble> > coords (neq);
1870  Array<OneD, NekDouble> coll (3);
1871  Array<OneD, DNekMatSharedPtr> I (3);
1872  Array<OneD, NekDouble> tmp (nq0);
1873 
1874  Mat = MemoryManager<DNekMat>::
1875  AllocateSharedPtr(neq,nq0*nq1*nq2);
1876  int cnt = 0;
1877  for(int i = 0; i < nq; ++i)
1878  {
1879  for(int j = 0; j < nq; ++j)
1880  {
1881  for(int k = 0; k < nq-i; ++k,++cnt)
1882  {
1883  coords[cnt] = Array<OneD, NekDouble>(3);
1884  coords[cnt][0] = -1.0 + 2*k/(NekDouble)(nq-1);
1885  coords[cnt][1] = -1.0 + 2*j/(NekDouble)(nq-1);
1886  coords[cnt][2] = -1.0 + 2*i/(NekDouble)(nq-1);
1887  }
1888  }
1889  }
1890 
1891  for(int i = 0; i < neq; ++i)
1892  {
1893  LocCoordToLocCollapsed(coords[i],coll);
1894 
1895  I[0] = m_base[0]->GetI(coll );
1896  I[1] = m_base[1]->GetI(coll+1);
1897  I[2] = m_base[2]->GetI(coll+2);
1898 
1899  // interpolate first coordinate direction
1900  NekDouble fac;
1901  for( int k = 0; k < nq2; ++k)
1902  {
1903  for (int j = 0; j < nq1; ++j)
1904  {
1905 
1906  fac = (I[1]->GetPtr())[j]*(I[2]->GetPtr())[k];
1907  Vmath::Smul(nq0,fac,I[0]->GetPtr(),1,tmp,1);
1908 
1909  Vmath::Vcopy(nq0, &tmp[0], 1,
1910  Mat->GetRawPtr() +
1911  k * nq0 * nq1 * neq +
1912  j * nq0 * neq + i,
1913  neq);
1914  }
1915  }
1916  }
1917  }
1918  break;
1919  default:
1920  {
1921  Mat = StdExpansion::CreateGeneralMatrix(mkey);
1922  }
1923  break;
1924  }
1925 
1926  return Mat;
1927  }
1928 
1929  DNekMatSharedPtr StdPrismExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
1930  {
1931  return v_GenMatrix(mkey);
1932  }
1933 
1934 
1935 
1936  /**
1937  * @brief Compute the local mode number in the expansion for a
1938  * particular tensorial combination.
1939  *
1940  * Modes are numbered with the r index travelling fastest, followed by
1941  * q and then p, and each q-r plane is of size (R+1-p). For example,
1942  * with P=1, Q=2, R=3, the indexing inside each q-r plane (with r
1943  * increasing upwards and q to the right) is:
1944  *
1945  * p = 0: p = 1:
1946  * -----------------------
1947  * 3 7 11
1948  * 2 6 10 14 17 20
1949  * 1 5 9 13 16 19
1950  * 0 4 8 12 15 18
1951  *
1952  * Note that in this element, we must have that \f$ P <= R \f$.
1953  */
1954  int StdPrismExp::GetMode(int p, int q, int r)
1955  {
1956  int Q = m_base[1]->GetNumModes() - 1;
1957  int R = m_base[2]->GetNumModes() - 1;
1958 
1959  return r + // Skip along stacks (r-direction)
1960  q*(R+1-p) + // Skip along columns (q-direction)
1961  (Q+1)*(p*R + 1-(p-2)*(p-1)/2); // Skip along rows (p-direction)
1962  }
1963 
1964  void StdPrismExp::v_MultiplyByStdQuadratureMetric(
1965  const Array<OneD, const NekDouble>& inarray,
1966  Array<OneD, NekDouble>& outarray)
1967  {
1968  int i, j;
1969  int nquad0 = m_base[0]->GetNumPoints();
1970  int nquad1 = m_base[1]->GetNumPoints();
1971  int nquad2 = m_base[2]->GetNumPoints();
1972 
1973  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1974  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1975  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
1976 
1977  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1978 
1979  // Multiply by integration constants in x-direction
1980  for(i = 0; i < nquad1*nquad2; ++i)
1981  {
1982  Vmath::Vmul(nquad0, inarray.get()+i*nquad0, 1,
1983  w0.get(), 1, outarray.get()+i*nquad0,1);
1984  }
1985 
1986  // Multiply by integration constants in y-direction
1987  for(j = 0; j < nquad2; ++j)
1988  {
1989  for(i = 0; i < nquad1; ++i)
1990  {
1991  Blas::Dscal(nquad0,w1[i], &outarray[0]+i*nquad0 +
1992  j*nquad0*nquad1,1);
1993  }
1994  }
1995 
1996  // Multiply by integration constants in z-direction; need to
1997  // incorporate factor (1-eta_3)/2 into weights, but only if using
1998  // GLL quadrature points.
1999  switch(m_base[2]->GetPointsType())
2000  {
2001  // (1,0) Jacobi inner product.
2002  case LibUtilities::eGaussRadauMAlpha1Beta0:
2003  for(i = 0; i < nquad2; ++i)
2004  {
2005  Blas::Dscal(nquad0*nquad1, 0.5*w2[i],
2006  &outarray[0]+i*nquad0*nquad1, 1);
2007  }
2008  break;
2009 
2010  default:
2011  for(i = 0; i < nquad2; ++i)
2012  {
2013  Blas::Dscal(nquad0*nquad1,0.5*(1-z2[i])*w2[i],
2014  &outarray[0]+i*nquad0*nquad1,1);
2015  }
2016  break;
2017  }
2018 
2019  }
2020 
2021  void StdPrismExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
2022  const StdMatrixKey &mkey)
2023  {
2024  // Generate an orthonogal expansion
2025  int qa = m_base[0]->GetNumPoints();
2026  int qb = m_base[1]->GetNumPoints();
2027  int qc = m_base[2]->GetNumPoints();
2028  int nmodes_a = m_base[0]->GetNumModes();
2029  int nmodes_b = m_base[1]->GetNumModes();
2030  int nmodes_c = m_base[2]->GetNumModes();
2031  // Declare orthogonal basis.
2032  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
2033  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
2034  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
2035 
2036  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
2037  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A,nmodes_b,pb);
2038  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_B,nmodes_c,pc);
2039  StdPrismExp OrthoExp(Ba,Bb,Bc);
2040 
2041  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2042  int i,j,k,cnt = 0;
2043 
2044  // project onto modal space.
2045  OrthoExp.FwdTrans(array,orthocoeffs);
2046 
2047  if(mkey.ConstFactorExists(eFactorSVVPowerKerDiffCoeff))
2048  {
2049  // Rodrigo's power kernel
2050  NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
2051  NekDouble SvvDiffCoeff =
2052  mkey.GetConstFactor(eFactorSVVPowerKerDiffCoeff)*
2053  mkey.GetConstFactor(eFactorSVVDiffCoeff);
2054 
2055  for(int i = 0; i < nmodes_a; ++i)
2056  {
2057  for(int j = 0; j < nmodes_b; ++j)
2058  {
2059  NekDouble fac1 = std::max(
2060  pow((1.0*i)/(nmodes_a-1),cutoff*nmodes_a),
2061  pow((1.0*j)/(nmodes_b-1),cutoff*nmodes_b));
2062 
2063  for(int k = 0; k < nmodes_c-i; ++k)
2064  {
2065  NekDouble fac = std::max(fac1,
2066  pow((1.0*k)/(nmodes_c-1),cutoff*nmodes_c));
2067 
2068  orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2069  cnt++;
2070  }
2071  }
2072  }
2073  }
2074  else if(mkey.ConstFactorExists(eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
2075  {
2076  NekDouble SvvDiffCoeff =
2077  mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff)*
2078  mkey.GetConstFactor(eFactorSVVDiffCoeff);
2079 
2080  int max_abc = max(nmodes_a-kSVVDGFiltermodesmin,
2081  nmodes_b-kSVVDGFiltermodesmin);
2082  max_abc = max(max_abc, nmodes_c-kSVVDGFiltermodesmin);
2083  // clamp max_abc
2084  max_abc = max(max_abc,0);
2085  max_abc = min(max_abc,kSVVDGFiltermodesmax-kSVVDGFiltermodesmin);
2086 
2087  for(int i = 0; i < nmodes_a; ++i)
2088  {
2089  for(int j = 0; j < nmodes_b; ++j)
2090  {
2091  int maxij = max(i,j);
2092 
2093  for(int k = 0; k < nmodes_c-i; ++k)
2094  {
2095  int maxijk = max(maxij,k);
2096  maxijk = min(maxijk,kSVVDGFiltermodesmax-1);
2097 
2098  orthocoeffs[cnt] *= SvvDiffCoeff *
2099  kSVVDGFilter[max_abc][maxijk];
2100  cnt++;
2101  }
2102  }
2103  }
2104  }
2105  else
2106  {
2107  // SVV filter paramaters (how much added diffusion relative
2108  // to physical one and fraction of modes from which you
2109  // start applying this added diffusion)
2110  //
2111  NekDouble SvvDiffCoeff = mkey.GetConstFactor(StdRegions::eFactorSVVDiffCoeff);
2112  NekDouble SVVCutOff = mkey.GetConstFactor(StdRegions::eFactorSVVCutoffRatio);
2113 
2114  //Defining the cut of mode
2115  int cutoff_a = (int) (SVVCutOff*nmodes_a);
2116  int cutoff_b = (int) (SVVCutOff*nmodes_b);
2117  int cutoff_c = (int) (SVVCutOff*nmodes_c);
2118  //To avoid the fac[j] from blowing up
2119  NekDouble epsilon = 1;
2120 
2121  int nmodes = min(min(nmodes_a,nmodes_b),nmodes_c);
2122  NekDouble cutoff = min(min(cutoff_a,cutoff_b),cutoff_c);
2123 
2124  //------"New" Version August 22nd '13--------------------
2125  for(i = 0; i < nmodes_a; ++i)//P
2126  {
2127  for(j = 0; j < nmodes_b; ++j) //Q
2128  {
2129  for(k = 0; k < nmodes_c-i; ++k) //R
2130  {
2131  if(j >= cutoff || i + k >= cutoff)
2132  {
2133  orthocoeffs[cnt] *=
2134  (SvvDiffCoeff*exp(-(i+k-nmodes)*(i+k-nmodes)/
2135  ((NekDouble)((i+k-cutoff+epsilon)*
2136  (i+k-cutoff+epsilon))))*
2137  exp(-(j-nmodes)*(j-nmodes)/
2138  ((NekDouble)((j-cutoff+epsilon)*
2139  (j-cutoff+epsilon)))));
2140  }
2141  else
2142  {
2143  orthocoeffs[cnt] *= 0.0;
2144  }
2145  cnt++;
2146  }
2147  }
2148  }
2149  }
2150 
2151  // backward transform to physical space
2152  OrthoExp.BwdTrans(orthocoeffs,array);
2153  }
2154 
2155 
2156 
2157  void StdPrismExp::v_ReduceOrderCoeffs(
2158  int numMin,
2159  const Array<OneD, const NekDouble> &inarray,
2160  Array<OneD, NekDouble> &outarray)
2161  {
2162  int nquad0 = m_base[0]->GetNumPoints();
2163  int nquad1 = m_base[1]->GetNumPoints();
2164  int nquad2 = m_base[2]->GetNumPoints();
2165  int nqtot = nquad0*nquad1*nquad2;
2166  int nmodes0 = m_base[0]->GetNumModes();
2167  int nmodes1 = m_base[1]->GetNumModes();
2168  int nmodes2 = m_base[2]->GetNumModes();
2169  int numMax = nmodes0;
2170 
2171  Array<OneD, NekDouble> coeff (m_ncoeffs);
2172  Array<OneD, NekDouble> coeff_tmp1(m_ncoeffs, 0.0);
2173  Array<OneD, NekDouble> phys_tmp (nqtot, 0.0);
2174  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
2175 
2176 
2177  const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
2178  const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
2179  const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
2180 
2181  LibUtilities::BasisKey bortho0(
2182  LibUtilities::eOrtho_A, nmodes0, Pkey0);
2183  LibUtilities::BasisKey bortho1(
2184  LibUtilities::eOrtho_A, nmodes1, Pkey1);
2185  LibUtilities::BasisKey bortho2(
2186  LibUtilities::eOrtho_B, nmodes2, Pkey2);
2187 
2188  int cnt = 0;
2189  int u = 0;
2190  int i = 0;
2191  StdRegions::StdPrismExpSharedPtr OrthoPrismExp;
2192 
2193  OrthoPrismExp = MemoryManager<StdRegions::StdPrismExp>
2194  ::AllocateSharedPtr(bortho0, bortho1, bortho2);
2195 
2196  BwdTrans(inarray,phys_tmp);
2197  OrthoPrismExp->FwdTrans(phys_tmp, coeff);
2198 
2199  // filtering
2200  for (u = 0; u < numMin; ++u)
2201  {
2202  for (i = 0; i < numMin; ++i)
2203  {
2204  Vmath::Vcopy(numMin - u, tmp = coeff + cnt, 1,
2205  tmp2 = coeff_tmp1 + cnt, 1);
2206  cnt += numMax - u;
2207  }
2208 
2209  for (i = numMin; i < numMax; ++i)
2210  {
2211  cnt += numMax - u;
2212  }
2213  }
2214 
2215  OrthoPrismExp->BwdTrans(coeff_tmp1, phys_tmp);
2216  StdPrismExp::FwdTrans(phys_tmp, outarray);
2217  }
2218  }//end namespace
2219 }//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:171
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: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::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:134
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::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:733
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::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_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_GetTraceBasisKey
virtual const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int k) const
Definition: StdPrismExp.cpp:1000
Nektar::StdRegions::StdPrismExp::v_GetTraceInteriorToElementMap
virtual void v_GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
Definition: StdPrismExp.cpp:1655
Nektar::StdRegions::StdPrismExp::v_ReduceOrderCoeffs
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdPrismExp.cpp:2157
Nektar::StdRegions::StdPrismExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Definition: StdPrismExp.cpp:1841
Nektar::StdRegions::StdPrismExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdPrismExp.cpp:1929
Nektar::StdRegions::StdPrismExp::v_GetTracePointsKey
virtual LibUtilities::PointsKey v_GetTracePointsKey(const int i, const int j) const
Definition: StdPrismExp.cpp:980
Nektar::StdRegions::StdPrismExp::v_GetTraceNumPoints
virtual int v_GetTraceNumPoints(const int i) const
Definition: StdPrismExp.cpp:959
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_GetNedges
virtual int v_GetNedges() const
Definition: StdPrismExp.cpp:847
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_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:713
Nektar::StdRegions::StdPrismExp::v_GetTraceIntNcoeffs
virtual int v_GetTraceIntNcoeffs(const int i) const
Definition: StdPrismExp.cpp:926
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_GetTraceToElementMap
virtual void v_GetTraceToElementMap(const int fid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation faceOrient, int P, int Q)
Definition: StdPrismExp.cpp:1225
Nektar::StdRegions::StdPrismExp::v_GetTotalTraceIntNcoeffs
virtual int v_GetTotalTraceIntNcoeffs() const
Definition: StdPrismExp.cpp:948
Nektar::StdRegions::StdPrismExp::v_GetVertexMap
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdPrismExp.cpp:1062
Nektar::StdRegions::StdPrismExp::v_SVVLaplacianFilter
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdPrismExp.cpp:2021
Nektar::StdRegions::StdPrismExp::v_GetTraceNcoeffs
virtual int v_GetTraceNcoeffs(const int i) const
Definition: StdPrismExp.cpp:908
Nektar::StdRegions::StdPrismExp::v_GetNtraces
virtual int v_GetNtraces() const
Definition: StdPrismExp.cpp:852
Nektar::StdRegions::StdPrismExp::v_PhysEvaluateBasis
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) final
Definition: StdPrismExp.cpp:744
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:1165
Nektar::StdRegions::StdPrismExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdPrismExp.cpp:1127
Nektar::StdRegions::StdPrismExp::v_GetTraceNumModes
virtual void v_GetTraceNumModes(const int fid, int &numModes0, int &numModes1, Orientation faceOrient=eDir1FwdDir1_Dir2FwdDir2)
Definition: StdPrismExp.cpp:778
Nektar::StdRegions::StdPrismExp::v_GetNverts
virtual int v_GetNverts() const
Definition: StdPrismExp.cpp:842
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:1954
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:737
Nektar::StdRegions::StdPrismExp::v_GetEdgeNcoeffs
virtual int v_GetEdgeNcoeffs(const int i) const
Definition: StdPrismExp.cpp:820
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:886
Nektar::StdRegions::StdPrismExp::v_CalcNumberOfCoefficients
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdPrismExp.cpp:1039
Nektar::StdRegions::StdPrismExp::v_MultiplyByStdQuadratureMetric
virtual void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdPrismExp.cpp:1964
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_GetEdgeInteriorToElementMap
virtual void v_GetEdgeInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
Definition: StdPrismExp.cpp:1541
Nektar::StdRegions::StdPrismExp::v_IsBoundaryInteriorExpansion
virtual bool v_IsBoundaryInteriorExpansion()
Definition: StdPrismExp.cpp:1051
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:866
Nektar::StdRegions::StdPrismExp::v_LocCollapsedToLocCoord
virtual void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi)
Definition: StdPrismExp.cpp:704
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:861
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:683
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::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:51
Nektar::StdRegions::StdPrismExpSharedPtr
std::shared_ptr< StdPrismExp > StdPrismExpSharedPtr
Definition: StdPrismExp.h:263
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::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::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:565
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
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267
Generated on Wed Nov 24 2021 10:33:29 for Nektar++ by doxygen 1.9.1