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