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  {
436  int nquad1 = m_base[1]->GetNumPoints();
437  int nquad2 = m_base[2]->GetNumPoints();
438  int order0 = m_base[0]->GetNumModes();
439  int order1 = m_base[1]->GetNumModes();
440 
441  Array<OneD, NekDouble> tmp(inarray.num_elements());
442  Array<OneD, NekDouble> wsp(order0*nquad2*(nquad1+order1));
443 
444  MultiplyByQuadratureMetric(inarray,tmp);
445 
446  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
447  m_base[1]->GetBdata(),
448  m_base[2]->GetBdata(),
449  tmp,outarray,wsp,
450  true,true,true);
451  }
452 
453  void StdPrismExp::v_IProductWRTBase_SumFacKernel(
454  const Array<OneD, const NekDouble>& base0,
455  const Array<OneD, const NekDouble>& base1,
456  const Array<OneD, const NekDouble>& base2,
457  const Array<OneD, const NekDouble>& inarray,
458  Array<OneD, NekDouble> &outarray,
459  Array<OneD, NekDouble> &wsp,
460  bool doCheckCollDir0,
461  bool doCheckCollDir1,
462  bool doCheckCollDir2)
463  {
464  // Interior prism implementation based on Spen's book page
465  // 119. and 608.
466  const int nquad0 = m_base[0]->GetNumPoints();
467  const int nquad1 = m_base[1]->GetNumPoints();
468  const int nquad2 = m_base[2]->GetNumPoints();
469  const int order0 = m_base[0]->GetNumModes ();
470  const int order1 = m_base[1]->GetNumModes ();
471  const int order2 = m_base[2]->GetNumModes ();
472 
473  int i, mode;
474 
475  ASSERTL1(wsp.num_elements() >= nquad1*nquad2*order0 +
476  nquad2*order0*order1,
477  "Insufficient workspace size");
478 
479  Array<OneD, NekDouble> tmp0 = wsp;
480  Array<OneD, NekDouble> tmp1 = wsp + nquad1*nquad2*order0;
481 
482  // Inner product with respect to the '0' direction
483  Blas::Dgemm('T', 'N', nquad1*nquad2, order0, nquad0,
484  1.0, inarray.get(), nquad0,
485  base0.get(), nquad0,
486  0.0, tmp0.get(), nquad1*nquad2);
487 
488  // Inner product with respect to the '1' direction
489  Blas::Dgemm('T', 'N', nquad2*order0, order1, nquad1,
490  1.0, tmp0.get(), nquad1,
491  base1.get(), nquad1,
492  0.0, tmp1.get(), nquad2*order0);
493 
494  // Inner product with respect to the '2' direction
495  for (mode=i=0; i < order0; ++i)
496  {
497  Blas::Dgemm('T', 'N', order2-i, order1, nquad2,
498  1.0, base2.get() + mode*nquad2, nquad2,
499  tmp1.get() + i*nquad2, nquad2*order0,
500  0.0, outarray.get()+mode*order1, order2-i);
501  mode += order2-i;
502  }
503 
504  // Fix top singular vertices; performs phi_{0,q,1} +=
505  // phi_1(xi_1)*phi_q(xi_2)*phi_{01}*phi_r(xi_2).
506  if (m_base[0]->GetBasisType() == LibUtilities::eModified_A)
507  {
508  for (i = 0; i < order1; ++i)
509  {
510  mode = GetMode(0,i,1);
511  outarray[mode] += Blas::Ddot(
512  nquad2, base2.get()+nquad2, 1,
513  tmp1.get()+i*order0*nquad2+nquad2, 1);
514  }
515  }
516  }
517 
518  /**
519  * \brief Inner product of \a inarray over region with respect to the
520  * object's default expansion basis; output in \a outarray.
521  */
522  void StdPrismExp::v_IProductWRTDerivBase(
523  const int dir,
524  const Array<OneD, const NekDouble>& inarray,
525  Array<OneD, NekDouble>& outarray)
526  {
527  v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
528  }
529 
530  void StdPrismExp::v_IProductWRTDerivBase_MatOp(
531  const int dir,
532  const Array<OneD, const NekDouble>& inarray,
533  Array<OneD, NekDouble>& outarray)
534  {
535  ASSERTL0(dir >= 0 && dir <= 2, "input dir is out of range");
536 
537  int nq = GetTotPoints();
538  MatrixType mtype;
539 
540  switch (dir)
541  {
542  case 0:
543  mtype = eIProductWRTDerivBase0;
544  break;
545  case 1:
546  mtype = eIProductWRTDerivBase1;
547  break;
548  case 2:
549  mtype = eIProductWRTDerivBase2;
550  break;
551  }
552 
553  StdMatrixKey iprodmatkey(mtype,DetShapeType(),*this);
554  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
555 
556  Blas::Dgemv('N',m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
557  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
558  }
559 
560  void StdPrismExp::v_IProductWRTDerivBase_SumFac(
561  const int dir,
562  const Array<OneD, const NekDouble>& inarray,
563  Array<OneD, NekDouble>& outarray)
564  {
565  ASSERTL0(dir >= 0 && dir <= 2, "input dir is out of range");
566 
567  int i;
568  int order0 = m_base[0]->GetNumModes ();
569  int order1 = m_base[1]->GetNumModes ();
570  int nquad0 = m_base[0]->GetNumPoints();
571  int nquad1 = m_base[1]->GetNumPoints();
572  int nquad2 = m_base[2]->GetNumPoints();
573 
574  const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
575  const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
576  Array<OneD, NekDouble> gfac0(nquad0);
577  Array<OneD, NekDouble> gfac2(nquad2);
578  Array<OneD, NekDouble> tmp0 (nquad0*nquad1*nquad2);
579  Array<OneD, NekDouble> wsp (order0*nquad2*(nquad1+order1));
580 
581  // set up geometric factor: (1+z0)/2
582  for (i = 0; i < nquad0; ++i)
583  {
584  gfac0[i] = 0.5*(1+z0[i]);
585  }
586 
587  // Set up geometric factor: 2/(1-z2)
588  for (i = 0; i < nquad2; ++i)
589  {
590  gfac2[i] = 2.0/(1-z2[i]);
591  }
592 
593  // Scale first derivative term by gfac2.
594  if (dir != 1)
595  {
596  for (i = 0; i < nquad2; ++i)
597  {
598  Vmath::Smul(nquad0*nquad1,gfac2[i],
599  &inarray[0]+i*nquad0*nquad1,1,
600  &tmp0 [0]+i*nquad0*nquad1,1);
601  }
602  MultiplyByQuadratureMetric(tmp0,tmp0);
603  }
604 
605  switch (dir)
606  {
607  case 0:
608  {
609  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
610  m_base[1]->GetBdata (),
611  m_base[2]->GetBdata (),
612  tmp0,outarray,wsp,
613  true,true,true);
614  break;
615  }
616 
617  case 1:
618  {
619  MultiplyByQuadratureMetric(inarray,tmp0);
620  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata (),
621  m_base[1]->GetDbdata(),
622  m_base[2]->GetBdata (),
623  tmp0,outarray,wsp,
624  true,true,true);
625  break;
626  }
627 
628  case 2:
629  {
630  Array<OneD, NekDouble> tmp1(m_ncoeffs);
631 
632  // Scale eta_1 derivative with gfac0.
633  for(i = 0; i < nquad1*nquad2; ++i)
634  {
635  Vmath::Vmul(nquad0,&gfac0[0],1,&tmp0[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
636  }
637 
638  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
639  m_base[1]->GetBdata(),
640  m_base[2]->GetBdata(),
641  tmp0,tmp1,wsp,
642  true,true,true);
643 
644  MultiplyByQuadratureMetric(inarray, tmp0);
645  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
646  m_base[1]->GetBdata(),
647  m_base[2]->GetDbdata(),
648  tmp0,outarray,wsp,
649  true,true,true);
650 
651  Vmath::Vadd(m_ncoeffs,&tmp1[0],1,&outarray[0],1,&outarray[0],1);
652  break;
653  }
654  }
655  }
656 
657 
658  //---------------------------------------
659  // Evaluation functions
660  //---------------------------------------
661 
662 
663 
664  void StdPrismExp::v_LocCoordToLocCollapsed(
665  const Array<OneD, const NekDouble>& xi,
666  Array<OneD, NekDouble>& eta)
667  {
668 
669  if( fabs(xi[2]-1.0) < NekConstants::kNekZeroTol)
670  {
671  // Very top point of the prism
672  eta[0] = -1.0;
673  eta[1] = xi[1];
674  eta[2] = 1.0;
675  }
676  else
677  {
678  // Third basis function collapsed to "pr" direction instead of
679  // "qr" direction
680  eta[2] = xi[2]; // eta_z = xi_z
681  eta[1] = xi[1]; //eta_y = xi_y
682  eta[0] = 2.0*(1.0 + xi[0])/(1.0 - xi[2]) - 1.0;
683  }
684  }
685 
686  void StdPrismExp::v_GetCoords(Array<OneD, NekDouble>& xi_x,
687  Array<OneD, NekDouble>& xi_y,
688  Array<OneD, NekDouble>& xi_z)
689  {
690  Array<OneD, const NekDouble> etaBar_x = m_base[0]->GetZ();
691  Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
692  Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
693  int Qx = GetNumPoints(0);
694  int Qy = GetNumPoints(1);
695  int Qz = GetNumPoints(2);
696 
697  // Convert collapsed coordinates into cartesian coordinates: eta --> xi
698  for (int k = 0; k < Qz; ++k) {
699  for (int j = 0; j < Qy; ++j) {
700  for (int i = 0; i < Qx; ++i) {
701  int s = i + Qx*(j + Qy*k);
702  xi_x[s] = (1.0 - eta_z[k])*(1.0 + etaBar_x[i]) / 2.0 - 1.0;
703  xi_y[s] = eta_y[j];
704  xi_z[s] = eta_z[k];
705  }
706  }
707  }
708  }
709 
710  void StdPrismExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
711  {
712  Array<OneD, NekDouble> tmp(m_ncoeffs,0.0);
713  tmp[mode] = 1.0;
714  StdPrismExp::v_BwdTrans(tmp, outarray);
715  }
716 
717 
718  //---------------------------------------
719  // Helper functions
720  //---------------------------------------
721 
722  int StdPrismExp::v_GetNverts() const
723  {
724  return 6;
725  }
726 
727  int StdPrismExp::v_GetNedges() const
728  {
729  return 9;
730  }
731 
732  int StdPrismExp::v_GetNfaces() const
733  {
734  return 5;
735  }
736 
737  /**
738  * \brief Return Shape of region, using ShapeType enum list;
739  * i.e. prism.
740  */
741  LibUtilities::ShapeType StdPrismExp::v_DetShapeType() const
742  {
743  return LibUtilities::ePrism;
744  }
745 
746  int StdPrismExp::v_NumBndryCoeffs() const
747  {
748  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
749  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
750  "BasisType is not a boundary interior form");
751  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
752  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
753  "BasisType is not a boundary interior form");
754  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_B ||
755  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
756  "BasisType is not a boundary interior form");
757 
758  int P = m_base[0]->GetNumModes();
759  int Q = m_base[1]->GetNumModes();
760  int R = m_base[2]->GetNumModes();
761 
762  return LibUtilities::StdPrismData::
763  getNumberOfBndCoefficients(P,Q,R);
764  }
765 
766  int StdPrismExp::v_NumDGBndryCoeffs() const
767  {
768  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
769  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
770  "BasisType is not a boundary interior form");
771  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
772  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
773  "BasisType is not a boundary interior form");
774  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_B ||
775  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
776  "BasisType is not a boundary interior form");
777 
778  int P = m_base[0]->GetNumModes()-1;
779  int Q = m_base[1]->GetNumModes()-1;
780  int R = m_base[2]->GetNumModes()-1;
781 
782  return (P+1)*(Q+1) // 1 rect. face on base
783  + 2*(Q+1)*(R+1) // other 2 rect. faces
784  + 2*(R+1) + P*(1 + 2*R - P); // 2 tri. faces
785  }
786 
787  int StdPrismExp::v_GetEdgeNcoeffs(const int i) const
788  {
789  ASSERTL2(i >= 0 && i <= 8, "edge id is out of range");
790 
791  if (i == 0 || i == 2)
792  {
793  return GetBasisNumModes(0);
794  }
795  else if (i == 1 || i == 3 || i == 8)
796  {
797  return GetBasisNumModes(1);
798  }
799  else
800  {
801  return GetBasisNumModes(2);
802  }
803  }
804 
805  int StdPrismExp::v_GetTotalEdgeIntNcoeffs() const
806  {
807  int P = GetBasisNumModes(0)-2;
808  int Q = GetBasisNumModes(1)-2;
809  int R = GetBasisNumModes(2)-2;
810 
811  return 2*P+3*Q+3*R;
812  }
813 
814  int StdPrismExp::v_GetFaceNcoeffs(const int i) const
815  {
816  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
817  if (i == 0)
818  {
819  return GetBasisNumModes(0)*GetBasisNumModes(1);
820  }
821  else if (i == 1 || i == 3)
822  {
823  int P = GetBasisNumModes(0)-1, Q = GetBasisNumModes(2)-1;
824  return Q+1 + (P*(1 + 2*Q - P))/2;
825  }
826  else
827  {
828  return GetBasisNumModes(1)*GetBasisNumModes(2);
829  }
830  }
831 
832  int StdPrismExp::v_GetFaceIntNcoeffs(const int i) const
833  {
834  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
835 
836  int Pi = GetBasisNumModes(0) - 2;
837  int Qi = GetBasisNumModes(1) - 2;
838  int Ri = GetBasisNumModes(2) - 2;
839 
840  if (i == 0)
841  {
842  return Pi * Qi;
843  }
844  else if (i == 1 || i == 3)
845  {
846  return Pi * (2*Ri - Pi - 1) / 2;
847  }
848  else
849  {
850  return Qi * Ri;
851  }
852  }
853 
854  int StdPrismExp::v_GetTotalFaceIntNcoeffs() const
855  {
856  int Pi = GetBasisNumModes(0) - 2;
857  int Qi = GetBasisNumModes(1) - 2;
858  int Ri = GetBasisNumModes(2) - 2;
859 
860  return Pi * Qi +
861  Pi * (2*Ri - Pi - 1) +
862  2* Qi * Ri;
863  }
864 
865  int StdPrismExp::v_GetFaceNumPoints(const int i) const
866  {
867  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
868 
869  if (i == 0)
870  {
871  return m_base[0]->GetNumPoints()*
872  m_base[1]->GetNumPoints();
873  }
874  else if (i == 1 || i == 3)
875  {
876  return m_base[0]->GetNumPoints()*
877  m_base[2]->GetNumPoints();
878  }
879  else
880  {
881  return m_base[1]->GetNumPoints()*
882  m_base[2]->GetNumPoints();
883  }
884  }
885 
886  LibUtilities::PointsKey StdPrismExp::v_GetFacePointsKey(
887  const int i, const int j) const
888  {
889  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
890  ASSERTL2(j == 0 || j == 1, "face direction is out of range");
891 
892  if (i == 0)
893  {
894  return m_base[j]->GetPointsKey();
895  }
896  else if (i == 1 || i == 3)
897  {
898  return m_base[2*j]->GetPointsKey();
899  }
900  else
901  {
902  return m_base[j+1]->GetPointsKey();
903  }
904  }
905 
906  const LibUtilities::BasisKey StdPrismExp::v_DetFaceBasisKey(
907  const int i, const int k) const
908  {
909  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
910  ASSERTL2(k >= 0 && k <= 1, "basis key id is out of range");
911 
912  switch(i)
913  {
914  case 0:
915  {
916  return EvaluateQuadFaceBasisKey(k,
917  m_base[k]->GetBasisType(),
918  m_base[k]->GetNumPoints(),
919  m_base[k]->GetNumModes());
920  }
921  case 2:
922  case 4:
923  {
924  return EvaluateQuadFaceBasisKey(k,
925  m_base[k+1]->GetBasisType(),
926  m_base[k+1]->GetNumPoints(),
927  m_base[k+1]->GetNumModes());
928  }
929  case 1:
930  case 3:
931  {
932  return EvaluateTriFaceBasisKey(k,
933  m_base[2*k]->GetBasisType(),
934  m_base[2*k]->GetNumPoints(),
935  m_base[2*k]->GetNumModes());
936 
937  }
938  break;
939  }
940 
941  // Should never get here.
942  return LibUtilities::NullBasisKey;
943  }
944 
945  int StdPrismExp::v_CalcNumberOfCoefficients(const std::vector<unsigned int> &nummodes,
946  int &modes_offset)
947  {
948  int nmodes = LibUtilities::StdPrismData::getNumberOfCoefficients(
949  nummodes[modes_offset],
950  nummodes[modes_offset+1],
951  nummodes[modes_offset+2]);
952 
953  modes_offset += 3;
954  return nmodes;
955  }
956 
957  LibUtilities::BasisType StdPrismExp::v_GetEdgeBasisType(const int i) const
958  {
959  ASSERTL2(i >= 0 && i <= 8, "edge id is out of range");
960  if (i == 0 || i == 2)
961  {
962  return GetBasisType(0);
963  }
964  else if (i == 1 || i == 3 || i == 8)
965  {
966  return GetBasisType(1);
967  }
968  else
969  {
970  return GetBasisType(2);
971  }
972  }
973 
974  bool StdPrismExp::v_IsBoundaryInteriorExpansion()
975  {
976  return (m_base[0]->GetBasisType() == LibUtilities::eModified_A) &&
977  (m_base[1]->GetBasisType() == LibUtilities::eModified_A) &&
978  (m_base[2]->GetBasisType() == LibUtilities::eModified_B);
979  }
980 
981  //---------------------------------------
982  // Mappings
983  //---------------------------------------
984 
985 
986  void StdPrismExp::v_GetFaceToElementMap(
987  const int fid,
988  const Orientation faceOrient,
989  Array<OneD, unsigned int> &maparray,
990  Array<OneD, int> &signarray,
991  int nummodesA,
992  int nummodesB)
993  {
994  ASSERTL1(GetEdgeBasisType(0) == GetEdgeBasisType(1),
995  "Method only implemented if BasisType is identical"
996  "in x and y directions");
997  ASSERTL1(GetEdgeBasisType(0) == LibUtilities::eModified_A &&
998  GetEdgeBasisType(4) == LibUtilities::eModified_B,
999  "Method only implemented for Modified_A BasisType"
1000  "(x and y direction) and Modified_B BasisType (z "
1001  "direction)");
1002 
1003  int i, j, p, q, r, nFaceCoeffs, idx = 0;
1004 
1005  if (nummodesA == -1)
1006  {
1007  switch (fid)
1008  {
1009  case 0:
1010  nummodesA = m_base[0]->GetNumModes();
1011  nummodesB = m_base[1]->GetNumModes();
1012  break;
1013  case 1:
1014  case 3:
1015  nummodesA = m_base[0]->GetNumModes();
1016  nummodesB = m_base[2]->GetNumModes();
1017  break;
1018  case 2:
1019  case 4:
1020  nummodesA = m_base[1]->GetNumModes();
1021  nummodesB = m_base[2]->GetNumModes();
1022  break;
1023  }
1024  nFaceCoeffs = GetFaceNcoeffs(fid);
1025  }
1026  else if (fid == 1 || fid == 3)
1027  {
1028  nFaceCoeffs = nummodesB + (nummodesA-1)*(1+2*(nummodesB-1)-(nummodesA-1))/2;
1029  }
1030  else
1031  {
1032  nFaceCoeffs = nummodesA*nummodesB;
1033  }
1034 
1035  // Allocate the map array and sign array; set sign array to ones (+)
1036  if (maparray.num_elements() != nFaceCoeffs)
1037  {
1038  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1039  }
1040 
1041  if (signarray.num_elements() != nFaceCoeffs)
1042  {
1043  signarray = Array<OneD, int>(nFaceCoeffs,1);
1044  }
1045  else
1046  {
1047  fill(signarray.get(), signarray.get() + nFaceCoeffs, 1);
1048  }
1049 
1050  // Set up an array indexing for quads, since the ordering may need
1051  // to be transposed.
1052  Array<OneD, int> arrayindx(nFaceCoeffs,-1);
1053 
1054  if (fid != 1 && fid != 3)
1055  {
1056  for (i = 0; i < nummodesB; i++)
1057  {
1058  for (j = 0; j < nummodesA; j++)
1059  {
1060  if (faceOrient < 9)
1061  {
1062  arrayindx[i*nummodesA+j] = i*nummodesA+j;
1063  }
1064  else
1065  {
1066  arrayindx[i*nummodesA+j] = j*nummodesB+i;
1067  }
1068  }
1069  }
1070  }
1071 
1072  // Set up ordering inside each 2D face. Also for triangular faces,
1073  // populate signarray.
1074  switch (fid)
1075  {
1076  case 0: // Bottom quad
1077  for (q = 0; q < nummodesB; ++q)
1078  {
1079  for (p = 0; p < nummodesA; ++p)
1080  {
1081  maparray[arrayindx[q*nummodesA+p]] = GetMode(p,q,0);
1082  }
1083  }
1084  break;
1085 
1086  case 1: // Left triangle
1087  for (p = 0; p < nummodesA; ++p)
1088  {
1089  for (r = 0; r < nummodesB-p; ++r)
1090  {
1091  if ((int)faceOrient == 7 && p > 1)
1092  {
1093  signarray[idx] = p % 2 ? -1 : 1;
1094  }
1095  maparray[idx++] = GetMode(p,0,r);
1096  }
1097  }
1098  break;
1099 
1100  case 2: // Slanted quad
1101  for (q = 0; q < nummodesA; ++q)
1102  {
1103  maparray[arrayindx[q]] = GetMode(1,q,0);
1104  }
1105  for (q = 0; q < nummodesA; ++q)
1106  {
1107  maparray[arrayindx[nummodesA+q]] = GetMode(0,q,1);
1108  }
1109  for (r = 1; r < nummodesB-1; ++r)
1110  {
1111  for (q = 0; q < nummodesA; ++q)
1112  {
1113  maparray[arrayindx[(r+1)*nummodesA+q]] = GetMode(1,q,r);
1114  }
1115  }
1116  break;
1117 
1118  case 3: // Right triangle
1119  for (p = 0; p < nummodesA; ++p)
1120  {
1121  for (r = 0; r < nummodesB-p; ++r)
1122  {
1123  if ((int)faceOrient == 7 && p > 1)
1124  {
1125  signarray[idx] = p % 2 ? -1 : 1;
1126  }
1127  maparray[idx++] = GetMode(p, 1, r);
1128  }
1129  }
1130  break;
1131 
1132  case 4: // Rear quad
1133  for (r = 0; r < nummodesB; ++r)
1134  {
1135  for (q = 0; q < nummodesA; ++q)
1136  {
1137  maparray[arrayindx[r*nummodesA+q]] = GetMode(0, q, r);
1138  }
1139  }
1140  break;
1141 
1142  default:
1143  ASSERTL0(false, "Face to element map unavailable.");
1144  }
1145 
1146  if (fid == 1 || fid == 3)
1147  {
1148  // Triangles only have one possible orientation (base
1149  // direction reversed); swap edge modes.
1150  if ((int)faceOrient == 7)
1151  {
1152  swap(maparray[0], maparray[nummodesA]);
1153  for (i = 1; i < nummodesA-1; ++i)
1154  {
1155  swap(maparray[i+1], maparray[nummodesA+i]);
1156  }
1157  }
1158  }
1159  else
1160  {
1161  // The code below is exactly the same as that taken from
1162  // StdHexExp and reverses the 'b' and 'a' directions as
1163  // appropriate (1st and 2nd if statements respectively) in
1164  // quadrilateral faces.
1165  if (faceOrient == 6 || faceOrient == 8 ||
1166  faceOrient == 11 || faceOrient == 12)
1167  {
1168  if (faceOrient < 9)
1169  {
1170  for (i = 3; i < nummodesB; i += 2)
1171  {
1172  for (j = 0; j < nummodesA; j++)
1173  {
1174  signarray[arrayindx[i*nummodesA+j]] *= -1;
1175  }
1176  }
1177 
1178  for (i = 0; i < nummodesA; i++)
1179  {
1180  swap(maparray [i], maparray [i+nummodesA]);
1181  swap(signarray[i], signarray[i+nummodesA]);
1182  }
1183  }
1184  else
1185  {
1186  for (i = 0; i < nummodesB; i++)
1187  {
1188  for (j = 3; j < nummodesA; j += 2)
1189  {
1190  signarray[arrayindx[i*nummodesA+j]] *= -1;
1191  }
1192  }
1193 
1194  for (i = 0; i < nummodesB; i++)
1195  {
1196  swap (maparray [i], maparray [i+nummodesB]);
1197  swap (signarray[i], signarray[i+nummodesB]);
1198  }
1199  }
1200  }
1201 
1202  if (faceOrient == 7 || faceOrient == 8 ||
1203  faceOrient == 10 || faceOrient == 12)
1204  {
1205  if (faceOrient < 9)
1206  {
1207  for (i = 0; i < nummodesB; i++)
1208  {
1209  for (j = 3; j < nummodesA; j += 2)
1210  {
1211  signarray[arrayindx[i*nummodesA+j]] *= -1;
1212  }
1213  }
1214 
1215  for(i = 0; i < nummodesB; i++)
1216  {
1217  swap(maparray [i*nummodesA], maparray [i*nummodesA+1]);
1218  swap(signarray[i*nummodesA], signarray[i*nummodesA+1]);
1219  }
1220  }
1221  else
1222  {
1223  for (i = 3; i < nummodesB; i += 2)
1224  {
1225  for (j = 0; j < nummodesA; j++)
1226  {
1227  signarray[arrayindx[i*nummodesA+j]] *= -1;
1228  }
1229  }
1230 
1231  for (i = 0; i < nummodesA; i++)
1232  {
1233  swap(maparray [i*nummodesB], maparray [i*nummodesB+1]);
1234  swap(signarray[i*nummodesB], signarray[i*nummodesB+1]);
1235  }
1236  }
1237  }
1238  }
1239  }
1240 
1241  int StdPrismExp::v_GetVertexMap(const int vId, bool useCoeffPacking)
1242  {
1243  ASSERTL0(GetEdgeBasisType(vId) == LibUtilities::eModified_A ||
1244  GetEdgeBasisType(vId) == LibUtilities::eModified_A ||
1245  GetEdgeBasisType(vId) == LibUtilities::eModified_B,
1246  "Mapping not defined for this type of basis");
1247 
1248  int l = 0;
1249 
1250  if(useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1251  {
1252  switch (vId)
1253  {
1254  case 0:
1255  l = GetMode(0,0,0);
1256  break;
1257  case 1:
1258  l = GetMode(0,0,1);
1259  break;
1260  case 2:
1261  l = GetMode(0,1,0);
1262  break;
1263  case 3:
1264  l = GetMode(0,1,1);
1265  break;
1266  case 4:
1267  l = GetMode(1,0,0);
1268  break;
1269  case 5:
1270  l = GetMode(1,1,0);
1271  break;
1272  default:
1273  ASSERTL0(false, "local vertex id must be between 0 and 5");
1274  }
1275  }
1276  else
1277  {
1278  switch (vId)
1279  {
1280  case 0:
1281  l = GetMode(0,0,0);
1282  break;
1283  case 1:
1284  l = GetMode(1,0,0);
1285  break;
1286  case 2:
1287  l = GetMode(1,1,0);
1288  break;
1289  case 3:
1290  l = GetMode(0,1,0);
1291  break;
1292  case 4:
1293  l = GetMode(0,0,1);
1294  break;
1295  case 5:
1296  l = GetMode(0,1,1);
1297  break;
1298  default:
1299  ASSERTL0(false, "local vertex id must be between 0 and 5");
1300  }
1301  }
1302 
1303  return l;
1304  }
1305 
1306  void StdPrismExp::v_GetEdgeInteriorMap(
1307  const int eid,
1308  const Orientation edgeOrient,
1309  Array<OneD, unsigned int> &maparray,
1310  Array<OneD, int> &signarray)
1311  {
1312  int i;
1313  bool signChange;
1314  const int P = m_base[0]->GetNumModes() - 1;
1315  const int Q = m_base[1]->GetNumModes() - 1;
1316  const int R = m_base[2]->GetNumModes() - 1;
1317  const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid) - 2;
1318 
1319  if (maparray.num_elements() != nEdgeIntCoeffs)
1320  {
1321  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1322  }
1323 
1324  if(signarray.num_elements() != nEdgeIntCoeffs)
1325  {
1326  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1327  }
1328  else
1329  {
1330  fill(signarray.get(), signarray.get()+nEdgeIntCoeffs, 1);
1331  }
1332 
1333  // If edge is oriented backwards, change sign of modes which have
1334  // degree 2n+1, n >= 1.
1335  signChange = edgeOrient == eBackwards;
1336 
1337  switch (eid)
1338  {
1339  case 0:
1340  for (i = 2; i <= P; ++i)
1341  {
1342  maparray[i-2] = GetMode(i,0,0);
1343  }
1344  break;
1345 
1346  case 1:
1347  for (i = 2; i <= Q; ++i)
1348  {
1349  maparray[i-2] = GetMode(1,i,0);
1350  }
1351  break;
1352 
1353  case 2:
1354  // Base quad; reverse direction.
1355  //signChange = !signChange;
1356  for (i = 2; i <= P; ++i)
1357  {
1358  maparray[i-2] = GetMode(i,1,0);
1359  }
1360  break;
1361 
1362  case 3:
1363  // Base quad; reverse direction.
1364  //signChange = !signChange;
1365  for (i = 2; i <= Q; ++i)
1366  {
1367  maparray[i-2] = GetMode(0,i,0);
1368  }
1369  break;
1370 
1371  case 4:
1372  for (i = 2; i <= R; ++i)
1373  {
1374  maparray[i-2] = GetMode(0,0,i);
1375  }
1376  break;
1377 
1378  case 5:
1379  for (i = 1; i <= R-1; ++i)
1380  {
1381  maparray[i-1] = GetMode(1,0,i);
1382  }
1383  break;
1384 
1385  case 6:
1386  for (i = 1; i <= R-1; ++i)
1387  {
1388  maparray[i-1] = GetMode(1,1,i);
1389  }
1390  break;
1391 
1392  case 7:
1393  for (i = 2; i <= R; ++i)
1394  {
1395  maparray[i-2] = GetMode(0,1,i);
1396  }
1397  break;
1398 
1399  case 8:
1400  for (i = 2; i <= Q; ++i)
1401  {
1402  maparray[i-2] = GetMode(0,i,1);
1403  }
1404  break;
1405 
1406  default:
1407  ASSERTL0(false, "Edge not defined.");
1408  break;
1409  }
1410 
1411  if (signChange)
1412  {
1413  for (i = 1; i < nEdgeIntCoeffs; i += 2)
1414  {
1415  signarray[i] = -1;
1416  }
1417  }
1418  }
1419 
1420  void StdPrismExp::v_GetFaceInteriorMap(
1421  const int fid,
1422  const Orientation faceOrient,
1423  Array<OneD, unsigned int> &maparray,
1424  Array<OneD, int> &signarray)
1425  {
1426  const int P = m_base[0]->GetNumModes() - 1;
1427  const int Q = m_base[1]->GetNumModes() - 1;
1428  const int R = m_base[2]->GetNumModes() - 1;
1429  const int nFaceIntCoeffs = v_GetFaceIntNcoeffs(fid);
1430  int p, q, r, idx = 0;
1431  int nummodesA = 0;
1432  int nummodesB = 0;
1433  int i = 0;
1434  int j = 0;
1435 
1436  if (maparray.num_elements() != nFaceIntCoeffs)
1437  {
1438  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1439  }
1440 
1441  if (signarray.num_elements() != nFaceIntCoeffs)
1442  {
1443  signarray = Array<OneD, int>(nFaceIntCoeffs, 1);
1444  }
1445  else
1446  {
1447  fill(signarray.get(), signarray.get()+nFaceIntCoeffs, 1);
1448  }
1449 
1450  // Set up an array indexing for quad faces, since the ordering may
1451  // need to be transposed depending on orientation.
1452  Array<OneD, int> arrayindx(nFaceIntCoeffs);
1453  if (fid != 1 && fid != 3)
1454  {
1455  if (fid == 0) // Base quad
1456  {
1457  nummodesA = P-1;
1458  nummodesB = Q-1;
1459  }
1460  else // front and back quad
1461  {
1462  nummodesA = Q-1;
1463  nummodesB = R-1;
1464  }
1465 
1466  for (i = 0; i < nummodesB; i++)
1467  {
1468  for (j = 0; j < nummodesA; j++)
1469  {
1470  if (faceOrient < 9)
1471  {
1472  arrayindx[i*nummodesA+j] = i*nummodesA+j;
1473  }
1474  else
1475  {
1476  arrayindx[i*nummodesA+j] = j*nummodesB+i;
1477  }
1478  }
1479  }
1480  }
1481 
1482  switch (fid)
1483  {
1484  case 0: // Bottom quad
1485  for (q = 2; q <= Q; ++q)
1486  {
1487  for (p = 2; p <= P; ++p)
1488  {
1489  maparray[arrayindx[(q-2)*nummodesA+(p-2)]] = GetMode(p,q,0);
1490  }
1491  }
1492  break;
1493 
1494  case 1: // Left triangle
1495  for (p = 2; p <= P; ++p)
1496  {
1497  for (r = 1; r <= R-p; ++r)
1498  {
1499  if ((int)faceOrient == 7)
1500  {
1501  signarray[idx] = p % 2 ? -1 : 1;
1502  }
1503  maparray[idx++] = GetMode(p,0,r);
1504  }
1505  }
1506  break;
1507 
1508  case 2: // Slanted quad
1509  for (r = 1; r <= R-1; ++r)
1510  {
1511  for (q = 2; q <= Q; ++q)
1512  {
1513  maparray[arrayindx[(r-1)*nummodesA+(q-2)]] = GetMode(1, q, r);
1514  }
1515  }
1516  break;
1517 
1518  case 3: // Right triangle
1519  for (p = 2; p <= P; ++p)
1520  {
1521  for (r = 1; r <= R-p; ++r)
1522  {
1523  if ((int)faceOrient == 7)
1524  {
1525  signarray[idx] = p % 2 ? -1 : 1;
1526  }
1527  maparray[idx++] = GetMode(p, 1, r);
1528  }
1529  }
1530  break;
1531 
1532  case 4: // Back quad
1533  for (r = 2; r <= R; ++r)
1534  {
1535  for (q = 2; q <= Q; ++q)
1536  {
1537  maparray[arrayindx[(r-2)*nummodesA+(q-2)]] = GetMode(0, q, r);
1538  }
1539  }
1540  break;
1541 
1542  default:
1543  ASSERTL0(false, "Face interior map not available.");
1544  }
1545 
1546  // Triangular faces are processed in the above switch loop; for
1547  // remaining quad faces, set up orientation if necessary.
1548  if (fid == 1 || fid == 3)
1549  return;
1550 
1551  if (faceOrient == 6 || faceOrient == 8 ||
1552  faceOrient == 11 || faceOrient == 12)
1553  {
1554  if (faceOrient < 9)
1555  {
1556  for (i = 1; i < nummodesB; i += 2)
1557  {
1558  for (j = 0; j < nummodesA; j++)
1559  {
1560  signarray[arrayindx[i*nummodesA+j]] *= -1;
1561  }
1562  }
1563  }
1564  else
1565  {
1566  for (i = 0; i < nummodesB; i++)
1567  {
1568  for (j = 1; j < nummodesA; j += 2)
1569  {
1570  signarray[arrayindx[i*nummodesA+j]] *= -1;
1571  }
1572  }
1573  }
1574  }
1575 
1576  if (faceOrient == 7 || faceOrient == 8 ||
1577  faceOrient == 10 || faceOrient == 12)
1578  {
1579  if (faceOrient < 9)
1580  {
1581  for (i = 0; i < nummodesB; i++)
1582  {
1583  for (j = 1; j < nummodesA; j += 2)
1584  {
1585  signarray[arrayindx[i*nummodesA+j]] *= -1;
1586  }
1587  }
1588  }
1589  else
1590  {
1591  for (i = 1; i < nummodesB; i += 2)
1592  {
1593  for (j = 0; j < nummodesA; j++)
1594  {
1595  signarray[arrayindx[i*nummodesA+j]] *= -1;
1596  }
1597  }
1598  }
1599  }
1600  }
1601 
1602  void StdPrismExp::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
1603  {
1604  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1605  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1606  "BasisType is not a boundary interior form");
1607  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1608  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1609  "BasisType is not a boundary interior form");
1610  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_B ||
1611  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1612  "BasisType is not a boundary interior form");
1613 
1614  int P = m_base[0]->GetNumModes() - 1, p;
1615  int Q = m_base[1]->GetNumModes() - 1, q;
1616  int R = m_base[2]->GetNumModes() - 1, r;
1617 
1618  int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1619 
1620  if(outarray.num_elements()!=nIntCoeffs)
1621  {
1622  outarray = Array<OneD, unsigned int>(nIntCoeffs);
1623  }
1624 
1625  int idx = 0;
1626 
1627  // Loop over all interior modes.
1628  for (p = 2; p <= P; ++p)
1629  {
1630  for (q = 2; q <= Q; ++q)
1631  {
1632  for (r = 1; r <= R-p; ++r)
1633  {
1634  outarray[idx++] = GetMode(p,q,r);
1635  }
1636  }
1637  }
1638  }
1639 
1640  void StdPrismExp::v_GetBoundaryMap(Array<OneD, unsigned int> &maparray)
1641  {
1642  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1643  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1644  "BasisType is not a boundary interior form");
1645  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1646  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1647  "BasisType is not a boundary interior form");
1648  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_B ||
1649  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1650  "BasisType is not a boundary interior form");
1651 
1652  int P = m_base[0]->GetNumModes() - 1, p;
1653  int Q = m_base[1]->GetNumModes() - 1, q;
1654  int R = m_base[2]->GetNumModes() - 1, r;
1655  int idx = 0;
1656 
1657  // Loop over all boundary modes (in ascending order).
1658  for (p = 0; p <= P; ++p)
1659  {
1660  // First two q-r planes are entirely boundary modes.
1661  if (p <= 1)
1662  {
1663  for (q = 0; q <= Q; ++q)
1664  {
1665  for (r = 0; r <= R-p; ++r)
1666  {
1667  maparray[idx++] = GetMode(p,q,r);
1668  }
1669  }
1670  }
1671  else
1672  {
1673  // Remaining q-r planes contain boundary modes on the two
1674  // left-hand sides and bottom edge.
1675  for (q = 0; q <= Q; ++q)
1676  {
1677  if (q <= 1)
1678  {
1679  for (r = 0; r <= R-p; ++r)
1680  {
1681  maparray[idx++] = GetMode(p,q,r);
1682  }
1683  }
1684  else
1685  {
1686  maparray[idx++] = GetMode(p,q,0);
1687  }
1688  }
1689  }
1690  }
1691  }
1692 
1693 
1694 
1695  //---------------------------------------
1696  // Wrapper functions
1697  //---------------------------------------
1698 
1699  DNekMatSharedPtr StdPrismExp::v_GenMatrix(const StdMatrixKey &mkey)
1700  {
1701 
1702  MatrixType mtype = mkey.GetMatrixType();
1703 
1704  DNekMatSharedPtr Mat;
1705 
1706  switch(mtype)
1707  {
1708  case ePhysInterpToEquiSpaced:
1709  {
1710  int nq0 = m_base[0]->GetNumPoints();
1711  int nq1 = m_base[1]->GetNumPoints();
1712  int nq2 = m_base[2]->GetNumPoints();
1713  int nq = max(nq0,max(nq1,nq2));
1714  int neq = LibUtilities::StdPrismData::
1715  getNumberOfCoefficients (nq,nq,nq);
1716  Array<OneD, Array<OneD, NekDouble> > coords (neq);
1717  Array<OneD, NekDouble> coll (3);
1718  Array<OneD, DNekMatSharedPtr> I (3);
1719  Array<OneD, NekDouble> tmp (nq0);
1720 
1721  Mat = MemoryManager<DNekMat>::
1722  AllocateSharedPtr(neq,nq0*nq1*nq2);
1723  int cnt = 0;
1724  for(int i = 0; i < nq; ++i)
1725  {
1726  for(int j = 0; j < nq; ++j)
1727  {
1728  for(int k = 0; k < nq-i; ++k,++cnt)
1729  {
1730  coords[cnt] = Array<OneD, NekDouble>(3);
1731  coords[cnt][0] = -1.0 + 2*k/(NekDouble)(nq-1);
1732  coords[cnt][1] = -1.0 + 2*j/(NekDouble)(nq-1);
1733  coords[cnt][2] = -1.0 + 2*i/(NekDouble)(nq-1);
1734  }
1735  }
1736  }
1737 
1738  for(int i = 0; i < neq; ++i)
1739  {
1740  LocCoordToLocCollapsed(coords[i],coll);
1741 
1742  I[0] = m_base[0]->GetI(coll );
1743  I[1] = m_base[1]->GetI(coll+1);
1744  I[2] = m_base[2]->GetI(coll+2);
1745 
1746  // interpolate first coordinate direction
1747  NekDouble fac;
1748  for( int k = 0; k < nq2; ++k)
1749  {
1750  for (int j = 0; j < nq1; ++j)
1751  {
1752 
1753  fac = (I[1]->GetPtr())[j]*(I[2]->GetPtr())[k];
1754  Vmath::Smul(nq0,fac,I[0]->GetPtr(),1,tmp,1);
1755 
1756  Vmath::Vcopy(nq0, &tmp[0], 1,
1757  Mat->GetRawPtr() +
1758  k * nq0 * nq1 * neq +
1759  j * nq0 * neq + i,
1760  neq);
1761  }
1762  }
1763  }
1764  }
1765  break;
1766  default:
1767  {
1768  Mat = StdExpansion::CreateGeneralMatrix(mkey);
1769  }
1770  break;
1771  }
1772 
1773  return Mat;
1774  }
1775 
1776  DNekMatSharedPtr StdPrismExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
1777  {
1778  return v_GenMatrix(mkey);
1779  }
1780 
1781 
1782 
1783  /**
1784  * @brief Compute the local mode number in the expansion for a
1785  * particular tensorial combination.
1786  *
1787  * Modes are numbered with the r index travelling fastest, followed by
1788  * q and then p, and each q-r plane is of size (R+1-p). For example,
1789  * with P=1, Q=2, R=3, the indexing inside each q-r plane (with r
1790  * increasing upwards and q to the right) is:
1791  *
1792  * p = 0: p = 1:
1793  * -----------------------
1794  * 3 7 11
1795  * 2 6 10 14 17 20
1796  * 1 5 9 13 16 19
1797  * 0 4 8 12 15 18
1798  *
1799  * Note that in this element, we must have that \f$ P <= R \f$.
1800  */
1801  int StdPrismExp::GetMode(int p, int q, int r)
1802  {
1803  int Q = m_base[1]->GetNumModes() - 1;
1804  int R = m_base[2]->GetNumModes() - 1;
1805 
1806  return r + // Skip along stacks (r-direction)
1807  q*(R+1-p) + // Skip along columns (q-direction)
1808  (Q+1)*(p*R + 1-(p-2)*(p-1)/2); // Skip along rows (p-direction)
1809  }
1810 
1811  void StdPrismExp::v_MultiplyByStdQuadratureMetric(
1812  const Array<OneD, const NekDouble>& inarray,
1813  Array<OneD, NekDouble>& outarray)
1814  {
1815  int i, j;
1816  int nquad0 = m_base[0]->GetNumPoints();
1817  int nquad1 = m_base[1]->GetNumPoints();
1818  int nquad2 = m_base[2]->GetNumPoints();
1819 
1820  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1821  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1822  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
1823 
1824  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1825 
1826  // Multiply by integration constants in x-direction
1827  for(i = 0; i < nquad1*nquad2; ++i)
1828  {
1829  Vmath::Vmul(nquad0, inarray.get()+i*nquad0, 1,
1830  w0.get(), 1, outarray.get()+i*nquad0,1);
1831  }
1832 
1833  // Multiply by integration constants in y-direction
1834  for(j = 0; j < nquad2; ++j)
1835  {
1836  for(i = 0; i < nquad1; ++i)
1837  {
1838  Blas::Dscal(nquad0,w1[i], &outarray[0]+i*nquad0 +
1839  j*nquad0*nquad1,1);
1840  }
1841  }
1842 
1843  // Multiply by integration constants in z-direction; need to
1844  // incorporate factor (1-eta_3)/2 into weights, but only if using
1845  // GLL quadrature points.
1846  switch(m_base[2]->GetPointsType())
1847  {
1848  // (1,0) Jacobi inner product.
1849  case LibUtilities::eGaussRadauMAlpha1Beta0:
1850  for(i = 0; i < nquad2; ++i)
1851  {
1852  Blas::Dscal(nquad0*nquad1, 0.5*w2[i],
1853  &outarray[0]+i*nquad0*nquad1, 1);
1854  }
1855  break;
1856 
1857  default:
1858  for(i = 0; i < nquad2; ++i)
1859  {
1860  Blas::Dscal(nquad0*nquad1,0.25*(1-z2[i])*w2[i],
1861  &outarray[0]+i*nquad0*nquad1,1);
1862  }
1863  break;
1864  }
1865 
1866  }
1867 
1868  void StdPrismExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
1869  const StdMatrixKey &mkey)
1870  {
1871  // Generate an orthonogal expansion
1872  int qa = m_base[0]->GetNumPoints();
1873  int qb = m_base[1]->GetNumPoints();
1874  int qc = m_base[2]->GetNumPoints();
1875  int nmodes_a = m_base[0]->GetNumModes();
1876  int nmodes_b = m_base[1]->GetNumModes();
1877  int nmodes_c = m_base[2]->GetNumModes();
1878  // Declare orthogonal basis.
1879  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
1880  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
1881  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
1882 
1883  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
1884  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A,nmodes_b,pb);
1885  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_B,nmodes_c,pc);
1886  StdPrismExp OrthoExp(Ba,Bb,Bc);
1887 
1888  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1889  int i,j,k,cnt = 0;
1890 
1891  //SVV filter paramaters (how much added diffusion relative to physical one
1892  // and fraction of modes from which you start applying this added diffusion)
1893  //
1894  NekDouble SvvDiffCoeff = mkey.GetConstFactor(StdRegions::eFactorSVVDiffCoeff);
1895  NekDouble SVVCutOff = mkey.GetConstFactor(StdRegions::eFactorSVVCutoffRatio);
1896 
1897  //Defining the cut of mode
1898  int cutoff_a = (int) (SVVCutOff*nmodes_a);
1899  int cutoff_b = (int) (SVVCutOff*nmodes_b);
1900  int cutoff_c = (int) (SVVCutOff*nmodes_c);
1901  //To avoid the fac[j] from blowing up
1902  NekDouble epsilon = 1;
1903 
1904  // project onto modal space.
1905  OrthoExp.FwdTrans(array,orthocoeffs);
1906  int nmodes = min(min(nmodes_a,nmodes_b),nmodes_c);
1907  NekDouble cutoff = min(min(cutoff_a,cutoff_b),cutoff_c);
1908 
1909  //------"New" Version August 22nd '13--------------------
1910  for(i = 0; i < nmodes_a; ++i)//P
1911  {
1912  for(j = 0; j < nmodes_b; ++j) //Q
1913  {
1914  for(k = 0; k < nmodes_c-i; ++k) //R
1915  {
1916  if(j >= cutoff || i + k >= cutoff)
1917  {
1918  orthocoeffs[cnt] *= (1.0+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)))));
1919  }
1920  cnt++;
1921  }
1922  }
1923  }
1924 
1925  // backward transform to physical space
1926  OrthoExp.BwdTrans(orthocoeffs,array);
1927  }
1928 
1929 
1930 
1931  void StdPrismExp::v_ReduceOrderCoeffs(
1932  int numMin,
1933  const Array<OneD, const NekDouble> &inarray,
1934  Array<OneD, NekDouble> &outarray)
1935  {
1936  int nquad0 = m_base[0]->GetNumPoints();
1937  int nquad1 = m_base[1]->GetNumPoints();
1938  int nquad2 = m_base[2]->GetNumPoints();
1939  int nqtot = nquad0*nquad1*nquad2;
1940  int nmodes0 = m_base[0]->GetNumModes();
1941  int nmodes1 = m_base[1]->GetNumModes();
1942  int nmodes2 = m_base[2]->GetNumModes();
1943  int numMax = nmodes0;
1944 
1945  Array<OneD, NekDouble> coeff (m_ncoeffs);
1946  Array<OneD, NekDouble> coeff_tmp1(m_ncoeffs, 0.0);
1947  Array<OneD, NekDouble> phys_tmp (nqtot, 0.0);
1948  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
1949 
1950 
1951  const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
1952  const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
1953  const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
1954 
1955  LibUtilities::BasisKey bortho0(
1956  LibUtilities::eOrtho_A, nmodes0, Pkey0);
1957  LibUtilities::BasisKey bortho1(
1958  LibUtilities::eOrtho_A, nmodes1, Pkey1);
1959  LibUtilities::BasisKey bortho2(
1960  LibUtilities::eOrtho_B, nmodes2, Pkey2);
1961 
1962  int cnt = 0;
1963  int u = 0;
1964  int i = 0;
1965  StdRegions::StdPrismExpSharedPtr OrthoPrismExp;
1966 
1967  OrthoPrismExp = MemoryManager<StdRegions::StdPrismExp>
1968  ::AllocateSharedPtr(bortho0, bortho1, bortho2);
1969 
1970  BwdTrans(inarray,phys_tmp);
1971  OrthoPrismExp->FwdTrans(phys_tmp, coeff);
1972 
1973  // filtering
1974  for (u = 0; u < numMin; ++u)
1975  {
1976  for (i = 0; i < numMin; ++i)
1977  {
1978  Vmath::Vcopy(numMin - u, tmp = coeff + cnt, 1,
1979  tmp2 = coeff_tmp1 + cnt, 1);
1980  cnt += numMax - u;
1981  }
1982 
1983  for (i = numMin; i < numMax; ++i)
1984  {
1985  cnt += numMax - u;
1986  }
1987  }
1988 
1989  OrthoPrismExp->BwdTrans(coeff_tmp1, phys_tmp);
1990  StdPrismExp::FwdTrans(phys_tmp, outarray);
1991  }
1992  }//end namespace
1993 }//end namespace
Generated on Wed Mar 4 2015 16:33:17 for Nektar++ by   doxygen 1.8.1.2