Nektar++
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_dxi0 = 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_dxi2 =
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  int i, j, mode, mode1, cnt;
379 
380  // Perform summation over '2' direction
381  mode = mode1 = cnt = 0;
382  for(i = 0; i < order0; ++i)
383  {
384  for(j = 0; j < order1-i; ++j, ++cnt)
385  {
386  Blas::Dgemv('N', nquad2, order2-i-j,
387  1.0, base2.get()+mode*nquad2, nquad2,
388  inarray.get()+mode1, 1,
389  0.0, tmp.get()+cnt*nquad2, 1);
390  mode += order2-i-j;
391  mode1 += order2-i-j;
392  }
393  //increment mode in case order1!=order2
394  for(j = order1-i; j < order2-i; ++j)
395  {
396  mode += order2-i-j;
397  }
398  }
399 
400  // fix for modified basis by adding split of top singular
401  // vertex mode - currently (1+c)/2 x (1-b)/2 x (1-a)/2
402  // component is evaluated
403  if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
404  {
405  // top singular vertex - (1+c)/2 x (1+b)/2 x (1-a)/2 component
406  Blas::Daxpy(nquad2,inarray[1],base2.get()+nquad2,1,
407  &tmp[0]+nquad2,1);
408 
409  // top singular vertex - (1+c)/2 x (1-b)/2 x (1+a)/2 component
410  Blas::Daxpy(nquad2,inarray[1],base2.get()+nquad2,1,
411  &tmp[0]+order1*nquad2,1);
412  }
413 
414  // Perform summation over '1' direction
415  mode = 0;
416  for(i = 0; i < order0; ++i)
417  {
418  Blas::Dgemm('N', 'T', nquad1, nquad2, order1-i,
419  1.0, base1.get()+mode*nquad1, nquad1,
420  tmp.get()+mode*nquad2, nquad2,
421  0.0, tmp1.get()+i*nquad1*nquad2, nquad1);
422  mode += order1-i;
423  }
424 
425  // fix for modified basis by adding additional split of
426  // top and base singular vertex modes as well as singular
427  // edge
428  if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
429  {
430  // use tmp to sort out singular vertices and
431  // singular edge components with (1+b)/2 (1+a)/2 form
432  for(i = 0; i < nquad2; ++i)
433  {
434  Blas::Daxpy(nquad1,tmp[nquad2+i], base1.get()+nquad1,1,
435  &tmp1[nquad1*nquad2]+i*nquad1,1);
436  }
437  }
438 
439  // Perform summation over '0' direction
440  Blas::Dgemm('N', 'T', nquad0, nquad1*nquad2, order0,
441  1.0, base0.get(), nquad0,
442  tmp1.get(), nquad1*nquad2,
443  0.0, outarray.get(), nquad0);
444  }
445 
446 
447  /**
448  * @param inarray array of physical quadrature points to be
449  * transformed.
450  * @param outarray updated array of expansion coefficients.
451  */
452  void StdTetExp::v_FwdTrans(const Array<OneD, const NekDouble>& inarray,
453  Array<OneD, NekDouble> &outarray)
454  {
455  v_IProductWRTBase(inarray,outarray);
456 
457  // get Mass matrix inverse
458  StdMatrixKey masskey(eInvMass,DetShapeType(),*this);
459  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
460 
461  // copy inarray in case inarray == outarray
462  DNekVec in (m_ncoeffs, outarray);
463  DNekVec out(m_ncoeffs, outarray, eWrapper);
464 
465  out = (*matsys)*in;
466  }
467 
468 
469  //---------------------------------------
470  // Inner product functions
471  //---------------------------------------
472 
473  /**
474  * \f$ \begin{array}{rcl} I_{pqr} = (\phi_{pqr}, u)_{\delta} & = &
475  * \sum_{i=0}^{nq_0} \sum_{j=0}^{nq_1} \sum_{k=0}^{nq_2} \psi_{p}^{a}
476  * (\eta_{1i}) \psi_{pq}^{b} (\eta_{2j}) \psi_{pqr}^{c} (\eta_{3k})
477  * w_i w_j w_k u(\eta_{1,i} \eta_{2,j} \eta_{3,k}) J_{i,j,k}\\ & = &
478  * \sum_{i=0}^{nq_0} \psi_p^a(\eta_{1,i}) \sum_{j=0}^{nq_1}
479  * \psi_{pq}^b(\eta_{2,j}) \sum_{k=0}^{nq_2} \psi_{pqr}^c
480  * u(\eta_{1i},\eta_{2j},\eta_{3k}) J_{i,j,k} \end{array} \f$ \n
481  *
482  * where
483  *
484  * \f$ \phi_{pqr} (\xi_1 , \xi_2 , \xi_3) = \psi_p^a (\eta_1)
485  * \psi_{pq}^b (\eta_2) \psi_{pqr}^c (\eta_3) \f$
486  *
487  * which can be implemented as \n \f$f_{pqr} (\xi_{3k}) =
488  * \sum_{k=0}^{nq_3} \psi_{pqr}^c u(\eta_{1i},\eta_{2j},\eta_{3k})
489  *
490  * J_{i,j,k} = {\bf B_3 U} \f$ \n
491  *
492  * \f$ g_{pq} (\xi_{3k}) = \sum_{j=0}^{nq_1} \psi_{pq}^b (\xi_{2j})
493  * f_{pqr} (\xi_{3k}) = {\bf B_2 F} \f$ \n
494  *
495  * \f$ (\phi_{pqr}, u)_{\delta} = \sum_{k=0}^{nq_0} \psi_{p}^a
496  * (\xi_{3k}) g_{pq} (\xi_{3k}) = {\bf B_1 G} \f$
497  *
498  * @param inarray Function evaluated at physical collocation
499  * points.
500  * @param outarray Inner product with respect to each basis
501  * function over the element.
502  */
503  void StdTetExp::v_IProductWRTBase(
504  const Array<OneD, const NekDouble>& inarray,
505  Array<OneD, NekDouble> & outarray)
506  {
507  ASSERTL1( (m_base[1]->GetBasisType() != LibUtilities::eOrtho_B) ||
508  (m_base[1]->GetBasisType() != LibUtilities::eModified_B),
509  "Basis[1] is not a general tensor type");
510 
511  ASSERTL1( (m_base[2]->GetBasisType() != LibUtilities::eOrtho_C) ||
512  (m_base[2]->GetBasisType() != LibUtilities::eModified_C),
513  "Basis[2] is not a general tensor type");
514 
515  if(m_base[0]->Collocation() && m_base[1]->Collocation())
516  {
517  MultiplyByQuadratureMetric(inarray,outarray);
518  }
519  else
520  {
521  StdTetExp::v_IProductWRTBase_SumFac(inarray,outarray);
522  }
523  }
524 
525 
526  void StdTetExp::v_IProductWRTBase_MatOp(
527  const Array<OneD, const NekDouble>& inarray,
528  Array<OneD, NekDouble>& outarray)
529  {
530  int nq = GetTotPoints();
531  StdMatrixKey iprodmatkey(eIProductWRTBase,DetShapeType(),*this);
532  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
533 
534  Blas::Dgemv('N',m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
535  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
536  }
537 
538 
539  /**
540  * @param inarray Function evaluated at physical collocation
541  * points.
542  * @param outarray Inner product with respect to each basis
543  * function over the element.
544  */
545  void StdTetExp::v_IProductWRTBase_SumFac(
546  const Array<OneD, const NekDouble>& inarray,
547  Array<OneD, NekDouble>& outarray,
548  bool multiplybyweights)
549  {
550  int nquad0 = m_base[0]->GetNumPoints();
551  int nquad1 = m_base[1]->GetNumPoints();
552  int nquad2 = m_base[2]->GetNumPoints();
553  int order0 = m_base[0]->GetNumModes();
554  int order1 = m_base[1]->GetNumModes();
555 
556  Array<OneD, NekDouble> wsp (nquad1*nquad2*order0 +
557  nquad2*order0*(2*order1-order0+1)/2);
558 
559  if(multiplybyweights)
560  {
561  Array<OneD, NekDouble> tmp (nquad0*nquad1*nquad2);
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  else
571  {
572  StdTetExp::IProductWRTBase_SumFacKernel(
573  m_base[0]->GetBdata(),
574  m_base[1]->GetBdata(),
575  m_base[2]->GetBdata(),
576  inarray, outarray, wsp, true, true, true);
577  }
578  }
579 
580 
581  void StdTetExp::v_IProductWRTBase_SumFacKernel(
582  const Array<OneD, const NekDouble>& base0,
583  const Array<OneD, const NekDouble>& base1,
584  const Array<OneD, const NekDouble>& base2,
585  const Array<OneD, const NekDouble>& inarray,
586  Array<OneD, NekDouble> &outarray,
587  Array<OneD, NekDouble> &wsp,
588  bool doCheckCollDir0,
589  bool doCheckCollDir1,
590  bool doCheckCollDir2)
591  {
592  int nquad0 = m_base[0]->GetNumPoints();
593  int nquad1 = m_base[1]->GetNumPoints();
594  int nquad2 = m_base[2]->GetNumPoints();
595 
596  int order0 = m_base[0]->GetNumModes();
597  int order1 = m_base[1]->GetNumModes();
598  int order2 = m_base[2]->GetNumModes();
599 
600  Array<OneD, NekDouble > tmp1 = wsp;
601  Array<OneD, NekDouble > tmp2 = wsp + nquad1*nquad2*order0;
602 
603  int i,j, mode,mode1, cnt;
604 
605  // Inner product with respect to the '0' direction
606  Blas::Dgemm('T', 'N', nquad1*nquad2, order0, nquad0,
607  1.0, inarray.get(), nquad0,
608  base0.get(), nquad0,
609  0.0, tmp1.get(), nquad1*nquad2);
610 
611  // Inner product with respect to the '1' direction
612  for(mode=i=0; i < order0; ++i)
613  {
614  Blas::Dgemm('T', 'N', nquad2, order1-i, nquad1,
615  1.0, tmp1.get()+i*nquad1*nquad2, nquad1,
616  base1.get()+mode*nquad1, nquad1,
617  0.0, tmp2.get()+mode*nquad2, nquad2);
618  mode += order1-i;
619  }
620 
621  // fix for modified basis for base singular vertex
622  if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
623  {
624  //base singular vertex and singular edge (1+b)/2
625  //(1+a)/2 components (makes tmp[nquad2] entry into (1+b)/2)
626  Blas::Dgemv('T', nquad1, nquad2,
627  1.0, tmp1.get()+nquad1*nquad2, nquad1,
628  base1.get()+nquad1, 1,
629  1.0, tmp2.get()+nquad2, 1);
630  }
631 
632  // Inner product with respect to the '2' direction
633  mode = mode1 = cnt = 0;
634  for(i = 0; i < order0; ++i)
635  {
636  for(j = 0; j < order1-i; ++j, ++cnt)
637  {
638  Blas::Dgemv('T', nquad2, order2-i-j,
639  1.0, base2.get()+mode*nquad2, nquad2,
640  tmp2.get()+cnt*nquad2, 1,
641  0.0, outarray.get()+mode1, 1);
642  mode += order2-i-j;
643  mode1 += order2-i-j;
644  }
645  //increment mode in case order1!=order2
646  for(j = order1-i; j < order2-i; ++j)
647  {
648  mode += order2-i-j;
649  }
650  }
651 
652  // fix for modified basis for top singular vertex component
653  // Already have evaluated (1+c)/2 (1-b)/2 (1-a)/2
654  if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
655  {
656  // add in (1+c)/2 (1+b)/2 component
657  outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
658  &tmp2[nquad2],1);
659 
660  // add in (1+c)/2 (1-b)/2 (1+a)/2 component
661  outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
662  &tmp2[nquad2*order1],1);
663  }
664  }
665 
666 
667  void StdTetExp::v_IProductWRTDerivBase(
668  const int dir,
669  const Array<OneD, const NekDouble>& inarray,
670  Array<OneD, NekDouble>& outarray)
671  {
672  StdTetExp::v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
673  }
674 
675 
676  void StdTetExp::v_IProductWRTDerivBase_MatOp(
677  const int dir,
678  const Array<OneD, const NekDouble>& inarray,
679  Array<OneD, NekDouble>& outarray)
680  {
681  ASSERTL0((dir==0)||(dir==1)||(dir==2),"input dir is out of range");
682 
683  int nq = GetTotPoints();
684  MatrixType mtype;
685 
686  switch (dir)
687  {
688  case 0:
689  mtype = eIProductWRTDerivBase0;
690  break;
691  case 1:
692  mtype = eIProductWRTDerivBase1;
693  break;
694  case 2:
695  mtype = eIProductWRTDerivBase2;
696  break;
697  }
698 
699  StdMatrixKey iprodmatkey(mtype,DetShapeType(),*this);
700  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
701 
702  Blas::Dgemv('N',m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
703  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
704  }
705 
706 
707  /**
708  * @param inarray Function evaluated at physical collocation
709  * points.
710  * @param outarray Inner product with respect to each basis
711  * function over the element.
712  */
713  void StdTetExp::v_IProductWRTDerivBase_SumFac(
714  const int dir,
715  const Array<OneD, const NekDouble>& inarray,
716  Array<OneD, NekDouble>& outarray)
717  {
718  int i;
719  int nquad0 = m_base[0]->GetNumPoints();
720  int nquad1 = m_base[1]->GetNumPoints();
721  int nquad2 = m_base[2]->GetNumPoints();
722  int nqtot = nquad0*nquad1*nquad2;
723  int nmodes0 = m_base[0]->GetNumModes();
724  int nmodes1 = m_base[1]->GetNumModes();
725  int wspsize = nquad0 + nquad1 + nquad2 + max(nqtot,m_ncoeffs)
726  + nquad1*nquad2*nmodes0 + nquad2*nmodes0*(2*nmodes1-nmodes0+1)/2;
727 
728  Array<OneD, NekDouble> gfac0(wspsize);
729  Array<OneD, NekDouble> gfac1(gfac0 + nquad0);
730  Array<OneD, NekDouble> gfac2(gfac1 + nquad1);
731  Array<OneD, NekDouble> tmp0 (gfac2 + nquad2);
732  Array<OneD, NekDouble> wsp(tmp0 + max(nqtot,m_ncoeffs));
733 
734  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
735  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
736  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
737 
738  // set up geometric factor: (1+z0)/2
739  for(i = 0; i < nquad0; ++i)
740  {
741  gfac0[i] = 0.5*(1+z0[i]);
742  }
743 
744  // set up geometric factor: 2/(1-z1)
745  for(i = 0; i < nquad1; ++i)
746  {
747  gfac1[i] = 2.0/(1-z1[i]);
748  }
749 
750  // Set up geometric factor: 2/(1-z2)
751  for(i = 0; i < nquad2; ++i)
752  {
753  gfac2[i] = 2.0/(1-z2[i]);
754  }
755 
756  // Derivative in first direction is always scaled as follows
757  for(i = 0; i < nquad1*nquad2; ++i)
758  {
759  Vmath::Smul(nquad0,gfac1[i%nquad1],&inarray[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
760  }
761  for(i = 0; i < nquad2; ++i)
762  {
763  Vmath::Smul(nquad0*nquad1,gfac2[i],&tmp0[0]+i*nquad0*nquad1,1,&tmp0[0]+i*nquad0*nquad1,1);
764  }
765 
766  MultiplyByQuadratureMetric(tmp0,tmp0);
767 
768  switch(dir)
769  {
770  case 0:
771  {
772  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
773  m_base[1]->GetBdata(),
774  m_base[2]->GetBdata(),
775  tmp0,outarray,wsp,
776  false, true, true);
777  }
778  break;
779  case 1:
780  {
781  Array<OneD, NekDouble> tmp3(m_ncoeffs);
782 
783  for(i = 0; i < nquad1*nquad2; ++i)
784  {
785  Vmath::Vmul(nquad0,&gfac0[0],1,&tmp0[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
786  }
787 
788  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
789  m_base[1]->GetBdata(),
790  m_base[2]->GetBdata(),
791  tmp0,tmp3,wsp,
792  false, true, true);
793 
794  for(i = 0; i < nquad2; ++i)
795  {
796  Vmath::Smul(nquad0*nquad1,gfac2[i],&inarray[0]+i*nquad0*nquad1,1,&tmp0[0]+i*nquad0*nquad1,1);
797  }
798  MultiplyByQuadratureMetric(tmp0,tmp0);
799  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
800  m_base[1]->GetDbdata(),
801  m_base[2]->GetBdata(),
802  tmp0,outarray,wsp,
803  true, false, true);
804  Vmath::Vadd(m_ncoeffs,&tmp3[0],1,&outarray[0],1,&outarray[0],1);
805  }
806  break;
807  case 2:
808  {
809  Array<OneD, NekDouble> tmp3(m_ncoeffs);
810  Array<OneD, NekDouble> tmp4(m_ncoeffs);
811  for(i = 0; i < nquad1; ++i)
812  {
813  gfac1[i] = (1+z1[i])/2;
814  }
815 
816  for(i = 0; i < nquad1*nquad2; ++i)
817  {
818  Vmath::Vmul(nquad0,&gfac0[0],1,&tmp0[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
819  }
820  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
821  m_base[1]->GetBdata(),
822  m_base[2]->GetBdata(),
823  tmp0,tmp3,wsp,
824  false, true, true);
825 
826  for(i = 0; i < nquad2; ++i)
827  {
828  Vmath::Smul(nquad0*nquad1,gfac2[i],&inarray[0]+i*nquad0*nquad1,1,&tmp0[0]+i*nquad0*nquad1,1);
829  }
830  for(i = 0; i < nquad1*nquad2; ++i)
831  {
832  Vmath::Smul(nquad0,gfac1[i%nquad1],&tmp0[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
833  }
834  MultiplyByQuadratureMetric(tmp0,tmp0);
835  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
836  m_base[1]->GetDbdata(),
837  m_base[2]->GetBdata(),
838  tmp0,tmp4,wsp,
839  true, false, true);
840 
841  MultiplyByQuadratureMetric(inarray,tmp0);
842  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
843  m_base[1]->GetBdata(),
844  m_base[2]->GetDbdata(),
845  tmp0,outarray,wsp,
846  true, true, false);
847 
848  Vmath::Vadd(m_ncoeffs,&tmp3[0],1,&outarray[0],1,&outarray[0],1);
849  Vmath::Vadd(m_ncoeffs,&tmp4[0],1,&outarray[0],1,&outarray[0],1);
850  }
851  break;
852  default:
853  {
854  ASSERTL1(false, "input dir is out of range");
855  }
856  break;
857  }
858  }
859 
860 
861  //---------------------------------------
862  // Evaluation functions
863  //---------------------------------------
864 
865 
866  void StdTetExp::v_LocCoordToLocCollapsed(
867  const Array<OneD, const NekDouble>& xi,
868  Array<OneD, NekDouble>& eta)
869  {
870  if( fabs(xi[2]-1.0) < NekConstants::kNekZeroTol)
871  {
872  // Very top point of the tetrahedron
873  eta[0] = -1.0;
874  eta[1] = -1.0;
875  eta[2] = xi[2];
876  }
877  else
878  {
879  if( fabs(xi[1]-1.0) < NekConstants::kNekZeroTol )
880  {
881  // Distant diagonal edge shared by all eta_x
882  // coordinate planes: the xi_y == -xi_z line
883  eta[0] = -1.0;
884  }
885  else if (fabs(xi[1] + xi[2]) < NekConstants::kNekZeroTol)
886  {
887  eta[0] = -1.0;
888  }
889  else
890  {
891  eta[0] = 2.0*(1.0+xi[0])/(-xi[1]-xi[2]) - 1.0;
892  }
893  eta[1] = 2.0*(1.0+xi[1])/(1.0-xi[2]) - 1.0;
894  eta[2] = xi[2];
895  }
896  }
897 
898  void StdTetExp::v_FillMode(
899  const int mode,
900  Array<OneD, NekDouble> &outarray)
901  {
902  Array<OneD, NekDouble> tmp(m_ncoeffs,0.0);
903  tmp[mode] = 1.0;
904  StdTetExp::v_BwdTrans(tmp, outarray);
905  }
906 
907 
908  //---------------------------
909  // Helper functions
910  //---------------------------
911 
912  int StdTetExp::v_GetNverts() const
913  {
914  return 4;
915  }
916 
917  int StdTetExp::v_GetNedges() const
918  {
919  return 6;
920  }
921 
922  int StdTetExp::v_GetNfaces() const
923  {
924  return 4;
925  }
926 
927  LibUtilities::ShapeType StdTetExp::v_DetShapeType() const
928  {
929  return DetShapeType();
930  }
931 
932  int StdTetExp::v_NumBndryCoeffs() const
933  {
934  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
935  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
936  "BasisType is not a boundary interior form");
937  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
938  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
939  "BasisType is not a boundary interior form");
940  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
941  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
942  "BasisType is not a boundary interior form");
943 
944  int P = m_base[0]->GetNumModes();
945  int Q = m_base[1]->GetNumModes();
946  int R = m_base[2]->GetNumModes();
947 
948  return LibUtilities::StdTetData::
949  getNumberOfBndCoefficients(P, Q, R);
950  }
951 
952  int StdTetExp::v_NumDGBndryCoeffs() const
953  {
954  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
955  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
956  "BasisType is not a boundary interior form");
957  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
958  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
959  "BasisType is not a boundary interior form");
960  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
961  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
962  "BasisType is not a boundary interior form");
963 
964  int P = m_base[0]->GetNumModes()-1;
965  int Q = m_base[1]->GetNumModes()-1;
966  int R = m_base[2]->GetNumModes()-1;
967 
968 
969  return (Q+1) + P*(1 + 2*Q - P)/2 // base face
970  + (R+1) + P*(1 + 2*R - P)/2 // front face
971  + 2*(R+1) + Q*(1 + 2*R - Q); // back two faces
972  }
973 
974  int StdTetExp::v_GetEdgeNcoeffs(const int i) const
975  {
976  ASSERTL2((i >= 0) && (i <= 5), "edge id is out of range");
977  int P = m_base[0]->GetNumModes();
978  int Q = m_base[1]->GetNumModes();
979  int R = m_base[2]->GetNumModes();
980 
981  if (i == 0)
982  {
983  return P;
984  }
985  else if (i == 1 || i == 2)
986  {
987  return Q;
988  }
989  else
990  {
991  return R;
992  }
993  }
994 
995  int StdTetExp::v_GetTotalEdgeIntNcoeffs() const
996  {
997  int P = m_base[0]->GetNumModes()-2;
998  int Q = m_base[1]->GetNumModes()-2;
999  int R = m_base[2]->GetNumModes()-2;
1000 
1001  return P+Q+4*R;
1002  }
1003 
1004  int StdTetExp::v_GetFaceNcoeffs(const int i) const
1005  {
1006  ASSERTL2((i >= 0) && (i <= 3), "face id is out of range");
1007  int nFaceCoeffs = 0;
1008  int nummodesA, nummodesB, P, Q;
1009  if (i == 0)
1010  {
1011  nummodesA = GetBasisNumModes(0);
1012  nummodesB = GetBasisNumModes(1);
1013  }
1014  else if ((i == 1) || (i == 2))
1015  {
1016  nummodesA = GetBasisNumModes(0);
1017  nummodesB = GetBasisNumModes(2);
1018  }
1019  else
1020  {
1021  nummodesA = GetBasisNumModes(1);
1022  nummodesB = GetBasisNumModes(2);
1023  }
1024  P = nummodesA - 1;
1025  Q = nummodesB - 1;
1026  nFaceCoeffs = Q+1 + (P*(1 + 2*Q - P))/2;
1027  return nFaceCoeffs;
1028  }
1029 
1030  int StdTetExp::v_GetFaceIntNcoeffs(const int i) const
1031  {
1032  ASSERTL2((i >= 0) && (i <= 3), "face id is out of range");
1033  int Pi = m_base[0]->GetNumModes() - 2;
1034  int Qi = m_base[1]->GetNumModes() - 2;
1035  int Ri = m_base[2]->GetNumModes() - 2;
1036 
1037  if((i == 0))
1038  {
1039  return Pi * (2*Qi - Pi - 1) / 2;
1040  }
1041  else if((i == 1))
1042  {
1043  return Pi * (2*Ri - Pi - 1) / 2;
1044  }
1045  else
1046  {
1047  return Qi * (2*Ri - Qi - 1) / 2;
1048  }
1049  }
1050 
1051  int StdTetExp::v_GetTotalFaceIntNcoeffs() const
1052  {
1053  int Pi = m_base[0]->GetNumModes() - 2;
1054  int Qi = m_base[1]->GetNumModes() - 2;
1055  int Ri = m_base[2]->GetNumModes() - 2;
1056 
1057  return Pi * (2*Qi - Pi - 1) / 2 +
1058  Pi * (2*Ri - Pi - 1) / 2 +
1059  Qi * (2*Ri - Qi - 1);
1060  }
1061 
1062  int StdTetExp::v_GetFaceNumPoints(const int i) const
1063  {
1064  ASSERTL2(i >= 0 && i <= 3, "face id is out of range");
1065 
1066  if (i == 0)
1067  {
1068  return m_base[0]->GetNumPoints()*
1069  m_base[1]->GetNumPoints();
1070  }
1071  else if (i == 1)
1072  {
1073  return m_base[0]->GetNumPoints()*
1074  m_base[2]->GetNumPoints();
1075  }
1076  else
1077  {
1078  return m_base[1]->GetNumPoints()*
1079  m_base[2]->GetNumPoints();
1080  }
1081  }
1082 
1083  LibUtilities::PointsKey StdTetExp::v_GetFacePointsKey(
1084  const int i, const int j) const
1085  {
1086  ASSERTL2(i >= 0 && i <= 3, "face id is out of range");
1087  ASSERTL2(j == 0 || j == 1, "face direction is out of range");
1088 
1089  if (i == 0)
1090  {
1091  return m_base[j]->GetPointsKey();
1092  }
1093  else if (i == 1)
1094  {
1095  return m_base[2*j]->GetPointsKey();
1096  }
1097  else
1098  {
1099  return m_base[j+1]->GetPointsKey();
1100  }
1101  }
1102 
1103  int StdTetExp::v_CalcNumberOfCoefficients(
1104  const std::vector<unsigned int>& nummodes,
1105  int & modes_offset)
1106  {
1107  int nmodes = LibUtilities::StdTetData::getNumberOfCoefficients(
1108  nummodes[modes_offset],
1109  nummodes[modes_offset+1],
1110  nummodes[modes_offset+2]);
1111  modes_offset += 3;
1112 
1113  return nmodes;
1114  }
1115 
1116  const LibUtilities::BasisKey StdTetExp::v_DetFaceBasisKey(
1117  const int i, const int k) const
1118  {
1119  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
1120  ASSERTL2(k == 0 || k == 1, "face direction out of range");
1121 
1122  int dir = k;
1123  switch(i)
1124  {
1125  case 0:
1126  dir = k;
1127  break;
1128  case 1:
1129  dir = 2*k;
1130  break;
1131  case 2:
1132  case 3:
1133  dir = k+1;
1134  break;
1135  }
1136 
1137  return EvaluateTriFaceBasisKey(k,
1138  m_base[dir]->GetBasisType(),
1139  m_base[dir]->GetNumPoints(),
1140  m_base[dir]->GetNumModes());
1141 
1142  // Should not get here.
1143  return LibUtilities::NullBasisKey;
1144  }
1145 
1146  LibUtilities::BasisType StdTetExp::v_GetEdgeBasisType(const int i) const
1147  {
1148  ASSERTL2(i >= 0 && i <= 5, "edge id is out of range");
1149 
1150  if (i == 0)
1151  {
1152  return GetBasisType(0);
1153  }
1154  else if (i == 1 || i == 2)
1155  {
1156  return GetBasisType(1);
1157  }
1158  else
1159  {
1160  return GetBasisType(2);
1161  }
1162  }
1163 
1164  void StdTetExp::v_GetCoords(
1165  Array<OneD, NekDouble> &xi_x,
1166  Array<OneD, NekDouble> &xi_y,
1167  Array<OneD, NekDouble> &xi_z)
1168  {
1169  Array<OneD, const NekDouble> eta_x = m_base[0]->GetZ();
1170  Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
1171  Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
1172  int Qx = GetNumPoints(0);
1173  int Qy = GetNumPoints(1);
1174  int Qz = GetNumPoints(2);
1175 
1176  // Convert collapsed coordinates into cartesian coordinates: eta
1177  // --> xi
1178  for( int k = 0; k < Qz; ++k ) {
1179  for( int j = 0; j < Qy; ++j ) {
1180  for( int i = 0; i < Qx; ++i ) {
1181  int s = i + Qx*(j + Qy*k);
1182  xi_x[s] = (eta_x[i] + 1.0) * (1.0 - eta_y[j]) * (1.0 - eta_z[k]) / 4 - 1.0;
1183  xi_y[s] = (eta_y[j] + 1.0) * (1.0 - eta_z[k]) / 2 - 1.0;
1184  xi_z[s] = eta_z[k];
1185  }
1186  }
1187  }
1188  }
1189 
1190  bool StdTetExp::v_IsBoundaryInteriorExpansion()
1191  {
1192  return (m_base[0]->GetBasisType() == LibUtilities::eModified_A) &&
1193  (m_base[1]->GetBasisType() == LibUtilities::eModified_B) &&
1194  (m_base[2]->GetBasisType() == LibUtilities::eModified_C);
1195  }
1196 
1197 
1198  //--------------------------
1199  // Mappings
1200  //--------------------------
1201 
1202  /**
1203  * Maps Expansion2D modes of a 2D face to the corresponding expansion
1204  * modes.
1205  */
1206  void StdTetExp::v_GetFaceToElementMap(
1207  const int fid,
1208  const Orientation faceOrient,
1209  Array<OneD, unsigned int> &maparray,
1210  Array<OneD, int> &signarray,
1211  int nummodesA,
1212  int nummodesB)
1213  {
1214  int P, Q, i, j, k, idx;
1215 
1216  ASSERTL1(v_IsBoundaryInteriorExpansion(),
1217  "Method only implemented for Modified_A BasisType (x "
1218  "direction), Modified_B BasisType (y direction), and "
1219  "Modified_C BasisType(z direction)");
1220 
1221  int nFaceCoeffs = 0;
1222 
1223  if (nummodesA == -1)
1224  {
1225  switch(fid)
1226  {
1227  case 0:
1228  nummodesA = m_base[0]->GetNumModes();
1229  nummodesB = m_base[1]->GetNumModes();
1230  break;
1231  case 1:
1232  nummodesA = m_base[0]->GetNumModes();
1233  nummodesB = m_base[2]->GetNumModes();
1234  break;
1235  case 2:
1236  case 3:
1237  nummodesA = m_base[1]->GetNumModes();
1238  nummodesB = m_base[2]->GetNumModes();
1239  break;
1240  }
1241  }
1242 
1243  P = nummodesA;
1244  Q = nummodesB;
1245 
1246  nFaceCoeffs = Q + ((P-1)*(1 + 2*(Q-1) - (P-1)))/2;
1247 
1248  // Allocate the map array and sign array; set sign array to ones (+)
1249  if(maparray.num_elements() != nFaceCoeffs)
1250  {
1251  maparray = Array<OneD, unsigned int>(nFaceCoeffs,1);
1252  }
1253 
1254  if(signarray.num_elements() != nFaceCoeffs)
1255  {
1256  signarray = Array<OneD, int>(nFaceCoeffs,1);
1257  }
1258  else
1259  {
1260  fill( signarray.get() , signarray.get()+nFaceCoeffs, 1 );
1261  }
1262 
1263  switch (fid)
1264  {
1265  case 0:
1266  idx = 0;
1267  for (i = 0; i < P; ++i)
1268  {
1269  for (j = 0; j < Q-i; ++j)
1270  {
1271  if ((int)faceOrient == 7 && i > 1)
1272  {
1273  signarray[idx] = (i%2 ? -1 : 1);
1274  }
1275  maparray[idx++] = GetMode(i,j,0);
1276  }
1277  }
1278  break;
1279  case 1:
1280  idx = 0;
1281  for (i = 0; i < P; ++i)
1282  {
1283  for (k = 0; k < Q-i; ++k)
1284  {
1285  if ((int)faceOrient == 7 && i > 1)
1286  {
1287  signarray[idx] = (i%2 ? -1: 1);
1288  }
1289  maparray[idx++] = GetMode(i,0,k);
1290  }
1291  }
1292  break;
1293  case 2:
1294  idx = 0;
1295  for (j = 0; j < P-1; ++j)
1296  {
1297  for (k = 0; k < Q-1-j; ++k)
1298  {
1299  if ((int)faceOrient == 7 && j > 1)
1300  {
1301  signarray[idx] = ((j+1)%2 ? -1: 1);
1302  }
1303  maparray[idx++] = GetMode(1,j,k);
1304  // Incorporate modes from zeroth plane where needed.
1305  if (j == 0 && k == 0) {
1306  maparray[idx++] = GetMode(0,0,1);
1307  }
1308  if (j == 0 && k == Q-2) {
1309  for (int r = 0; r < Q-1; ++r) {
1310  maparray[idx++] = GetMode(0,1,r);
1311  }
1312  }
1313  }
1314  }
1315  break;
1316  case 3:
1317  idx = 0;
1318  for (j = 0; j < P; ++j)
1319  {
1320  for (k = 0; k < Q-j; ++k)
1321  {
1322  if ((int)faceOrient == 7 && j > 1)
1323  {
1324  signarray[idx] = (j%2 ? -1: 1);
1325  }
1326  maparray[idx++] = GetMode(0,j,k);
1327  }
1328  }
1329  break;
1330  default:
1331  ASSERTL0(false, "Element map not available.");
1332  }
1333 
1334  if ((int)faceOrient == 7)
1335  {
1336  swap(maparray[0], maparray[Q]);
1337 
1338  for (i = 1; i < Q-1; ++i)
1339  {
1340  swap(maparray[i+1], maparray[Q+i]);
1341  }
1342  }
1343  }
1344 
1345  int StdTetExp::v_GetVertexMap(const int localVertexId, bool useCoeffPacking)
1346  {
1347  ASSERTL0((GetEdgeBasisType(localVertexId)==LibUtilities::eModified_A)||
1348  (GetEdgeBasisType(localVertexId)==LibUtilities::eModified_B)||
1349  (GetEdgeBasisType(localVertexId)==LibUtilities::eModified_C),
1350  "Mapping not defined for this type of basis");
1351 
1352  int localDOF = 0;
1353  if(useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1354  {
1355  switch(localVertexId)
1356  {
1357  case 0:
1358  {
1359  localDOF = GetMode(0,0,0);
1360  break;
1361  }
1362  case 1:
1363  {
1364  localDOF = GetMode(0,0,1);
1365  break;
1366  }
1367  case 2:
1368  {
1369  localDOF = GetMode(0,1,0);
1370  break;
1371  }
1372  case 3:
1373  {
1374  localDOF = GetMode(1,0,0);
1375  break;
1376  }
1377  default:
1378  {
1379  ASSERTL0(false,"Vertex ID must be between 0 and 3");
1380  break;
1381  }
1382  }
1383  }
1384  else
1385  {
1386  switch(localVertexId)
1387  {
1388  case 0:
1389  {
1390  localDOF = GetMode(0,0,0);
1391  break;
1392  }
1393  case 1:
1394  {
1395  localDOF = GetMode(1,0,0);
1396  break;
1397  }
1398  case 2:
1399  {
1400  localDOF = GetMode(0,1,0);
1401  break;
1402  }
1403  case 3:
1404  {
1405  localDOF = GetMode(0,0,1);
1406  break;
1407  }
1408  default:
1409  {
1410  ASSERTL0(false,"Vertex ID must be between 0 and 3");
1411  break;
1412  }
1413  }
1414 
1415  }
1416 
1417  return localDOF;
1418  }
1419 
1420  /**
1421  * Maps interior modes of an edge to the elemental modes.
1422  */
1423  void StdTetExp::v_GetEdgeInteriorMap(
1424  const int eid,
1425  const Orientation edgeOrient,
1426  Array<OneD, unsigned int> &maparray,
1427  Array<OneD, int> &signarray)
1428  {
1429  int i;
1430  const int P = m_base[0]->GetNumModes();
1431  const int Q = m_base[1]->GetNumModes();
1432  const int R = m_base[2]->GetNumModes();
1433 
1434  const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid)-2;
1435 
1436  if(maparray.num_elements() != nEdgeIntCoeffs)
1437  {
1438  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1439  }
1440  else
1441  {
1442  fill( maparray.get(), maparray.get() + nEdgeIntCoeffs, 0);
1443  }
1444 
1445  if(signarray.num_elements() != nEdgeIntCoeffs)
1446  {
1447  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1448  }
1449  else
1450  {
1451  fill( signarray.get() , signarray.get()+nEdgeIntCoeffs, 1 );
1452  }
1453 
1454  switch (eid)
1455  {
1456  case 0:
1457  for (i = 0; i < P-2; ++i)
1458  {
1459  maparray[i] = GetMode(i+2, 0, 0);
1460  }
1461  if(edgeOrient==eBackwards)
1462  {
1463  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1464  {
1465  signarray[i] = -1;
1466  }
1467  }
1468  break;
1469  case 1:
1470  for (i = 0; i < Q-2; ++i)
1471  {
1472  maparray[i] = GetMode(1, i+1, 0);
1473  }
1474  if(edgeOrient==eBackwards)
1475  {
1476  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1477  {
1478  signarray[i] = -1;
1479  }
1480  }
1481  break;
1482  case 2:
1483  for (i = 0; i < Q-2; ++i)
1484  {
1485  maparray[i] = GetMode(0, i+2, 0);
1486  }
1487  if(edgeOrient==eBackwards)
1488  {
1489  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1490  {
1491  signarray[i] = -1;
1492  }
1493  }
1494  break;
1495  case 3:
1496  for (i = 0; i < R-2; ++i)
1497  {
1498  maparray[i] = GetMode(0, 0, i+2);
1499  }
1500  if(edgeOrient==eBackwards)
1501  {
1502  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1503  {
1504  signarray[i] = -1;
1505  }
1506  }
1507  break;
1508  case 4:
1509  for (i = 0; i < R - 2; ++i)
1510  {
1511  maparray[i] = GetMode(1, 0, i+1);
1512  }
1513  if(edgeOrient==eBackwards)
1514  {
1515  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1516  {
1517  signarray[i] = -1;
1518  }
1519  }
1520  break;
1521  case 5:
1522  for (i = 0; i < R-2; ++i)
1523  {
1524  maparray[i] = GetMode(0, 1, i+1);
1525  }
1526  if(edgeOrient==eBackwards)
1527  {
1528  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1529  {
1530  signarray[i] = -1;
1531  }
1532  }
1533  break;
1534  default:
1535  ASSERTL0(false, "Edge not defined.");
1536  break;
1537  }
1538  }
1539 
1540  void StdTetExp::v_GetFaceInteriorMap(
1541  const int fid,
1542  const Orientation faceOrient,
1543  Array<OneD, unsigned int> &maparray,
1544  Array<OneD, int> &signarray)
1545  {
1546  int i, j, idx, k;
1547  const int P = m_base[0]->GetNumModes();
1548  const int Q = m_base[1]->GetNumModes();
1549  const int R = m_base[2]->GetNumModes();
1550 
1551  const int nFaceIntCoeffs = v_GetFaceIntNcoeffs(fid);
1552 
1553  if(maparray.num_elements() != nFaceIntCoeffs)
1554  {
1555  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1556  }
1557 
1558  if(signarray.num_elements() != nFaceIntCoeffs)
1559  {
1560  signarray = Array<OneD, int>(nFaceIntCoeffs,1);
1561  }
1562  else
1563  {
1564  fill( signarray.get() , signarray.get()+nFaceIntCoeffs, 1 );
1565  }
1566 
1567  switch (fid)
1568  {
1569  case 0:
1570  idx = 0;
1571  for (i = 2; i < P-1; ++i)
1572  {
1573  for (j = 1; j < Q-i; ++j)
1574  {
1575  if ((int)faceOrient == 7)
1576  {
1577  signarray[idx] = (i%2 ? -1 : 1);
1578  }
1579  maparray[idx++] = GetMode(i,j,0);
1580  }
1581  }
1582  break;
1583  case 1:
1584  idx = 0;
1585  for (i = 2; i < P; ++i)
1586  {
1587  for (k = 1; k < R-i; ++k)
1588  {
1589  if ((int)faceOrient == 7)
1590  {
1591  signarray[idx] = (i%2 ? -1: 1);
1592  }
1593  maparray[idx++] = GetMode(i,0,k);
1594  }
1595  }
1596  break;
1597  case 2:
1598  idx = 0;
1599  for (j = 1; j < Q-2; ++j)
1600  {
1601  for (k = 1; k < R-1-j; ++k)
1602  {
1603  if ((int)faceOrient == 7)
1604  {
1605  signarray[idx] = ((j+1)%2 ? -1: 1);
1606  }
1607  maparray[idx++] = GetMode(1,j,k);
1608  }
1609  }
1610  break;
1611  case 3:
1612  idx = 0;
1613  for (j = 2; j < Q-1; ++j)
1614  {
1615  for (k = 1; k < R-j; ++k)
1616  {
1617  if ((int)faceOrient == 7)
1618  {
1619  signarray[idx] = (j%2 ? -1: 1);
1620  }
1621  maparray[idx++] = GetMode(0,j,k);
1622  }
1623  }
1624  break;
1625  default:
1626  ASSERTL0(false, "Face interior map not available.");
1627  break;
1628  }
1629  }
1630 
1631  /**
1632  * List of all interior modes in the expansion.
1633  */
1634  void StdTetExp::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
1635  {
1636  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1637  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1638  "BasisType is not a boundary interior form");
1639  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
1640  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1641  "BasisType is not a boundary interior form");
1642  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1643  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1644  "BasisType is not a boundary interior form");
1645 
1646  int P = m_base[0]->GetNumModes();
1647  int Q = m_base[1]->GetNumModes();
1648  int R = m_base[2]->GetNumModes();
1649 
1650  int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1651 
1652  if(outarray.num_elements() != nIntCoeffs)
1653  {
1654  outarray = Array<OneD, unsigned int>(nIntCoeffs);
1655  }
1656 
1657  int idx = 0;
1658  for (int i = 2; i < P-2; ++i)
1659  {
1660  for (int j = 1; j < Q-i-1; ++j)
1661  {
1662  for (int k = 1; k < R-i-j; ++k)
1663  {
1664  outarray[idx++] = GetMode(i,j,k);
1665  }
1666  }
1667  }
1668  }
1669 
1670  /**
1671  * List of all boundary modes in the the expansion.
1672  */
1673  void StdTetExp::v_GetBoundaryMap(Array<OneD, unsigned int>& outarray)
1674  {
1675  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1676  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1677  "BasisType is not a boundary interior form");
1678  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
1679  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1680  "BasisType is not a boundary interior form");
1681  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1682  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1683  "BasisType is not a boundary interior form");
1684 
1685  int P = m_base[0]->GetNumModes();
1686  int Q = m_base[1]->GetNumModes();
1687  int R = m_base[2]->GetNumModes();
1688 
1689  int i,j,k;
1690  int idx = 0;
1691 
1692  int nBnd = NumBndryCoeffs();
1693 
1694  if (outarray.num_elements() != nBnd)
1695  {
1696  outarray = Array<OneD, unsigned int>(nBnd);
1697  }
1698 
1699  for (i = 0; i < P; ++i)
1700  {
1701  // First two Q-R planes are entirely boundary modes
1702  if (i < 2)
1703  {
1704  for (j = 0; j < Q-i; j++)
1705  {
1706  for (k = 0; k < R-i-j; ++k)
1707  {
1708  outarray[idx++] = GetMode(i,j,k);
1709  }
1710  }
1711  }
1712  // Remaining Q-R planes contain boundary modes on bottom and
1713  // left edge.
1714  else
1715  {
1716  for (k = 0; k < R-i; ++k)
1717  {
1718  outarray[idx++] = GetMode(i,0,k);
1719  }
1720  for (j = 1; j < Q-i; ++j)
1721  {
1722  outarray[idx++] = GetMode(i,j,0);
1723  }
1724  }
1725  }
1726  }
1727 
1728 
1729  //---------------------------------------
1730  // Wrapper functions
1731  //---------------------------------------
1732  DNekMatSharedPtr StdTetExp::v_GenMatrix(const StdMatrixKey &mkey)
1733  {
1734 
1735  MatrixType mtype = mkey.GetMatrixType();
1736 
1737  DNekMatSharedPtr Mat;
1738 
1739  switch(mtype)
1740  {
1741  case ePhysInterpToEquiSpaced:
1742  {
1743  int nq0 = m_base[0]->GetNumPoints();
1744  int nq1 = m_base[1]->GetNumPoints();
1745  int nq2 = m_base[2]->GetNumPoints();
1746  int nq = max(nq0,max(nq1,nq2));
1747  int neq = LibUtilities::StdTetData::
1748  getNumberOfCoefficients(nq,nq,nq);
1749  Array<OneD, Array<OneD, NekDouble> > coords(neq);
1750  Array<OneD, NekDouble> coll(3);
1751  Array<OneD, DNekMatSharedPtr> I(3);
1752  Array<OneD, NekDouble> tmp(nq0);
1753 
1754  Mat = MemoryManager<DNekMat>::
1755  AllocateSharedPtr(neq, nq0 * nq1 * nq2);
1756  int cnt = 0;
1757 
1758  for(int i = 0; i < nq; ++i)
1759  {
1760  for(int j = 0; j < nq-i; ++j)
1761  {
1762  for(int k = 0; k < nq-i-j; ++k,++cnt)
1763  {
1764  coords[cnt] = Array<OneD, NekDouble>(3);
1765  coords[cnt][0] = -1.0 + 2*k/(NekDouble)(nq-1);
1766  coords[cnt][1] = -1.0 + 2*j/(NekDouble)(nq-1);
1767  coords[cnt][2] = -1.0 + 2*i/(NekDouble)(nq-1);
1768  }
1769  }
1770  }
1771 
1772  for(int i = 0; i < neq; ++i)
1773  {
1774  LocCoordToLocCollapsed(coords[i],coll);
1775 
1776  I[0] = m_base[0]->GetI(coll);
1777  I[1] = m_base[1]->GetI(coll+1);
1778  I[2] = m_base[2]->GetI(coll+2);
1779 
1780  // interpolate first coordinate direction
1781  NekDouble fac;
1782  for( int k = 0; k < nq2; ++k)
1783  {
1784  for (int j = 0; j < nq1; ++j)
1785  {
1786 
1787  fac = (I[1]->GetPtr())[j]*(I[2]->GetPtr())[k];
1788  Vmath::Smul(nq0,fac,I[0]->GetPtr(),1,tmp,1);
1789 
1790  Vmath::Vcopy(nq0, &tmp[0], 1,
1791  Mat->GetRawPtr()+
1792  k*nq0*nq1*neq+
1793  j*nq0*neq+i,neq);
1794  }
1795  }
1796  }
1797  }
1798  break;
1799  default:
1800  {
1801  Mat = StdExpansion::CreateGeneralMatrix(mkey);
1802  }
1803  break;
1804  }
1805 
1806  return Mat;
1807  }
1808 
1809  DNekMatSharedPtr StdTetExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
1810  {
1811  return v_GenMatrix(mkey);
1812  }
1813 
1814 
1815  //---------------------------------------
1816  // Private helper functions
1817  //---------------------------------------
1818 
1819  /**
1820  * @brief Compute the mode number in the expansion for a particular
1821  * tensorial combination.
1822  *
1823  * Modes are numbered with the r index travelling fastest, followed by
1824  * q and then p, and each q-r plane is of size
1825  * (Q+1)*(Q+2)/2+max(0,R-Q-p)*Q. For example, when P=2, Q=3 and R=4
1826  * the indexing inside each q-r plane (with r increasing upwards and q
1827  * to the right) is:
1828  *
1829  * p = 0: p = 1: p = 2:
1830  * ----------------------------------
1831  * 4
1832  * 3 8 17
1833  * 2 7 11 16 20 26
1834  * 1 6 10 13 15 19 22 25 28
1835  * 0 5 9 12 14 18 21 23 24 27 29
1836  *
1837  * Note that in this element, we must have that \f$ P \leq Q \leq
1838  * R\f$.
1839  */
1840  int StdTetExp::GetMode(const int I, const int J, const int K)
1841  {
1842  const int Q = m_base[1]->GetNumModes();
1843  const int R = m_base[2]->GetNumModes();
1844 
1845  int i,j,q_hat,k_hat;
1846  int cnt = 0;
1847 
1848  // Traverse to q-r plane number I
1849  for (i = 0; i < I; ++i)
1850  {
1851  // Size of triangle part
1852  q_hat = min(Q,R-i);
1853  // Size of rectangle part
1854  k_hat = max(R-Q-i,0);
1855  cnt += q_hat*(q_hat+1)/2 + k_hat*Q;
1856  }
1857 
1858  // Traverse to q column J
1859  q_hat = R-I;
1860  for (j = 0; j < J; ++j)
1861  {
1862  cnt += q_hat;
1863  q_hat--;
1864  }
1865 
1866  // Traverse up stacks to K
1867  cnt += K;
1868 
1869  return cnt;
1870  }
1871 
1872  void StdTetExp::v_MultiplyByStdQuadratureMetric(
1873  const Array<OneD, const NekDouble>& inarray,
1874  Array<OneD, NekDouble>& outarray)
1875  {
1876  int i, j;
1877 
1878  int nquad0 = m_base[0]->GetNumPoints();
1879  int nquad1 = m_base[1]->GetNumPoints();
1880  int nquad2 = m_base[2]->GetNumPoints();
1881 
1882  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1883  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1884  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
1885 
1886  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
1887  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1888 
1889  // multiply by integration constants
1890  for(i = 0; i < nquad1*nquad2; ++i)
1891  {
1892  Vmath::Vmul(nquad0,(NekDouble*)&inarray[0]+i*nquad0,1,
1893  w0.get(),1, &outarray[0]+i*nquad0,1);
1894  }
1895 
1896  switch(m_base[1]->GetPointsType())
1897  {
1898  // (1,0) Jacobi Inner product.
1899  case LibUtilities::eGaussRadauMAlpha1Beta0:
1900  for(j = 0; j < nquad2; ++j)
1901  {
1902  for(i = 0; i < nquad1; ++i)
1903  {
1904  Blas::Dscal(nquad0,0.5*w1[i], &outarray[0]+i*nquad0+
1905  j*nquad0*nquad1,1);
1906  }
1907  }
1908  break;
1909 
1910  default:
1911  for(j = 0; j < nquad2; ++j)
1912  {
1913  for(i = 0; i < nquad1; ++i)
1914  {
1915  Blas::Dscal(nquad0,
1916  0.5*(1-z1[i])*w1[i],
1917  &outarray[0]+i*nquad0 + j*nquad0*nquad1,
1918  1 );
1919  }
1920  }
1921  break;
1922  }
1923 
1924  switch(m_base[2]->GetPointsType())
1925  {
1926  // (2,0) Jacobi inner product.
1927  case LibUtilities::eGaussRadauMAlpha2Beta0:
1928  for(i = 0; i < nquad2; ++i)
1929  {
1930  Blas::Dscal(nquad0*nquad1, 0.25*w2[i],
1931  &outarray[0]+i*nquad0*nquad1, 1);
1932  }
1933  break;
1934 
1935  default:
1936  for(i = 0; i < nquad2; ++i)
1937  {
1938  Blas::Dscal(nquad0*nquad1,0.25*(1-z2[i])*(1-z2[i])*w2[i],
1939  &outarray[0]+i*nquad0*nquad1,1);
1940  }
1941  break;
1942  }
1943  }
1944 
1945  void StdTetExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
1946  const StdMatrixKey &mkey)
1947  {
1948  //To do : 1) add a test to ensure 0 \leq SvvCutoff \leq 1.
1949  // 2) check if the transfer function needs an analytical
1950  // Fourier transform.
1951  // 3) if it doesn't : find a transfer function that renders
1952  // the if( cutoff_a ...) useless to reduce computational
1953  // cost.
1954  // 4) add SVVDiffCoef to both models!!
1955 
1956  int qa = m_base[0]->GetNumPoints();
1957  int qb = m_base[1]->GetNumPoints();
1958  int qc = m_base[2]->GetNumPoints();
1959  int nmodes_a = m_base[0]->GetNumModes();
1960  int nmodes_b = m_base[1]->GetNumModes();
1961  int nmodes_c = m_base[2]->GetNumModes();
1962 
1963  // Declare orthogonal basis.
1964  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
1965  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
1966  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
1967 
1968  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
1969  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B,nmodes_b,pb);
1970  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_C,nmodes_c,pc);
1971 
1972  StdTetExp OrthoExp(Ba,Bb,Bc);
1973 
1974 
1975  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1976  int i,j,k,cnt = 0;
1977 
1978  //SVV filter paramaters (how much added diffusion relative to physical one
1979  // and fraction of modes from which you start applying this added diffusion)
1980  //
1981  NekDouble SvvDiffCoeff = mkey.GetConstFactor(StdRegions::eFactorSVVDiffCoeff);
1982  NekDouble SVVCutOff = mkey.GetConstFactor(StdRegions::eFactorSVVCutoffRatio);
1983 
1984 
1985  //Defining the cut of mode
1986  int cutoff_a = (int) (SVVCutOff*nmodes_a);
1987  int cutoff_b = (int) (SVVCutOff*nmodes_b);
1988  int cutoff_c = (int) (SVVCutOff*nmodes_c);
1989  int nmodes = min(min(nmodes_a,nmodes_b),nmodes_c);
1990  NekDouble cutoff = min(min(cutoff_a,cutoff_b),cutoff_c);
1991  NekDouble epsilon = 1;
1992 
1993  // project onto physical space.
1994  OrthoExp.FwdTrans(array,orthocoeffs);
1995 
1996  //------"New" Version August 22nd '13--------------------
1997  for(i = 0; i < nmodes_a; ++i)
1998  {
1999  for(j = 0; j < nmodes_b-i; ++j)
2000  {
2001  for(k = 0; k < nmodes_c-i-j; ++k)
2002  {
2003  if(i + j + k >= cutoff)
2004  {
2005  orthocoeffs[cnt] *= ((SvvDiffCoeff)*exp(-(i+j+k-nmodes)*(i+j+k-nmodes)/((NekDouble)((i+j+k-cutoff+epsilon)*(i+j+k-cutoff+epsilon)))));
2006  }
2007  else
2008  {
2009  orthocoeffs[cnt] *= 0.0;
2010  }
2011  cnt++;
2012  }
2013  }
2014  }
2015  // backward transform to physical space
2016  OrthoExp.BwdTrans(orthocoeffs,array);
2017  }
2018 
2019 
2020  void StdTetExp::v_ReduceOrderCoeffs(
2021  int numMin,
2022  const Array<OneD, const NekDouble> &inarray,
2023  Array<OneD, NekDouble> &outarray)
2024  {
2025  int nquad0 = m_base[0]->GetNumPoints();
2026  int nquad1 = m_base[1]->GetNumPoints();
2027  int nquad2 = m_base[2]->GetNumPoints();
2028  int nqtot = nquad0 * nquad1 * nquad2;
2029  int nmodes0 = m_base[0]->GetNumModes();
2030  int nmodes1 = m_base[1]->GetNumModes();
2031  int nmodes2 = m_base[2]->GetNumModes();
2032  int numMax = nmodes0;
2033 
2034  Array<OneD, NekDouble> coeff (m_ncoeffs);
2035  Array<OneD, NekDouble> coeff_tmp1(m_ncoeffs, 0.0);
2036  Array<OneD, NekDouble> coeff_tmp2(m_ncoeffs, 0.0);
2037  Array<OneD, NekDouble> phys_tmp (nqtot, 0.0);
2038  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
2039 
2040  Vmath::Vcopy(m_ncoeffs,inarray,1,coeff_tmp2,1);
2041 
2042  const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
2043  const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
2044  const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
2045 
2046  LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A,
2047  nmodes0, Pkey0);
2048  LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_B,
2049  nmodes1, Pkey1);
2050  LibUtilities::BasisKey bortho2(LibUtilities::eOrtho_C,
2051  nmodes2, Pkey2);
2052 
2053  Vmath::Zero(m_ncoeffs, coeff_tmp2, 1);
2054 
2055  StdRegions::StdTetExpSharedPtr OrthoTetExp;
2056  OrthoTetExp = MemoryManager<StdRegions::StdTetExp>
2057  ::AllocateSharedPtr(bortho0, bortho1, bortho2);
2058 
2059  BwdTrans(inarray,phys_tmp);
2060  OrthoTetExp->FwdTrans(phys_tmp, coeff);
2061 
2062  Vmath::Zero(m_ncoeffs,outarray,1);
2063 
2064  // filtering
2065  int cnt = 0;
2066  for (int u = 0; u < numMin; ++u)
2067  {
2068  for (int i = 0; i < numMin-u; ++i)
2069  {
2070  Vmath::Vcopy(numMin - u - i, tmp = coeff + cnt, 1,
2071  tmp2 = coeff_tmp1 + cnt, 1);
2072  cnt += numMax - u - i;
2073  }
2074  for (int i = numMin; i < numMax-u; ++i)
2075  {
2076  cnt += numMax - u - i;
2077  }
2078  }
2079 
2080  OrthoTetExp->BwdTrans(coeff_tmp1,phys_tmp);
2081  FwdTrans(phys_tmp, outarray);
2082  }
2083 
2084 
2085  void StdTetExp::v_GetSimplexEquiSpacedConnectivity(
2086  Array<OneD, int> &conn,
2087  bool standard)
2088  {
2089  int np0 = m_base[0]->GetNumPoints();
2090  int np1 = m_base[1]->GetNumPoints();
2091  int np2 = m_base[2]->GetNumPoints();
2092  int np = max(np0,max(np1,np2));
2093 
2094 
2095  conn = Array<OneD, int>(4*(np-1)*(np-1)*(np-1));
2096 
2097  int row = 0;
2098  int rowp1 = 0;
2099  int plane = 0;
2100  int row1 = 0;
2101  int row1p1 = 0;
2102  int planep1= 0;
2103  int cnt = 0;
2104  for(int i = 0; i < np-1; ++i)
2105  {
2106  planep1 += (np-i)*(np-i+1)/2;
2107  row = 0; // current plane row offset
2108  rowp1 = 0; // current plane row plus one offset
2109  row1 = 0; // next plane row offset
2110  row1p1 = 0; // nex plane row plus one offset
2111  for(int j = 0; j < np-i-1; ++j)
2112  {
2113  rowp1 += np-i-j;
2114  row1p1 += np-i-j-1;
2115  for(int k = 0; k < np-i-j-2; ++k)
2116  {
2117  conn[cnt++] = plane + row +k+1;
2118  conn[cnt++] = plane + row +k;
2119  conn[cnt++] = plane + rowp1 +k;
2120  conn[cnt++] = planep1 + row1 +k;
2121 
2122  conn[cnt++] = plane + row +k+1;
2123  conn[cnt++] = plane + rowp1 +k+1;
2124  conn[cnt++] = planep1 + row1 +k+1;
2125  conn[cnt++] = planep1 + row1 +k;
2126 
2127  conn[cnt++] = plane + rowp1 +k+1;
2128  conn[cnt++] = plane + row +k+1;
2129  conn[cnt++] = plane + rowp1 +k;
2130  conn[cnt++] = planep1 + row1 +k;
2131 
2132  conn[cnt++] = planep1 + row1 +k;
2133  conn[cnt++] = planep1 + row1p1+k;
2134  conn[cnt++] = plane + rowp1 +k;
2135  conn[cnt++] = plane + rowp1 +k+1;
2136 
2137  conn[cnt++] = planep1 + row1 +k;
2138  conn[cnt++] = planep1 + row1p1+k;
2139  conn[cnt++] = planep1 + row1 +k+1;
2140  conn[cnt++] = plane + rowp1 +k+1;
2141 
2142  if(k < np-i-j-3)
2143  {
2144  conn[cnt++] = plane + rowp1 +k+1;
2145  conn[cnt++] = planep1 + row1p1 +k+1;
2146  conn[cnt++] = planep1 + row1 +k+1;
2147  conn[cnt++] = planep1 + row1p1 +k;
2148  }
2149  }
2150 
2151  conn[cnt++] = plane + row + np-i-j-1;
2152  conn[cnt++] = plane + row + np-i-j-2;
2153  conn[cnt++] = plane + rowp1 + np-i-j-2;
2154  conn[cnt++] = planep1 + row1 + np-i-j-2;
2155 
2156  row += np-i-j;
2157  row1 += np-i-j-1;
2158  }
2159  plane += (np-i)*(np-i+1)/2;
2160  }
2161  }
2162 
2163  }//end namespace
2164 }//end namespace
Nektar::StdRegions::StdTetExp::GetMode
int GetMode(const int i, const int j, const int k)
Compute the mode number in the expansion for a particular tensorial combination.
Definition: StdTetExp.cpp:1840
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::StdRegions::StdTetExp::v_NumBndryCoeffs
virtual int v_NumBndryCoeffs() const
Definition: StdTetExp.cpp:932
Nektar::StdRegions::StdTetExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Definition: StdTetExp.cpp:1732
Nektar::StdRegions::StdTetExp::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: StdTetExp.cpp:251
Nektar::LibUtilities::eModified_C
Principle Modified Functions .
Definition: BasisType.h:51
Nektar::StdRegions::StdMatrixKey::GetMatrixType
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
Nektar::StdRegions::StdTetExp::v_GetSimplexEquiSpacedConnectivity
virtual void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true)
Definition: StdTetExp.cpp:2085
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::StdRegions::StdTetExp::v_GetEdgeBasisType
virtual LibUtilities::BasisType v_GetEdgeBasisType(const int i) const
Definition: StdTetExp.cpp:1146
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:784
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::StdExpansion::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:902
Nektar::StdRegions::StdTetExp::v_GetEdgeNcoeffs
virtual int v_GetEdgeNcoeffs(const int i) const
Definition: StdTetExp.cpp:974
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::EvaluateTriFaceBasisKey
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Definition: StdExpansion3D.cpp:401
Nektar::LibUtilities::eModified_A
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::StdRegions::StdTetExp::v_LocCoordToLocCollapsed
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Definition: StdTetExp.cpp:866
Nektar::StdRegions::StdTetExp::v_GetCoords
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_x, Array< OneD, NekDouble > &coords_y, Array< OneD, NekDouble > &coords_z)
Definition: StdTetExp.cpp:1164
Nektar::StdRegions::StdTetExp::v_IProductWRTBase_SumFac
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: StdTetExp.cpp:545
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:1133
Nektar::StdRegions::StdTetExp::v_PhysDeriv
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_dx, Array< OneD, NekDouble > &out_dy, Array< OneD, NekDouble > &out_dz)
Calculate the derivative of the physical points.
Definition: StdTetExp.cpp:106
Nektar::StdRegions::StdTetExp::v_GetNverts
virtual int v_GetNverts() const
Definition: StdTetExp.cpp:912
Nektar::StdRegions::StdTetExp::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: StdTetExp.cpp:356
Nektar::LibUtilities::StdTetData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:186
Nektar::StdRegions::StdTetExp::v_CalcNumberOfCoefficients
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdTetExp.cpp:1103
Nektar::DNekMatSharedPtr
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
Nektar::StdRegions::StdTetExp::v_GetFaceToElementMap
virtual void v_GetFaceToElementMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int nummodesA=-1, int nummodesB=-1)
Definition: StdTetExp.cpp:1206
Nektar::StdRegions::StdTetExp::v_NumDGBndryCoeffs
virtual int v_NumDGBndryCoeffs() const
Definition: StdTetExp.cpp:952
Nektar::StdRegions::StdTetExp::v_GetEdgeInteriorMap
virtual void v_GetEdgeInteriorMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
Definition: StdTetExp.cpp:1423
Nektar::StdRegions::StdTetExp::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: StdTetExp.cpp:581
Nektar::StdRegions::StdTetExp::v_FwdTrans
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:452
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:684
Nektar::StdRegions::StdTetExp::v_GetFaceIntNcoeffs
virtual int v_GetFaceIntNcoeffs(const int i) const
Definition: StdTetExp.cpp:1030
Nektar::StdRegions::StdTetExp::DetShapeType
LibUtilities::ShapeType DetShapeType() const
Definition: StdTetExp.h:70
Nektar::LibUtilities::eGaussRadauMAlpha1Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
Nektar::StdRegions::StdTetExp::v_DetShapeType
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdTetExp.cpp:927
Nektar::StdRegions::StdTetExp::v_GetTotalFaceIntNcoeffs
virtual int v_GetTotalFaceIntNcoeffs() const
Definition: StdTetExp.cpp:1051
Nektar::LibUtilities::eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:47
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::StdTetExp::v_FillMode
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:898
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:397
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::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
Nektar::StdRegions::StdTetExp::v_GetFacePointsKey
virtual LibUtilities::PointsKey v_GetFacePointsKey(const int i, const int j) const
Definition: StdTetExp.cpp:1083
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
Nektar::StdRegions::StdTetExp::v_GetTotalEdgeIntNcoeffs
virtual int v_GetTotalEdgeIntNcoeffs() const
Definition: StdTetExp.cpp:995
Nektar::StdRegions::StdTetExp::v_BwdTrans_SumFac
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:323
Nektar::StdRegions::StdTetExp::v_IProductWRTDerivBase_SumFac
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:713
Nektar::StdRegions::StdTetExp::v_GetBoundaryMap
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdTetExp.cpp:1673
Nektar::LibUtilities::eModified_B
Principle Modified Functions .
Definition: BasisType.h:50
Nektar::StdRegions::StdTetExp::v_IProductWRTDerivBase
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:667
Nektar::StdRegions::StdTetExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdTetExp.cpp:1809
Nektar::StdRegions::StdTetExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdTetExp.cpp:1634
Nektar::LibUtilities::eOrtho_C
Principle Orthogonal Functions .
Definition: BasisType.h:48
Nektar::StdRegions::StdTetExp::v_MultiplyByStdQuadratureMetric
virtual void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:1872
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::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:301
Nektar::StdRegions::StdTetExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdTetExp.cpp:1945
Vmath::Ddot
T Ddot(int n, const Array< OneD, const T > &w, const int incw, const Array< OneD, const T > &x, const int incx, const Array< OneD, const int > &y, const int incy)
Definition: VmathArray.hpp:434
Nektar::StdRegions::StdTetExp::v_GetFaceInteriorMap
virtual void v_GetFaceInteriorMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
Definition: StdTetExp.cpp:1540
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::StdTetExp::v_IsBoundaryInteriorExpansion
virtual bool v_IsBoundaryInteriorExpansion()
Definition: StdTetExp.cpp:1190
Nektar::StdRegions::StdTetExp::v_GetNfaces
virtual int v_GetNfaces() const
Definition: StdTetExp.cpp:922
Nektar::StdRegions::StdTetExp::v_IProductWRTDerivBase_MatOp
virtual void v_IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:676
Nektar::StdRegions::StdTetExp::v_IProductWRTBase_MatOp
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:526
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:213
Nektar::StdRegions::StdTetExp::v_GetFaceNcoeffs
virtual int v_GetFaceNcoeffs(const int i) const
Definition: StdTetExp.cpp:1004
Nektar::LibUtilities::eGaussRadauMAlpha2Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
Nektar::StdRegions::StdTetExp::v_DetFaceBasisKey
virtual const LibUtilities::BasisKey v_DetFaceBasisKey(const int i, const int k) const
Definition: StdTetExp.cpp:1116
Nektar::StdRegions::StdTetExp::v_GetVertexMap
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdTetExp.cpp:1345
Nektar::StdRegions::StdTetExp::v_BwdTrans
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:293
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:509
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::StdTetExp::v_GetFaceNumPoints
virtual int v_GetFaceNumPoints(const int i) const
Definition: StdTetExp.cpp:1062
Nektar::StdRegions::StdTetExpSharedPtr
boost::shared_ptr< StdTetExp > StdTetExpSharedPtr
Definition: StdTetExp.h:268
Nektar::LibUtilities::BasisKey::GetNumModes
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
Nektar::LibUtilities::eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:53
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Nektar::StdRegions::StdExpansion::GetPointsType
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:216
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1356
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
Nektar::StdRegions::StdTetExp::v_GetNedges
virtual int v_GetNedges() const
Definition: StdTetExp.cpp:917
Nektar::StdRegions::StdTetExp::v_ReduceOrderCoeffs
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:2020
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::StdRegions::StdTetExp::v_IProductWRTBase
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:503
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
Nektar::LibUtilities::StdTetData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:209
StdTetExp.h
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:1823
Generated on Tue Dec 8 2015 23:45:30 for Nektar++ by   doxygen 1.8.9.1