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