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