Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 P,
1004  int Q)
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, k, p, q, r, nFaceCoeffs, idx = 0;
1016  int nummodesA, nummodesB;
1017 
1018  switch (fid)
1019  {
1020  case 0:
1021  nummodesA = m_base[0]->GetNumModes();
1022  nummodesB = m_base[1]->GetNumModes();
1023  break;
1024  case 1:
1025  case 3:
1026  nummodesA = m_base[0]->GetNumModes();
1027  nummodesB = m_base[2]->GetNumModes();
1028  break;
1029  case 2:
1030  case 4:
1031  nummodesA = m_base[1]->GetNumModes();
1032  nummodesB = m_base[2]->GetNumModes();
1033  break;
1034  }
1035 
1036  bool CheckForZeroedModes = false;
1037 
1038  if (P == -1)
1039  {
1040  P = nummodesA;
1041  Q = nummodesB;
1042  nFaceCoeffs = GetFaceNcoeffs(fid);
1043  }
1044  else if (fid == 1 || fid == 3)
1045  {
1046  nFaceCoeffs = P*(2*Q-P+1)/2;
1047  CheckForZeroedModes = true;
1048  }
1049  else
1050  {
1051  nFaceCoeffs = P*Q;
1052  CheckForZeroedModes = true;
1053  }
1054 
1055  // Allocate the map array and sign array; set sign array to ones (+)
1056  if (maparray.num_elements() != nFaceCoeffs)
1057  {
1058  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1059  }
1060 
1061  if (signarray.num_elements() != nFaceCoeffs)
1062  {
1063  signarray = Array<OneD, int>(nFaceCoeffs,1);
1064  }
1065  else
1066  {
1067  fill(signarray.get(), signarray.get() + nFaceCoeffs, 1);
1068  }
1069 
1070  // Set up an array indexing for quads, since the ordering may need
1071  // to be transposed.
1072  Array<OneD, int> arrayindx(nFaceCoeffs,-1);
1073 
1074  if (fid != 1 && fid != 3)
1075  {
1076  for (i = 0; i < Q; i++)
1077  {
1078  for (j = 0; j < P; j++)
1079  {
1080  if (faceOrient < 9)
1081  {
1082  arrayindx[i*P+j] = i*P+j;
1083  }
1084  else
1085  {
1086  arrayindx[i*P+j] = j*Q+i;
1087  }
1088  }
1089  }
1090  }
1091 
1092  // Set up ordering inside each 2D face. Also for triangular faces,
1093  // populate signarray.
1094  switch (fid)
1095  {
1096  case 0: // Bottom quad
1097  for (q = 0; q < Q; ++q)
1098  {
1099  for (p = 0; p < P; ++p)
1100  {
1101  maparray[arrayindx[q*P+p]] = GetMode(p,q,0);
1102  }
1103  }
1104  break;
1105 
1106  case 1: // Left triangle
1107  for (p = 0; p < P; ++p)
1108  {
1109  for (r = 0; r < Q-p; ++r)
1110  {
1111  if ((int)faceOrient == 7 && p > 1)
1112  {
1113  signarray[idx] = p % 2 ? -1 : 1;
1114  }
1115  maparray[idx++] = GetMode(p,0,r);
1116  }
1117  }
1118  break;
1119 
1120  case 2: // Slanted quad
1121  for (q = 0; q < P; ++q)
1122  {
1123  maparray[arrayindx[q]] = GetMode(1,q,0);
1124  }
1125  for (q = 0; q < P; ++q)
1126  {
1127  maparray[arrayindx[P+q]] = GetMode(0,q,1);
1128  }
1129  for (r = 1; r < Q-1; ++r)
1130  {
1131  for (q = 0; q < P; ++q)
1132  {
1133  maparray[arrayindx[(r+1)*P+q]] = GetMode(1,q,r);
1134  }
1135  }
1136  break;
1137 
1138  case 3: // Right triangle
1139  for (p = 0; p < P; ++p)
1140  {
1141  for (r = 0; r < Q-p; ++r)
1142  {
1143  if ((int)faceOrient == 7 && p > 1)
1144  {
1145  signarray[idx] = p % 2 ? -1 : 1;
1146  }
1147  maparray[idx++] = GetMode(p, 1, r);
1148  }
1149  }
1150  break;
1151 
1152  case 4: // Rear quad
1153  for (r = 0; r < Q; ++r)
1154  {
1155  for (q = 0; q < P; ++q)
1156  {
1157  maparray[arrayindx[r*P+q]] = GetMode(0, q, r);
1158  }
1159  }
1160  break;
1161 
1162  default:
1163  ASSERTL0(false, "Face to element map unavailable.");
1164  }
1165 
1166  if (fid == 1 || fid == 3)
1167  {
1168  if(CheckForZeroedModes)
1169  {
1170  // zero signmap and set maparray to zero if elemental
1171  // modes are not as large as face modesl
1172  idx = 0;
1173  for (j = 0; j < nummodesA; ++j)
1174  {
1175  idx += nummodesB-j;
1176  for (k = nummodesB-j; k < Q-j; ++k)
1177  {
1178  signarray[idx] = 0.0;
1179  maparray[idx++] = maparray[0];
1180  }
1181  }
1182 
1183  for (j = nummodesA; j < P; ++j)
1184  {
1185  for (k = 0; k < Q-j; ++k)
1186  {
1187  signarray[idx] = 0.0;
1188  maparray[idx++] = maparray[0];
1189  }
1190  }
1191  }
1192 
1193 
1194  // Triangles only have one possible orientation (base
1195  // direction reversed); swap edge modes.
1196  if ((int)faceOrient == 7)
1197  {
1198  swap(maparray[0], maparray[P]);
1199  for (i = 1; i < P-1; ++i)
1200  {
1201  swap(maparray[i+1], maparray[P+i]);
1202  }
1203  }
1204  }
1205  else
1206  {
1207 
1208  if(CheckForZeroedModes)
1209  {
1210  // zero signmap and set maparray to zero if elemental
1211  // modes are not as large as face modesl
1212  for (j = 0; j < nummodesA; ++j)
1213  {
1214  for (k = nummodesB; k < Q; ++k)
1215  {
1216  signarray[arrayindx[j+k*P]] = 0.0;
1217  maparray[arrayindx[j+k*P]] = maparray[0];
1218  }
1219  }
1220 
1221  for (j = nummodesA; j < P; ++j)
1222  {
1223  for (k = 0; k < Q; ++k)
1224  {
1225  signarray[arrayindx[j+k*P]] = 0.0;
1226  maparray[arrayindx[j+k*P]] = maparray[0];
1227  }
1228  }
1229  }
1230 
1231  // The code below is exactly the same as that taken from
1232  // StdHexExp and reverses the 'b' and 'a' directions as
1233  // appropriate (1st and 2nd if statements respectively) in
1234  // quadrilateral faces.
1235  if (faceOrient == 6 || faceOrient == 8 ||
1236  faceOrient == 11 || faceOrient == 12)
1237  {
1238  if (faceOrient < 9)
1239  {
1240  for (i = 3; i < Q; i += 2)
1241  {
1242  for (j = 0; j < P; j++)
1243  {
1244  signarray[arrayindx[i*P+j]] *= -1;
1245  }
1246  }
1247 
1248  for (i = 0; i < P; i++)
1249  {
1250  swap(maparray [i], maparray [i+P]);
1251  swap(signarray[i], signarray[i+P]);
1252  }
1253  }
1254  else
1255  {
1256  for (i = 0; i < Q; i++)
1257  {
1258  for (j = 3; j < P; j += 2)
1259  {
1260  signarray[arrayindx[i*P+j]] *= -1;
1261  }
1262  }
1263 
1264  for (i = 0; i < Q; i++)
1265  {
1266  swap (maparray [i], maparray [i+Q]);
1267  swap (signarray[i], signarray[i+Q]);
1268  }
1269  }
1270  }
1271 
1272  if (faceOrient == 7 || faceOrient == 8 ||
1273  faceOrient == 10 || faceOrient == 12)
1274  {
1275  if (faceOrient < 9)
1276  {
1277  for (i = 0; i < Q; i++)
1278  {
1279  for (j = 3; j < P; j += 2)
1280  {
1281  signarray[arrayindx[i*P+j]] *= -1;
1282  }
1283  }
1284 
1285  for(i = 0; i < Q; i++)
1286  {
1287  swap(maparray [i*P], maparray [i*P+1]);
1288  swap(signarray[i*P], signarray[i*P+1]);
1289  }
1290  }
1291  else
1292  {
1293  for (i = 3; i < Q; i += 2)
1294  {
1295  for (j = 0; j < P; j++)
1296  {
1297  signarray[arrayindx[i*P+j]] *= -1;
1298  }
1299  }
1300 
1301  for (i = 0; i < P; i++)
1302  {
1303  swap(maparray [i*Q], maparray [i*Q+1]);
1304  swap(signarray[i*Q], signarray[i*Q+1]);
1305  }
1306  }
1307  }
1308  }
1309  }
1310 
1311  int StdPrismExp::v_GetVertexMap(const int vId, bool useCoeffPacking)
1312  {
1313  ASSERTL0(GetEdgeBasisType(vId) == LibUtilities::eModified_A ||
1314  GetEdgeBasisType(vId) == LibUtilities::eModified_A ||
1315  GetEdgeBasisType(vId) == LibUtilities::eModified_B,
1316  "Mapping not defined for this type of basis");
1317 
1318  int l = 0;
1319 
1320  if(useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1321  {
1322  switch (vId)
1323  {
1324  case 0:
1325  l = GetMode(0,0,0);
1326  break;
1327  case 1:
1328  l = GetMode(0,0,1);
1329  break;
1330  case 2:
1331  l = GetMode(0,1,0);
1332  break;
1333  case 3:
1334  l = GetMode(0,1,1);
1335  break;
1336  case 4:
1337  l = GetMode(1,0,0);
1338  break;
1339  case 5:
1340  l = GetMode(1,1,0);
1341  break;
1342  default:
1343  ASSERTL0(false, "local vertex id must be between 0 and 5");
1344  }
1345  }
1346  else
1347  {
1348  switch (vId)
1349  {
1350  case 0:
1351  l = GetMode(0,0,0);
1352  break;
1353  case 1:
1354  l = GetMode(1,0,0);
1355  break;
1356  case 2:
1357  l = GetMode(1,1,0);
1358  break;
1359  case 3:
1360  l = GetMode(0,1,0);
1361  break;
1362  case 4:
1363  l = GetMode(0,0,1);
1364  break;
1365  case 5:
1366  l = GetMode(0,1,1);
1367  break;
1368  default:
1369  ASSERTL0(false, "local vertex id must be between 0 and 5");
1370  }
1371  }
1372 
1373  return l;
1374  }
1375 
1376  void StdPrismExp::v_GetEdgeInteriorMap(
1377  const int eid,
1378  const Orientation edgeOrient,
1379  Array<OneD, unsigned int> &maparray,
1380  Array<OneD, int> &signarray)
1381  {
1382  int i;
1383  bool signChange;
1384  const int P = m_base[0]->GetNumModes() - 1;
1385  const int Q = m_base[1]->GetNumModes() - 1;
1386  const int R = m_base[2]->GetNumModes() - 1;
1387  const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid) - 2;
1388 
1389  if (maparray.num_elements() != nEdgeIntCoeffs)
1390  {
1391  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1392  }
1393 
1394  if(signarray.num_elements() != nEdgeIntCoeffs)
1395  {
1396  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1397  }
1398  else
1399  {
1400  fill(signarray.get(), signarray.get()+nEdgeIntCoeffs, 1);
1401  }
1402 
1403  // If edge is oriented backwards, change sign of modes which have
1404  // degree 2n+1, n >= 1.
1405  signChange = edgeOrient == eBackwards;
1406 
1407  switch (eid)
1408  {
1409  case 0:
1410  for (i = 2; i <= P; ++i)
1411  {
1412  maparray[i-2] = GetMode(i,0,0);
1413  }
1414  break;
1415 
1416  case 1:
1417  for (i = 2; i <= Q; ++i)
1418  {
1419  maparray[i-2] = GetMode(1,i,0);
1420  }
1421  break;
1422 
1423  case 2:
1424  // Base quad; reverse direction.
1425  //signChange = !signChange;
1426  for (i = 2; i <= P; ++i)
1427  {
1428  maparray[i-2] = GetMode(i,1,0);
1429  }
1430  break;
1431 
1432  case 3:
1433  // Base quad; reverse direction.
1434  //signChange = !signChange;
1435  for (i = 2; i <= Q; ++i)
1436  {
1437  maparray[i-2] = GetMode(0,i,0);
1438  }
1439  break;
1440 
1441  case 4:
1442  for (i = 2; i <= R; ++i)
1443  {
1444  maparray[i-2] = GetMode(0,0,i);
1445  }
1446  break;
1447 
1448  case 5:
1449  for (i = 1; i <= R-1; ++i)
1450  {
1451  maparray[i-1] = GetMode(1,0,i);
1452  }
1453  break;
1454 
1455  case 6:
1456  for (i = 1; i <= R-1; ++i)
1457  {
1458  maparray[i-1] = GetMode(1,1,i);
1459  }
1460  break;
1461 
1462  case 7:
1463  for (i = 2; i <= R; ++i)
1464  {
1465  maparray[i-2] = GetMode(0,1,i);
1466  }
1467  break;
1468 
1469  case 8:
1470  for (i = 2; i <= Q; ++i)
1471  {
1472  maparray[i-2] = GetMode(0,i,1);
1473  }
1474  break;
1475 
1476  default:
1477  ASSERTL0(false, "Edge not defined.");
1478  break;
1479  }
1480 
1481  if (signChange)
1482  {
1483  for (i = 1; i < nEdgeIntCoeffs; i += 2)
1484  {
1485  signarray[i] = -1;
1486  }
1487  }
1488  }
1489 
1490  void StdPrismExp::v_GetFaceInteriorMap(
1491  const int fid,
1492  const Orientation faceOrient,
1493  Array<OneD, unsigned int> &maparray,
1494  Array<OneD, int> &signarray)
1495  {
1496  const int P = m_base[0]->GetNumModes() - 1;
1497  const int Q = m_base[1]->GetNumModes() - 1;
1498  const int R = m_base[2]->GetNumModes() - 1;
1499  const int nFaceIntCoeffs = v_GetFaceIntNcoeffs(fid);
1500  int p, q, r, idx = 0;
1501  int nummodesA = 0;
1502  int nummodesB = 0;
1503  int i = 0;
1504  int j = 0;
1505 
1506  if (maparray.num_elements() != nFaceIntCoeffs)
1507  {
1508  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1509  }
1510 
1511  if (signarray.num_elements() != nFaceIntCoeffs)
1512  {
1513  signarray = Array<OneD, int>(nFaceIntCoeffs, 1);
1514  }
1515  else
1516  {
1517  fill(signarray.get(), signarray.get()+nFaceIntCoeffs, 1);
1518  }
1519 
1520  // Set up an array indexing for quad faces, since the ordering may
1521  // need to be transposed depending on orientation.
1522  Array<OneD, int> arrayindx(nFaceIntCoeffs);
1523  if (fid != 1 && fid != 3)
1524  {
1525  if (fid == 0) // Base quad
1526  {
1527  nummodesA = P-1;
1528  nummodesB = Q-1;
1529  }
1530  else // front and back quad
1531  {
1532  nummodesA = Q-1;
1533  nummodesB = R-1;
1534  }
1535 
1536  for (i = 0; i < nummodesB; i++)
1537  {
1538  for (j = 0; j < nummodesA; j++)
1539  {
1540  if (faceOrient < 9)
1541  {
1542  arrayindx[i*nummodesA+j] = i*nummodesA+j;
1543  }
1544  else
1545  {
1546  arrayindx[i*nummodesA+j] = j*nummodesB+i;
1547  }
1548  }
1549  }
1550  }
1551 
1552  switch (fid)
1553  {
1554  case 0: // Bottom quad
1555  for (q = 2; q <= Q; ++q)
1556  {
1557  for (p = 2; p <= P; ++p)
1558  {
1559  maparray[arrayindx[(q-2)*nummodesA+(p-2)]] = GetMode(p,q,0);
1560  }
1561  }
1562  break;
1563 
1564  case 1: // Left triangle
1565  for (p = 2; p <= P; ++p)
1566  {
1567  for (r = 1; r <= R-p; ++r)
1568  {
1569  if ((int)faceOrient == 7)
1570  {
1571  signarray[idx] = p % 2 ? -1 : 1;
1572  }
1573  maparray[idx++] = GetMode(p,0,r);
1574  }
1575  }
1576  break;
1577 
1578  case 2: // Slanted quad
1579  for (r = 1; r <= R-1; ++r)
1580  {
1581  for (q = 2; q <= Q; ++q)
1582  {
1583  maparray[arrayindx[(r-1)*nummodesA+(q-2)]] = GetMode(1, q, r);
1584  }
1585  }
1586  break;
1587 
1588  case 3: // Right triangle
1589  for (p = 2; p <= P; ++p)
1590  {
1591  for (r = 1; r <= R-p; ++r)
1592  {
1593  if ((int)faceOrient == 7)
1594  {
1595  signarray[idx] = p % 2 ? -1 : 1;
1596  }
1597  maparray[idx++] = GetMode(p, 1, r);
1598  }
1599  }
1600  break;
1601 
1602  case 4: // Back quad
1603  for (r = 2; r <= R; ++r)
1604  {
1605  for (q = 2; q <= Q; ++q)
1606  {
1607  maparray[arrayindx[(r-2)*nummodesA+(q-2)]] = GetMode(0, q, r);
1608  }
1609  }
1610  break;
1611 
1612  default:
1613  ASSERTL0(false, "Face interior map not available.");
1614  }
1615 
1616  // Triangular faces are processed in the above switch loop; for
1617  // remaining quad faces, set up orientation if necessary.
1618  if (fid == 1 || fid == 3)
1619  return;
1620 
1621  if (faceOrient == 6 || faceOrient == 8 ||
1622  faceOrient == 11 || faceOrient == 12)
1623  {
1624  if (faceOrient < 9)
1625  {
1626  for (i = 1; i < nummodesB; i += 2)
1627  {
1628  for (j = 0; j < nummodesA; j++)
1629  {
1630  signarray[arrayindx[i*nummodesA+j]] *= -1;
1631  }
1632  }
1633  }
1634  else
1635  {
1636  for (i = 0; i < nummodesB; i++)
1637  {
1638  for (j = 1; j < nummodesA; j += 2)
1639  {
1640  signarray[arrayindx[i*nummodesA+j]] *= -1;
1641  }
1642  }
1643  }
1644  }
1645 
1646  if (faceOrient == 7 || faceOrient == 8 ||
1647  faceOrient == 10 || faceOrient == 12)
1648  {
1649  if (faceOrient < 9)
1650  {
1651  for (i = 0; i < nummodesB; i++)
1652  {
1653  for (j = 1; j < nummodesA; j += 2)
1654  {
1655  signarray[arrayindx[i*nummodesA+j]] *= -1;
1656  }
1657  }
1658  }
1659  else
1660  {
1661  for (i = 1; i < nummodesB; i += 2)
1662  {
1663  for (j = 0; j < nummodesA; j++)
1664  {
1665  signarray[arrayindx[i*nummodesA+j]] *= -1;
1666  }
1667  }
1668  }
1669  }
1670  }
1671 
1672  void StdPrismExp::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
1673  {
1674  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1675  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1676  "BasisType is not a boundary interior form");
1677  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1678  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1679  "BasisType is not a boundary interior form");
1680  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_B ||
1681  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1682  "BasisType is not a boundary interior form");
1683 
1684  int P = m_base[0]->GetNumModes() - 1, p;
1685  int Q = m_base[1]->GetNumModes() - 1, q;
1686  int R = m_base[2]->GetNumModes() - 1, r;
1687 
1688  int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1689 
1690  if(outarray.num_elements()!=nIntCoeffs)
1691  {
1692  outarray = Array<OneD, unsigned int>(nIntCoeffs);
1693  }
1694 
1695  int idx = 0;
1696 
1697  // Loop over all interior modes.
1698  for (p = 2; p <= P; ++p)
1699  {
1700  for (q = 2; q <= Q; ++q)
1701  {
1702  for (r = 1; r <= R-p; ++r)
1703  {
1704  outarray[idx++] = GetMode(p,q,r);
1705  }
1706  }
1707  }
1708  }
1709 
1710  void StdPrismExp::v_GetBoundaryMap(Array<OneD, unsigned int> &maparray)
1711  {
1712  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1713  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1714  "BasisType is not a boundary interior form");
1715  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1716  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1717  "BasisType is not a boundary interior form");
1718  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_B ||
1719  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1720  "BasisType is not a boundary interior form");
1721 
1722  int P = m_base[0]->GetNumModes() - 1, p;
1723  int Q = m_base[1]->GetNumModes() - 1, q;
1724  int R = m_base[2]->GetNumModes() - 1, r;
1725  int idx = 0;
1726 
1727  int nBnd = NumBndryCoeffs();
1728 
1729  if (maparray.num_elements() != nBnd)
1730  {
1731  maparray = Array<OneD, unsigned int>(nBnd);
1732  }
1733 
1734  // Loop over all boundary modes (in ascending order).
1735  for (p = 0; p <= P; ++p)
1736  {
1737  // First two q-r planes are entirely boundary modes.
1738  if (p <= 1)
1739  {
1740  for (q = 0; q <= Q; ++q)
1741  {
1742  for (r = 0; r <= R-p; ++r)
1743  {
1744  maparray[idx++] = GetMode(p,q,r);
1745  }
1746  }
1747  }
1748  else
1749  {
1750  // Remaining q-r planes contain boundary modes on the two
1751  // left-hand sides and bottom edge.
1752  for (q = 0; q <= Q; ++q)
1753  {
1754  if (q <= 1)
1755  {
1756  for (r = 0; r <= R-p; ++r)
1757  {
1758  maparray[idx++] = GetMode(p,q,r);
1759  }
1760  }
1761  else
1762  {
1763  maparray[idx++] = GetMode(p,q,0);
1764  }
1765  }
1766  }
1767  }
1768  }
1769 
1770 
1771 
1772  //---------------------------------------
1773  // Wrapper functions
1774  //---------------------------------------
1775 
1776  DNekMatSharedPtr StdPrismExp::v_GenMatrix(const StdMatrixKey &mkey)
1777  {
1778 
1779  MatrixType mtype = mkey.GetMatrixType();
1780 
1781  DNekMatSharedPtr Mat;
1782 
1783  switch(mtype)
1784  {
1785  case ePhysInterpToEquiSpaced:
1786  {
1787  int nq0 = m_base[0]->GetNumPoints();
1788  int nq1 = m_base[1]->GetNumPoints();
1789  int nq2 = m_base[2]->GetNumPoints();
1790  int nq;
1791 
1792  // take definition from key
1793  if(mkey.ConstFactorExists(eFactorConst))
1794  {
1795  nq = (int) mkey.GetConstFactor(eFactorConst);
1796  }
1797  else
1798  {
1799  nq = max(nq0,max(nq1,nq2));
1800  }
1801 
1802  int neq = LibUtilities::StdPrismData::
1803  getNumberOfCoefficients (nq,nq,nq);
1804  Array<OneD, Array<OneD, NekDouble> > coords (neq);
1805  Array<OneD, NekDouble> coll (3);
1806  Array<OneD, DNekMatSharedPtr> I (3);
1807  Array<OneD, NekDouble> tmp (nq0);
1808 
1809  Mat = MemoryManager<DNekMat>::
1810  AllocateSharedPtr(neq,nq0*nq1*nq2);
1811  int cnt = 0;
1812  for(int i = 0; i < nq; ++i)
1813  {
1814  for(int j = 0; j < nq; ++j)
1815  {
1816  for(int k = 0; k < nq-i; ++k,++cnt)
1817  {
1818  coords[cnt] = Array<OneD, NekDouble>(3);
1819  coords[cnt][0] = -1.0 + 2*k/(NekDouble)(nq-1);
1820  coords[cnt][1] = -1.0 + 2*j/(NekDouble)(nq-1);
1821  coords[cnt][2] = -1.0 + 2*i/(NekDouble)(nq-1);
1822  }
1823  }
1824  }
1825 
1826  for(int i = 0; i < neq; ++i)
1827  {
1828  LocCoordToLocCollapsed(coords[i],coll);
1829 
1830  I[0] = m_base[0]->GetI(coll );
1831  I[1] = m_base[1]->GetI(coll+1);
1832  I[2] = m_base[2]->GetI(coll+2);
1833 
1834  // interpolate first coordinate direction
1835  NekDouble fac;
1836  for( int k = 0; k < nq2; ++k)
1837  {
1838  for (int j = 0; j < nq1; ++j)
1839  {
1840 
1841  fac = (I[1]->GetPtr())[j]*(I[2]->GetPtr())[k];
1842  Vmath::Smul(nq0,fac,I[0]->GetPtr(),1,tmp,1);
1843 
1844  Vmath::Vcopy(nq0, &tmp[0], 1,
1845  Mat->GetRawPtr() +
1846  k * nq0 * nq1 * neq +
1847  j * nq0 * neq + i,
1848  neq);
1849  }
1850  }
1851  }
1852  }
1853  break;
1854  default:
1855  {
1856  Mat = StdExpansion::CreateGeneralMatrix(mkey);
1857  }
1858  break;
1859  }
1860 
1861  return Mat;
1862  }
1863 
1864  DNekMatSharedPtr StdPrismExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
1865  {
1866  return v_GenMatrix(mkey);
1867  }
1868 
1869 
1870 
1871  /**
1872  * @brief Compute the local mode number in the expansion for a
1873  * particular tensorial combination.
1874  *
1875  * Modes are numbered with the r index travelling fastest, followed by
1876  * q and then p, and each q-r plane is of size (R+1-p). For example,
1877  * with P=1, Q=2, R=3, the indexing inside each q-r plane (with r
1878  * increasing upwards and q to the right) is:
1879  *
1880  * p = 0: p = 1:
1881  * -----------------------
1882  * 3 7 11
1883  * 2 6 10 14 17 20
1884  * 1 5 9 13 16 19
1885  * 0 4 8 12 15 18
1886  *
1887  * Note that in this element, we must have that \f$ P <= R \f$.
1888  */
1889  int StdPrismExp::GetMode(int p, int q, int r)
1890  {
1891  int Q = m_base[1]->GetNumModes() - 1;
1892  int R = m_base[2]->GetNumModes() - 1;
1893 
1894  return r + // Skip along stacks (r-direction)
1895  q*(R+1-p) + // Skip along columns (q-direction)
1896  (Q+1)*(p*R + 1-(p-2)*(p-1)/2); // Skip along rows (p-direction)
1897  }
1898 
1899  void StdPrismExp::v_MultiplyByStdQuadratureMetric(
1900  const Array<OneD, const NekDouble>& inarray,
1901  Array<OneD, NekDouble>& outarray)
1902  {
1903  int i, j;
1904  int nquad0 = m_base[0]->GetNumPoints();
1905  int nquad1 = m_base[1]->GetNumPoints();
1906  int nquad2 = m_base[2]->GetNumPoints();
1907 
1908  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1909  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1910  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
1911 
1912  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1913 
1914  // Multiply by integration constants in x-direction
1915  for(i = 0; i < nquad1*nquad2; ++i)
1916  {
1917  Vmath::Vmul(nquad0, inarray.get()+i*nquad0, 1,
1918  w0.get(), 1, outarray.get()+i*nquad0,1);
1919  }
1920 
1921  // Multiply by integration constants in y-direction
1922  for(j = 0; j < nquad2; ++j)
1923  {
1924  for(i = 0; i < nquad1; ++i)
1925  {
1926  Blas::Dscal(nquad0,w1[i], &outarray[0]+i*nquad0 +
1927  j*nquad0*nquad1,1);
1928  }
1929  }
1930 
1931  // Multiply by integration constants in z-direction; need to
1932  // incorporate factor (1-eta_3)/2 into weights, but only if using
1933  // GLL quadrature points.
1934  switch(m_base[2]->GetPointsType())
1935  {
1936  // (1,0) Jacobi inner product.
1937  case LibUtilities::eGaussRadauMAlpha1Beta0:
1938  for(i = 0; i < nquad2; ++i)
1939  {
1940  Blas::Dscal(nquad0*nquad1, 0.5*w2[i],
1941  &outarray[0]+i*nquad0*nquad1, 1);
1942  }
1943  break;
1944 
1945  default:
1946  for(i = 0; i < nquad2; ++i)
1947  {
1948  Blas::Dscal(nquad0*nquad1,0.25*(1-z2[i])*w2[i],
1949  &outarray[0]+i*nquad0*nquad1,1);
1950  }
1951  break;
1952  }
1953 
1954  }
1955 
1956  void StdPrismExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
1957  const StdMatrixKey &mkey)
1958  {
1959  // Generate an orthonogal expansion
1960  int qa = m_base[0]->GetNumPoints();
1961  int qb = m_base[1]->GetNumPoints();
1962  int qc = m_base[2]->GetNumPoints();
1963  int nmodes_a = m_base[0]->GetNumModes();
1964  int nmodes_b = m_base[1]->GetNumModes();
1965  int nmodes_c = m_base[2]->GetNumModes();
1966  // Declare orthogonal basis.
1967  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
1968  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
1969  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
1970 
1971  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
1972  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A,nmodes_b,pb);
1973  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_B,nmodes_c,pc);
1974  StdPrismExp OrthoExp(Ba,Bb,Bc);
1975 
1976  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1977  int i,j,k,cnt = 0;
1978 
1979  //SVV filter paramaters (how much added diffusion relative to physical one
1980  // and fraction of modes from which you start applying this added diffusion)
1981  //
1982  NekDouble SvvDiffCoeff = mkey.GetConstFactor(StdRegions::eFactorSVVDiffCoeff);
1983  NekDouble SVVCutOff = mkey.GetConstFactor(StdRegions::eFactorSVVCutoffRatio);
1984 
1985  //Defining the cut of mode
1986  int cutoff_a = (int) (SVVCutOff*nmodes_a);
1987  int cutoff_b = (int) (SVVCutOff*nmodes_b);
1988  int cutoff_c = (int) (SVVCutOff*nmodes_c);
1989  //To avoid the fac[j] from blowing up
1990  NekDouble epsilon = 1;
1991 
1992  // project onto modal space.
1993  OrthoExp.FwdTrans(array,orthocoeffs);
1994  int nmodes = min(min(nmodes_a,nmodes_b),nmodes_c);
1995  NekDouble cutoff = min(min(cutoff_a,cutoff_b),cutoff_c);
1996 
1997  //------"New" Version August 22nd '13--------------------
1998  for(i = 0; i < nmodes_a; ++i)//P
1999  {
2000  for(j = 0; j < nmodes_b; ++j) //Q
2001  {
2002  for(k = 0; k < nmodes_c-i; ++k) //R
2003  {
2004  if(j >= cutoff || i + k >= cutoff)
2005  {
2006  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)))));
2007  }
2008  else
2009  {
2010  orthocoeffs[cnt] *= 0.0;
2011  }
2012  cnt++;
2013  }
2014  }
2015  }
2016 
2017  // backward transform to physical space
2018  OrthoExp.BwdTrans(orthocoeffs,array);
2019  }
2020 
2021 
2022 
2023  void StdPrismExp::v_ReduceOrderCoeffs(
2024  int numMin,
2025  const Array<OneD, const NekDouble> &inarray,
2026  Array<OneD, NekDouble> &outarray)
2027  {
2028  int nquad0 = m_base[0]->GetNumPoints();
2029  int nquad1 = m_base[1]->GetNumPoints();
2030  int nquad2 = m_base[2]->GetNumPoints();
2031  int nqtot = nquad0*nquad1*nquad2;
2032  int nmodes0 = m_base[0]->GetNumModes();
2033  int nmodes1 = m_base[1]->GetNumModes();
2034  int nmodes2 = m_base[2]->GetNumModes();
2035  int numMax = nmodes0;
2036 
2037  Array<OneD, NekDouble> coeff (m_ncoeffs);
2038  Array<OneD, NekDouble> coeff_tmp1(m_ncoeffs, 0.0);
2039  Array<OneD, NekDouble> phys_tmp (nqtot, 0.0);
2040  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
2041 
2042 
2043  const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
2044  const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
2045  const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
2046 
2047  LibUtilities::BasisKey bortho0(
2048  LibUtilities::eOrtho_A, nmodes0, Pkey0);
2049  LibUtilities::BasisKey bortho1(
2050  LibUtilities::eOrtho_A, nmodes1, Pkey1);
2051  LibUtilities::BasisKey bortho2(
2052  LibUtilities::eOrtho_B, nmodes2, Pkey2);
2053 
2054  int cnt = 0;
2055  int u = 0;
2056  int i = 0;
2057  StdRegions::StdPrismExpSharedPtr OrthoPrismExp;
2058 
2059  OrthoPrismExp = MemoryManager<StdRegions::StdPrismExp>
2060  ::AllocateSharedPtr(bortho0, bortho1, bortho2);
2061 
2062  BwdTrans(inarray,phys_tmp);
2063  OrthoPrismExp->FwdTrans(phys_tmp, coeff);
2064 
2065  // filtering
2066  for (u = 0; u < numMin; ++u)
2067  {
2068  for (i = 0; i < numMin; ++i)
2069  {
2070  Vmath::Vcopy(numMin - u, tmp = coeff + cnt, 1,
2071  tmp2 = coeff_tmp1 + cnt, 1);
2072  cnt += numMax - u;
2073  }
2074 
2075  for (i = numMin; i < numMax; ++i)
2076  {
2077  cnt += numMax - u;
2078  }
2079  }
2080 
2081  OrthoPrismExp->BwdTrans(coeff_tmp1, phys_tmp);
2082  StdPrismExp::FwdTrans(phys_tmp, outarray);
2083  }
2084  }//end namespace
2085 }//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:1376
Nektar::StdRegions::StdPrismExp::v_SVVLaplacianFilter
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdPrismExp.cpp:1956
StdPrismExp.h
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:470
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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:942
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:2023
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:406
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:1173
Nektar::StdRegions::StdPrismExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdPrismExp.cpp:1672
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:700
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::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:131
Nektar::StdRegions::StdPrismExp::v_GetFaceToElementMap
virtual void v_GetFaceToElementMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int P=-1, int Q=-1)
Definition: StdPrismExp.cpp:998
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:413
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:1776
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:349
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:1490
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
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:1710
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:1899
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:525
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:1864
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:1889
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:1422
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
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:1311
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:1905
Generated on Wed Mar 23 2016 20:14:23 for Nektar++ by   doxygen 1.8.8