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_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*(2*order1-order0+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*(2*order1-order0+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 P,
1212  int Q)
1213  {
1214  int nummodesA,nummodesB, 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  switch(fid)
1224  {
1225  case 0:
1226  nummodesA = m_base[0]->GetNumModes();
1227  nummodesB = m_base[1]->GetNumModes();
1228  break;
1229  case 1:
1230  nummodesA = m_base[0]->GetNumModes();
1231  nummodesB = m_base[2]->GetNumModes();
1232  break;
1233  case 2:
1234  case 3:
1235  nummodesA = m_base[1]->GetNumModes();
1236  nummodesB = m_base[2]->GetNumModes();
1237  break;
1238  }
1239 
1240  bool CheckForZeroedModes = false;
1241  if(P == -1)
1242  {
1243  P = nummodesA;
1244  Q = nummodesB;
1245  }
1246  else
1247  {
1248  CheckForZeroedModes = true;
1249  }
1250 
1251  nFaceCoeffs = P*(2*Q-P+1)/2;
1252 
1253  // Allocate the map array and sign array; set sign array to ones (+)
1254  if(maparray.num_elements() != nFaceCoeffs)
1255  {
1256  maparray = Array<OneD, unsigned int>(nFaceCoeffs,1);
1257  }
1258 
1259  if(signarray.num_elements() != nFaceCoeffs)
1260  {
1261  signarray = Array<OneD, int>(nFaceCoeffs,1);
1262  }
1263  else
1264  {
1265  fill(signarray.get(),signarray.get()+nFaceCoeffs, 1 );
1266  }
1267 
1268  switch (fid)
1269  {
1270  case 0:
1271  idx = 0;
1272  for (i = 0; i < P; ++i)
1273  {
1274  for (j = 0; j < Q-i; ++j)
1275  {
1276  if ((int)faceOrient == 7 && i > 1)
1277  {
1278  signarray[idx] = (i%2 ? -1 : 1);
1279  }
1280  maparray[idx++] = GetMode(i,j,0);
1281  }
1282  }
1283  break;
1284  case 1:
1285  idx = 0;
1286  for (i = 0; i < P; ++i)
1287  {
1288  for (k = 0; k < Q-i; ++k)
1289  {
1290  if ((int)faceOrient == 7 && i > 1)
1291  {
1292  signarray[idx] = (i%2 ? -1: 1);
1293  }
1294  maparray[idx++] = GetMode(i,0,k);
1295  }
1296  }
1297  break;
1298  case 2:
1299  idx = 0;
1300  for (j = 0; j < P-1; ++j)
1301  {
1302  for (k = 0; k < Q-1-j; ++k)
1303  {
1304  if ((int)faceOrient == 7 && j > 1)
1305  {
1306  signarray[idx] = ((j+1)%2 ? -1: 1);
1307  }
1308  maparray[idx++] = GetMode(1,j,k);
1309  // Incorporate modes from zeroth plane where needed.
1310  if (j == 0 && k == 0)
1311  {
1312  maparray[idx++] = GetMode(0,0,1);
1313  }
1314  if (j == 0 && k == Q-2)
1315  {
1316  for (int r = 0; r < Q-1; ++r)
1317  {
1318  maparray[idx++] = GetMode(0,1,r);
1319  }
1320  }
1321  }
1322  }
1323  break;
1324  case 3:
1325  idx = 0;
1326  for (j = 0; j < P; ++j)
1327  {
1328  for (k = 0; k < Q-j; ++k)
1329  {
1330  if ((int)faceOrient == 7 && j > 1)
1331  {
1332  signarray[idx] = (j%2 ? -1: 1);
1333  }
1334  maparray[idx++] = GetMode(0,j,k);
1335  }
1336  }
1337  break;
1338  default:
1339  ASSERTL0(false, "Element map not available.");
1340  }
1341 
1342  if ((int)faceOrient == 7)
1343  {
1344  swap(maparray[0], maparray[Q]);
1345 
1346  for (i = 1; i < Q-1; ++i)
1347  {
1348  swap(maparray[i+1], maparray[Q+i]);
1349  }
1350  }
1351 
1352  if(CheckForZeroedModes)
1353  {
1354  // zero signmap and set maparray to zero if elemental
1355  // modes are not as large as face modesl
1356  idx = 0;
1357  for (j = 0; j < nummodesA; ++j)
1358  {
1359  idx += nummodesB-j;
1360  for (k = nummodesB-j; k < Q-j; ++k)
1361  {
1362  signarray[idx] = 0.0;
1363  maparray[idx++] = maparray[0];
1364  }
1365  }
1366 
1367  for (j = nummodesA; j < P; ++j)
1368  {
1369  for (k = 0; k < Q-j; ++k)
1370  {
1371  signarray[idx] = 0.0;
1372  maparray[idx++] = maparray[0];
1373  }
1374  }
1375  }
1376  }
1377 
1378  int StdTetExp::v_GetVertexMap(const int localVertexId, bool useCoeffPacking)
1379  {
1380  ASSERTL0((GetEdgeBasisType(localVertexId)==LibUtilities::eModified_A)||
1381  (GetEdgeBasisType(localVertexId)==LibUtilities::eModified_B)||
1382  (GetEdgeBasisType(localVertexId)==LibUtilities::eModified_C),
1383  "Mapping not defined for this type of basis");
1384 
1385  int localDOF = 0;
1386  if(useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1387  {
1388  switch(localVertexId)
1389  {
1390  case 0:
1391  {
1392  localDOF = GetMode(0,0,0);
1393  break;
1394  }
1395  case 1:
1396  {
1397  localDOF = GetMode(0,0,1);
1398  break;
1399  }
1400  case 2:
1401  {
1402  localDOF = GetMode(0,1,0);
1403  break;
1404  }
1405  case 3:
1406  {
1407  localDOF = GetMode(1,0,0);
1408  break;
1409  }
1410  default:
1411  {
1412  ASSERTL0(false,"Vertex ID must be between 0 and 3");
1413  break;
1414  }
1415  }
1416  }
1417  else
1418  {
1419  switch(localVertexId)
1420  {
1421  case 0:
1422  {
1423  localDOF = GetMode(0,0,0);
1424  break;
1425  }
1426  case 1:
1427  {
1428  localDOF = GetMode(1,0,0);
1429  break;
1430  }
1431  case 2:
1432  {
1433  localDOF = GetMode(0,1,0);
1434  break;
1435  }
1436  case 3:
1437  {
1438  localDOF = GetMode(0,0,1);
1439  break;
1440  }
1441  default:
1442  {
1443  ASSERTL0(false,"Vertex ID must be between 0 and 3");
1444  break;
1445  }
1446  }
1447 
1448  }
1449 
1450  return localDOF;
1451  }
1452 
1453  /**
1454  * Maps interior modes of an edge to the elemental modes.
1455  */
1456  void StdTetExp::v_GetEdgeInteriorMap(
1457  const int eid,
1458  const Orientation edgeOrient,
1459  Array<OneD, unsigned int> &maparray,
1460  Array<OneD, int> &signarray)
1461  {
1462  int i;
1463  const int P = m_base[0]->GetNumModes();
1464  const int Q = m_base[1]->GetNumModes();
1465  const int R = m_base[2]->GetNumModes();
1466 
1467  const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid)-2;
1468 
1469  if(maparray.num_elements() != nEdgeIntCoeffs)
1470  {
1471  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1472  }
1473  else
1474  {
1475  fill( maparray.get(), maparray.get() + nEdgeIntCoeffs, 0);
1476  }
1477 
1478  if(signarray.num_elements() != nEdgeIntCoeffs)
1479  {
1480  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1481  }
1482  else
1483  {
1484  fill( signarray.get() , signarray.get()+nEdgeIntCoeffs, 1 );
1485  }
1486 
1487  switch (eid)
1488  {
1489  case 0:
1490  for (i = 0; i < P-2; ++i)
1491  {
1492  maparray[i] = GetMode(i+2, 0, 0);
1493  }
1494  if(edgeOrient==eBackwards)
1495  {
1496  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1497  {
1498  signarray[i] = -1;
1499  }
1500  }
1501  break;
1502  case 1:
1503  for (i = 0; i < Q-2; ++i)
1504  {
1505  maparray[i] = GetMode(1, i+1, 0);
1506  }
1507  if(edgeOrient==eBackwards)
1508  {
1509  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1510  {
1511  signarray[i] = -1;
1512  }
1513  }
1514  break;
1515  case 2:
1516  for (i = 0; i < Q-2; ++i)
1517  {
1518  maparray[i] = GetMode(0, i+2, 0);
1519  }
1520  if(edgeOrient==eBackwards)
1521  {
1522  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1523  {
1524  signarray[i] = -1;
1525  }
1526  }
1527  break;
1528  case 3:
1529  for (i = 0; i < R-2; ++i)
1530  {
1531  maparray[i] = GetMode(0, 0, i+2);
1532  }
1533  if(edgeOrient==eBackwards)
1534  {
1535  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1536  {
1537  signarray[i] = -1;
1538  }
1539  }
1540  break;
1541  case 4:
1542  for (i = 0; i < R - 2; ++i)
1543  {
1544  maparray[i] = GetMode(1, 0, i+1);
1545  }
1546  if(edgeOrient==eBackwards)
1547  {
1548  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1549  {
1550  signarray[i] = -1;
1551  }
1552  }
1553  break;
1554  case 5:
1555  for (i = 0; i < R-2; ++i)
1556  {
1557  maparray[i] = GetMode(0, 1, i+1);
1558  }
1559  if(edgeOrient==eBackwards)
1560  {
1561  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1562  {
1563  signarray[i] = -1;
1564  }
1565  }
1566  break;
1567  default:
1568  ASSERTL0(false, "Edge not defined.");
1569  break;
1570  }
1571  }
1572 
1573  void StdTetExp::v_GetFaceInteriorMap(
1574  const int fid,
1575  const Orientation faceOrient,
1576  Array<OneD, unsigned int> &maparray,
1577  Array<OneD, int> &signarray)
1578  {
1579  int i, j, idx, k;
1580  const int P = m_base[0]->GetNumModes();
1581  const int Q = m_base[1]->GetNumModes();
1582  const int R = m_base[2]->GetNumModes();
1583 
1584  const int nFaceIntCoeffs = v_GetFaceIntNcoeffs(fid);
1585 
1586  if(maparray.num_elements() != nFaceIntCoeffs)
1587  {
1588  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1589  }
1590 
1591  if(signarray.num_elements() != nFaceIntCoeffs)
1592  {
1593  signarray = Array<OneD, int>(nFaceIntCoeffs,1);
1594  }
1595  else
1596  {
1597  fill( signarray.get() , signarray.get()+nFaceIntCoeffs, 1 );
1598  }
1599 
1600  switch (fid)
1601  {
1602  case 0:
1603  idx = 0;
1604  for (i = 2; i < P-1; ++i)
1605  {
1606  for (j = 1; j < Q-i; ++j)
1607  {
1608  if ((int)faceOrient == 7)
1609  {
1610  signarray[idx] = (i%2 ? -1 : 1);
1611  }
1612  maparray[idx++] = GetMode(i,j,0);
1613  }
1614  }
1615  break;
1616  case 1:
1617  idx = 0;
1618  for (i = 2; i < P; ++i)
1619  {
1620  for (k = 1; k < R-i; ++k)
1621  {
1622  if ((int)faceOrient == 7)
1623  {
1624  signarray[idx] = (i%2 ? -1: 1);
1625  }
1626  maparray[idx++] = GetMode(i,0,k);
1627  }
1628  }
1629  break;
1630  case 2:
1631  idx = 0;
1632  for (j = 1; j < Q-2; ++j)
1633  {
1634  for (k = 1; k < R-1-j; ++k)
1635  {
1636  if ((int)faceOrient == 7)
1637  {
1638  signarray[idx] = ((j+1)%2 ? -1: 1);
1639  }
1640  maparray[idx++] = GetMode(1,j,k);
1641  }
1642  }
1643  break;
1644  case 3:
1645  idx = 0;
1646  for (j = 2; j < Q-1; ++j)
1647  {
1648  for (k = 1; k < R-j; ++k)
1649  {
1650  if ((int)faceOrient == 7)
1651  {
1652  signarray[idx] = (j%2 ? -1: 1);
1653  }
1654  maparray[idx++] = GetMode(0,j,k);
1655  }
1656  }
1657  break;
1658  default:
1659  ASSERTL0(false, "Face interior map not available.");
1660  break;
1661  }
1662  }
1663 
1664  /**
1665  * List of all interior modes in the expansion.
1666  */
1667  void StdTetExp::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
1668  {
1669  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1670  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1671  "BasisType is not a boundary interior form");
1672  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
1673  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1674  "BasisType is not a boundary interior form");
1675  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1676  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1677  "BasisType is not a boundary interior form");
1678 
1679  int P = m_base[0]->GetNumModes();
1680  int Q = m_base[1]->GetNumModes();
1681  int R = m_base[2]->GetNumModes();
1682 
1683  int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1684 
1685  if(outarray.num_elements() != nIntCoeffs)
1686  {
1687  outarray = Array<OneD, unsigned int>(nIntCoeffs);
1688  }
1689 
1690  int idx = 0;
1691  for (int i = 2; i < P-2; ++i)
1692  {
1693  for (int j = 1; j < Q-i-1; ++j)
1694  {
1695  for (int k = 1; k < R-i-j; ++k)
1696  {
1697  outarray[idx++] = GetMode(i,j,k);
1698  }
1699  }
1700  }
1701  }
1702 
1703  /**
1704  * List of all boundary modes in the the expansion.
1705  */
1706  void StdTetExp::v_GetBoundaryMap(Array<OneD, unsigned int>& outarray)
1707  {
1708  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1709  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1710  "BasisType is not a boundary interior form");
1711  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
1712  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1713  "BasisType is not a boundary interior form");
1714  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1715  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1716  "BasisType is not a boundary interior form");
1717 
1718  int P = m_base[0]->GetNumModes();
1719  int Q = m_base[1]->GetNumModes();
1720  int R = m_base[2]->GetNumModes();
1721 
1722  int i,j,k;
1723  int idx = 0;
1724 
1725  int nBnd = NumBndryCoeffs();
1726 
1727  if (outarray.num_elements() != nBnd)
1728  {
1729  outarray = Array<OneD, unsigned int>(nBnd);
1730  }
1731 
1732  for (i = 0; i < P; ++i)
1733  {
1734  // First two Q-R planes are entirely boundary modes
1735  if (i < 2)
1736  {
1737  for (j = 0; j < Q-i; j++)
1738  {
1739  for (k = 0; k < R-i-j; ++k)
1740  {
1741  outarray[idx++] = GetMode(i,j,k);
1742  }
1743  }
1744  }
1745  // Remaining Q-R planes contain boundary modes on bottom and
1746  // left edge.
1747  else
1748  {
1749  for (k = 0; k < R-i; ++k)
1750  {
1751  outarray[idx++] = GetMode(i,0,k);
1752  }
1753  for (j = 1; j < Q-i; ++j)
1754  {
1755  outarray[idx++] = GetMode(i,j,0);
1756  }
1757  }
1758  }
1759  }
1760 
1761 
1762  //---------------------------------------
1763  // Wrapper functions
1764  //---------------------------------------
1765  DNekMatSharedPtr StdTetExp::v_GenMatrix(const StdMatrixKey &mkey)
1766  {
1767 
1768  MatrixType mtype = mkey.GetMatrixType();
1769 
1770  DNekMatSharedPtr Mat;
1771 
1772  switch(mtype)
1773  {
1774  case ePhysInterpToEquiSpaced:
1775  {
1776  int nq0 = m_base[0]->GetNumPoints();
1777  int nq1 = m_base[1]->GetNumPoints();
1778  int nq2 = m_base[2]->GetNumPoints();
1779  int nq;
1780 
1781  // take definition from key
1782  if(mkey.ConstFactorExists(eFactorConst))
1783  {
1784  nq = (int) mkey.GetConstFactor(eFactorConst);
1785  }
1786  else
1787  {
1788  nq = max(nq0,max(nq1,nq2));
1789  }
1790 
1791  int neq = LibUtilities::StdTetData::
1792  getNumberOfCoefficients(nq,nq,nq);
1793  Array<OneD, Array<OneD, NekDouble> > coords(neq);
1794  Array<OneD, NekDouble> coll(3);
1795  Array<OneD, DNekMatSharedPtr> I(3);
1796  Array<OneD, NekDouble> tmp(nq0);
1797 
1798  Mat = MemoryManager<DNekMat>::
1799  AllocateSharedPtr(neq, nq0 * nq1 * nq2);
1800  int cnt = 0;
1801 
1802  for(int i = 0; i < nq; ++i)
1803  {
1804  for(int j = 0; j < nq-i; ++j)
1805  {
1806  for(int k = 0; k < nq-i-j; ++k,++cnt)
1807  {
1808  coords[cnt] = Array<OneD, NekDouble>(3);
1809  coords[cnt][0] = -1.0 + 2*k/(NekDouble)(nq-1);
1810  coords[cnt][1] = -1.0 + 2*j/(NekDouble)(nq-1);
1811  coords[cnt][2] = -1.0 + 2*i/(NekDouble)(nq-1);
1812  }
1813  }
1814  }
1815 
1816  for(int i = 0; i < neq; ++i)
1817  {
1818  LocCoordToLocCollapsed(coords[i],coll);
1819 
1820  I[0] = m_base[0]->GetI(coll);
1821  I[1] = m_base[1]->GetI(coll+1);
1822  I[2] = m_base[2]->GetI(coll+2);
1823 
1824  // interpolate first coordinate direction
1825  NekDouble fac;
1826  for( int k = 0; k < nq2; ++k)
1827  {
1828  for (int j = 0; j < nq1; ++j)
1829  {
1830 
1831  fac = (I[1]->GetPtr())[j]*(I[2]->GetPtr())[k];
1832  Vmath::Smul(nq0,fac,I[0]->GetPtr(),1,tmp,1);
1833 
1834  Vmath::Vcopy(nq0, &tmp[0], 1,
1835  Mat->GetRawPtr()+
1836  k*nq0*nq1*neq+
1837  j*nq0*neq+i,neq);
1838  }
1839  }
1840  }
1841  }
1842  break;
1843  default:
1844  {
1845  Mat = StdExpansion::CreateGeneralMatrix(mkey);
1846  }
1847  break;
1848  }
1849 
1850  return Mat;
1851  }
1852 
1853  DNekMatSharedPtr StdTetExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
1854  {
1855  return v_GenMatrix(mkey);
1856  }
1857 
1858 
1859  //---------------------------------------
1860  // Private helper functions
1861  //---------------------------------------
1862 
1863  /**
1864  * @brief Compute the mode number in the expansion for a particular
1865  * tensorial combination.
1866  *
1867  * Modes are numbered with the r index travelling fastest, followed by
1868  * q and then p, and each q-r plane is of size
1869  * (Q+1)*(Q+2)/2+max(0,R-Q-p)*Q. For example, when P=2, Q=3 and R=4
1870  * the indexing inside each q-r plane (with r increasing upwards and q
1871  * to the right) is:
1872  *
1873  * p = 0: p = 1: p = 2:
1874  * ----------------------------------
1875  * 4
1876  * 3 8 17
1877  * 2 7 11 16 20 26
1878  * 1 6 10 13 15 19 22 25 28
1879  * 0 5 9 12 14 18 21 23 24 27 29
1880  *
1881  * Note that in this element, we must have that \f$ P \leq Q \leq
1882  * R\f$.
1883  */
1884  int StdTetExp::GetMode(const int I, const int J, const int K)
1885  {
1886  const int Q = m_base[1]->GetNumModes();
1887  const int R = m_base[2]->GetNumModes();
1888 
1889  int i,j,q_hat,k_hat;
1890  int cnt = 0;
1891 
1892  // Traverse to q-r plane number I
1893  for (i = 0; i < I; ++i)
1894  {
1895  // Size of triangle part
1896  q_hat = min(Q,R-i);
1897  // Size of rectangle part
1898  k_hat = max(R-Q-i,0);
1899  cnt += q_hat*(q_hat+1)/2 + k_hat*Q;
1900  }
1901 
1902  // Traverse to q column J
1903  q_hat = R-I;
1904  for (j = 0; j < J; ++j)
1905  {
1906  cnt += q_hat;
1907  q_hat--;
1908  }
1909 
1910  // Traverse up stacks to K
1911  cnt += K;
1912 
1913  return cnt;
1914  }
1915 
1916  void StdTetExp::v_MultiplyByStdQuadratureMetric(
1917  const Array<OneD, const NekDouble>& inarray,
1918  Array<OneD, NekDouble>& outarray)
1919  {
1920  int i, j;
1921 
1922  int nquad0 = m_base[0]->GetNumPoints();
1923  int nquad1 = m_base[1]->GetNumPoints();
1924  int nquad2 = m_base[2]->GetNumPoints();
1925 
1926  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1927  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1928  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
1929 
1930  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
1931  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1932 
1933  // multiply by integration constants
1934  for(i = 0; i < nquad1*nquad2; ++i)
1935  {
1936  Vmath::Vmul(nquad0,(NekDouble*)&inarray[0]+i*nquad0,1,
1937  w0.get(),1, &outarray[0]+i*nquad0,1);
1938  }
1939 
1940  switch(m_base[1]->GetPointsType())
1941  {
1942  // (1,0) Jacobi Inner product.
1943  case LibUtilities::eGaussRadauMAlpha1Beta0:
1944  for(j = 0; j < nquad2; ++j)
1945  {
1946  for(i = 0; i < nquad1; ++i)
1947  {
1948  Blas::Dscal(nquad0,0.5*w1[i], &outarray[0]+i*nquad0+
1949  j*nquad0*nquad1,1);
1950  }
1951  }
1952  break;
1953 
1954  default:
1955  for(j = 0; j < nquad2; ++j)
1956  {
1957  for(i = 0; i < nquad1; ++i)
1958  {
1959  Blas::Dscal(nquad0,
1960  0.5*(1-z1[i])*w1[i],
1961  &outarray[0]+i*nquad0 + j*nquad0*nquad1,
1962  1 );
1963  }
1964  }
1965  break;
1966  }
1967 
1968  switch(m_base[2]->GetPointsType())
1969  {
1970  // (2,0) Jacobi inner product.
1971  case LibUtilities::eGaussRadauMAlpha2Beta0:
1972  for(i = 0; i < nquad2; ++i)
1973  {
1974  Blas::Dscal(nquad0*nquad1, 0.25*w2[i],
1975  &outarray[0]+i*nquad0*nquad1, 1);
1976  }
1977  break;
1978  // (1,0) Jacobi inner product.
1979  case LibUtilities::eGaussRadauMAlpha1Beta0:
1980  for(i = 0; i < nquad2; ++i)
1981  {
1982  Blas::Dscal(nquad0*nquad1, 0.25*(1-z2[i])*w2[i],
1983  &outarray[0]+i*nquad0*nquad1, 1);
1984  }
1985  break;
1986  default:
1987  for(i = 0; i < nquad2; ++i)
1988  {
1989  Blas::Dscal(nquad0*nquad1,0.25*(1-z2[i])*(1-z2[i])*w2[i],
1990  &outarray[0]+i*nquad0*nquad1,1);
1991  }
1992  break;
1993  }
1994  }
1995 
1996  void StdTetExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
1997  const StdMatrixKey &mkey)
1998  {
1999  //To do : 1) add a test to ensure 0 \leq SvvCutoff \leq 1.
2000  // 2) check if the transfer function needs an analytical
2001  // Fourier transform.
2002  // 3) if it doesn't : find a transfer function that renders
2003  // the if( cutoff_a ...) useless to reduce computational
2004  // cost.
2005  // 4) add SVVDiffCoef to both models!!
2006 
2007  int qa = m_base[0]->GetNumPoints();
2008  int qb = m_base[1]->GetNumPoints();
2009  int qc = m_base[2]->GetNumPoints();
2010  int nmodes_a = m_base[0]->GetNumModes();
2011  int nmodes_b = m_base[1]->GetNumModes();
2012  int nmodes_c = m_base[2]->GetNumModes();
2013 
2014  // Declare orthogonal basis.
2015  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
2016  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
2017  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
2018 
2019  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
2020  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B,nmodes_b,pb);
2021  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_C,nmodes_c,pc);
2022 
2023  StdTetExp OrthoExp(Ba,Bb,Bc);
2024 
2025 
2026  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2027  int i,j,k,cnt = 0;
2028 
2029  //SVV filter paramaters (how much added diffusion relative to physical one
2030  // and fraction of modes from which you start applying this added diffusion)
2031  //
2032  NekDouble SvvDiffCoeff = mkey.GetConstFactor(StdRegions::eFactorSVVDiffCoeff);
2033  NekDouble SVVCutOff = mkey.GetConstFactor(StdRegions::eFactorSVVCutoffRatio);
2034 
2035 
2036  //Defining the cut of mode
2037  int cutoff_a = (int) (SVVCutOff*nmodes_a);
2038  int cutoff_b = (int) (SVVCutOff*nmodes_b);
2039  int cutoff_c = (int) (SVVCutOff*nmodes_c);
2040  int nmodes = min(min(nmodes_a,nmodes_b),nmodes_c);
2041  NekDouble cutoff = min(min(cutoff_a,cutoff_b),cutoff_c);
2042  NekDouble epsilon = 1;
2043 
2044  // project onto physical space.
2045  OrthoExp.FwdTrans(array,orthocoeffs);
2046 
2047  //------"New" Version August 22nd '13--------------------
2048  for(i = 0; i < nmodes_a; ++i)
2049  {
2050  for(j = 0; j < nmodes_b-i; ++j)
2051  {
2052  for(k = 0; k < nmodes_c-i-j; ++k)
2053  {
2054  if(i + j + k >= cutoff)
2055  {
2056  orthocoeffs[cnt] *= ((SvvDiffCoeff)*exp(-(i+j+k-nmodes)*(i+j+k-nmodes)/((NekDouble)((i+j+k-cutoff+epsilon)*(i+j+k-cutoff+epsilon)))));
2057  }
2058  else
2059  {
2060  orthocoeffs[cnt] *= 0.0;
2061  }
2062  cnt++;
2063  }
2064  }
2065  }
2066  // backward transform to physical space
2067  OrthoExp.BwdTrans(orthocoeffs,array);
2068  }
2069 
2070 
2071  void StdTetExp::v_ReduceOrderCoeffs(
2072  int numMin,
2073  const Array<OneD, const NekDouble> &inarray,
2074  Array<OneD, NekDouble> &outarray)
2075  {
2076  int nquad0 = m_base[0]->GetNumPoints();
2077  int nquad1 = m_base[1]->GetNumPoints();
2078  int nquad2 = m_base[2]->GetNumPoints();
2079  int nqtot = nquad0 * nquad1 * nquad2;
2080  int nmodes0 = m_base[0]->GetNumModes();
2081  int nmodes1 = m_base[1]->GetNumModes();
2082  int nmodes2 = m_base[2]->GetNumModes();
2083  int numMax = nmodes0;
2084 
2085  Array<OneD, NekDouble> coeff (m_ncoeffs);
2086  Array<OneD, NekDouble> coeff_tmp1(m_ncoeffs, 0.0);
2087  Array<OneD, NekDouble> coeff_tmp2(m_ncoeffs, 0.0);
2088  Array<OneD, NekDouble> phys_tmp (nqtot, 0.0);
2089  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
2090 
2091  Vmath::Vcopy(m_ncoeffs,inarray,1,coeff_tmp2,1);
2092 
2093  const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
2094  const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
2095  const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
2096 
2097  LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A,
2098  nmodes0, Pkey0);
2099  LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_B,
2100  nmodes1, Pkey1);
2101  LibUtilities::BasisKey bortho2(LibUtilities::eOrtho_C,
2102  nmodes2, Pkey2);
2103 
2104  Vmath::Zero(m_ncoeffs, coeff_tmp2, 1);
2105 
2106  StdRegions::StdTetExpSharedPtr OrthoTetExp;
2107  OrthoTetExp = MemoryManager<StdRegions::StdTetExp>
2108  ::AllocateSharedPtr(bortho0, bortho1, bortho2);
2109 
2110  BwdTrans(inarray,phys_tmp);
2111  OrthoTetExp->FwdTrans(phys_tmp, coeff);
2112 
2113  Vmath::Zero(m_ncoeffs,outarray,1);
2114 
2115  // filtering
2116  int cnt = 0;
2117  for (int u = 0; u < numMin; ++u)
2118  {
2119  for (int i = 0; i < numMin-u; ++i)
2120  {
2121  Vmath::Vcopy(numMin - u - i, tmp = coeff + cnt, 1,
2122  tmp2 = coeff_tmp1 + cnt, 1);
2123  cnt += numMax - u - i;
2124  }
2125  for (int i = numMin; i < numMax-u; ++i)
2126  {
2127  cnt += numMax - u - i;
2128  }
2129  }
2130 
2131  OrthoTetExp->BwdTrans(coeff_tmp1,phys_tmp);
2132  FwdTrans(phys_tmp, outarray);
2133  }
2134 
2135 
2136  void StdTetExp::v_GetSimplexEquiSpacedConnectivity(
2137  Array<OneD, int> &conn,
2138  bool standard)
2139  {
2140  int np0 = m_base[0]->GetNumPoints();
2141  int np1 = m_base[1]->GetNumPoints();
2142  int np2 = m_base[2]->GetNumPoints();
2143  int np = max(np0,max(np1,np2));
2144 
2145 
2146  conn = Array<OneD, int>(4*(np-1)*(np-1)*(np-1));
2147 
2148  int row = 0;
2149  int rowp1 = 0;
2150  int plane = 0;
2151  int row1 = 0;
2152  int row1p1 = 0;
2153  int planep1= 0;
2154  int cnt = 0;
2155  for(int i = 0; i < np-1; ++i)
2156  {
2157  planep1 += (np-i)*(np-i+1)/2;
2158  row = 0; // current plane row offset
2159  rowp1 = 0; // current plane row plus one offset
2160  row1 = 0; // next plane row offset
2161  row1p1 = 0; // nex plane row plus one offset
2162  for(int j = 0; j < np-i-1; ++j)
2163  {
2164  rowp1 += np-i-j;
2165  row1p1 += np-i-j-1;
2166  for(int k = 0; k < np-i-j-2; ++k)
2167  {
2168  conn[cnt++] = plane + row +k+1;
2169  conn[cnt++] = plane + row +k;
2170  conn[cnt++] = plane + rowp1 +k;
2171  conn[cnt++] = planep1 + row1 +k;
2172 
2173  conn[cnt++] = plane + row +k+1;
2174  conn[cnt++] = plane + rowp1 +k+1;
2175  conn[cnt++] = planep1 + row1 +k+1;
2176  conn[cnt++] = planep1 + row1 +k;
2177 
2178  conn[cnt++] = plane + rowp1 +k+1;
2179  conn[cnt++] = plane + row +k+1;
2180  conn[cnt++] = plane + rowp1 +k;
2181  conn[cnt++] = planep1 + row1 +k;
2182 
2183  conn[cnt++] = planep1 + row1 +k;
2184  conn[cnt++] = planep1 + row1p1+k;
2185  conn[cnt++] = plane + rowp1 +k;
2186  conn[cnt++] = plane + rowp1 +k+1;
2187 
2188  conn[cnt++] = planep1 + row1 +k;
2189  conn[cnt++] = planep1 + row1p1+k;
2190  conn[cnt++] = planep1 + row1 +k+1;
2191  conn[cnt++] = plane + rowp1 +k+1;
2192 
2193  if(k < np-i-j-3)
2194  {
2195  conn[cnt++] = plane + rowp1 +k+1;
2196  conn[cnt++] = planep1 + row1p1 +k+1;
2197  conn[cnt++] = planep1 + row1 +k+1;
2198  conn[cnt++] = planep1 + row1p1 +k;
2199  }
2200  }
2201 
2202  conn[cnt++] = plane + row + np-i-j-1;
2203  conn[cnt++] = plane + row + np-i-j-2;
2204  conn[cnt++] = plane + rowp1 + np-i-j-2;
2205  conn[cnt++] = planep1 + row1 + np-i-j-2;
2206 
2207  row += np-i-j;
2208  row1 += np-i-j-1;
2209  }
2210  plane += (np-i)*(np-i+1)/2;
2211  }
2212  }
2213 
2214  }//end namespace
2215 }//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:1884
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:1765
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:2136
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:942
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:406
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:1173
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_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:1456
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:700
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::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:131
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:413
Nektar::StdRegions::StdExpansion3D::PhysTensorDeriv
void PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d1, Array< OneD, NekDouble > &outarray_d2, Array< OneD, NekDouble > &outarray_d3)
Calculate the 3D derivative in the local tensor/collapsed coordinate at the physical points...
Definition: StdExpansion3D.cpp:68
Nektar::StdRegions::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:1706
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:1853
Nektar::StdRegions::StdTetExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdTetExp.cpp:1667
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:1916
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:44
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_GetFaceToElementMap
virtual void v_GetFaceToElementMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int P=-1, int Q=-1)
Definition: StdTetExp.cpp:1206
Nektar::StdRegions::StdTetExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdTetExp.cpp:1996
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:1573
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:1378
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:525
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:1422
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
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:2071
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:1905
Generated on Wed Mar 23 2016 20:14:23 for Nektar++ by   doxygen 1.8.8