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