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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Header field for tetrahedral routines built upon
32 // StdExpansion3D
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <boost/core/ignore_unused.hpp>
37 
38 #include <StdRegions/StdTetExp.h>
39 
40 using namespace std;
41 
42 namespace Nektar
43 {
44  namespace StdRegions
45  {
46  StdTetExp::StdTetExp()
47  {
48  }
49 
50 
51  StdTetExp::StdTetExp(const LibUtilities::BasisKey &Ba,
52  const LibUtilities::BasisKey &Bb,
53  const LibUtilities::BasisKey &Bc):
54  StdExpansion(LibUtilities::StdTetData::getNumberOfCoefficients(
55  Ba.GetNumModes(),
56  Bb.GetNumModes(),
57  Bc.GetNumModes()),
58  3, Ba, Bb, Bc),
59  StdExpansion3D(LibUtilities::StdTetData::getNumberOfCoefficients(
60  Ba.GetNumModes(),
61  Bb.GetNumModes(),
62  Bc.GetNumModes()),
63  Ba, Bb, Bc)
64  {
65  ASSERTL0(Ba.GetNumModes() <= Bb.GetNumModes(),
66  "order in 'a' direction is higher than order "
67  "in 'b' direction");
68  ASSERTL0(Ba.GetNumModes() <= Bc.GetNumModes(),
69  "order in 'a' direction is higher than order "
70  "in 'c' direction");
71  ASSERTL0(Bb.GetNumModes() <= Bc.GetNumModes(),
72  "order in 'b' direction is higher than order "
73  "in 'c' direction");
74  }
75 
76  StdTetExp::StdTetExp(const StdTetExp &T):
77  StdExpansion(T),
78  StdExpansion3D(T)
79  {
80  }
81 
82 
83  StdTetExp::~StdTetExp()
84  {
85  }
86 
87  //----------------------------
88  // Differentiation Methods
89  //----------------------------
90 
91  /**
92  * \brief Calculate the derivative of the physical points
93  *
94  * The derivative is evaluated at the nodal physical points.
95  * Derivatives with respect to the local Cartesian coordinates
96  *
97  * \f$\begin{Bmatrix} \frac {\partial} {\partial \xi_1} \\ \frac
98  * {\partial} {\partial \xi_2} \\ \frac {\partial} {\partial \xi_3}
99  * \end{Bmatrix} = \begin{Bmatrix} \frac 4 {(1-\eta_2)(1-\eta_3)}
100  * \frac \partial {\partial \eta_1} \ \ \frac {2(1+\eta_1)}
101  * {(1-\eta_2)(1-\eta_3)} \frac \partial {\partial \eta_1} + \frac 2
102  * {1-\eta_3} \frac \partial {\partial \eta_3} \\ \frac {2(1 +
103  * \eta_1)} {2(1 - \eta_2)(1-\eta_3)} \frac \partial {\partial \eta_1}
104  * + \frac {1 + \eta_2} {1 - \eta_3} \frac \partial {\partial \eta_2}
105  * + \frac \partial {\partial \eta_3} \end{Bmatrix}\f$
106  **/
107  void StdTetExp::v_PhysDeriv(
108  const Array<OneD, const NekDouble>& inarray,
109  Array<OneD, NekDouble>& out_dxi0,
110  Array<OneD, NekDouble>& out_dxi1,
111  Array<OneD, NekDouble>& out_dxi2 )
112  {
113  int Q0 = m_base[0]->GetNumPoints();
114  int Q1 = m_base[1]->GetNumPoints();
115  int Q2 = m_base[2]->GetNumPoints();
116  int Qtot = Q0*Q1*Q2;
117 
118  // Compute the physical derivative
119  Array<OneD, NekDouble> out_dEta0(3*Qtot,0.0);
120  Array<OneD, NekDouble> out_dEta1 = out_dEta0 + Qtot;
121  Array<OneD, NekDouble> out_dEta2 = out_dEta1 + Qtot;
122 
123  bool Do_2 = (out_dxi2.size() > 0)? true:false;
124  bool Do_1 = (out_dxi1.size() > 0)? true:false;
125 
126  if(Do_2) // Need all local derivatives
127  {
128  PhysTensorDeriv(inarray, out_dEta0, out_dEta1, out_dEta2);
129  }
130  else if (Do_1) // Need 0 and 1 derivatives
131  {
132  PhysTensorDeriv(inarray, out_dEta0, out_dEta1, NullNekDouble1DArray);
133  }
134  else // Only need Eta0 derivaitve
135  {
136  PhysTensorDeriv(inarray, out_dEta0, NullNekDouble1DArray,
137  NullNekDouble1DArray);
138  }
139 
140  Array<OneD, const NekDouble> eta_0, eta_1, eta_2;
141  eta_0 = m_base[0]->GetZ();
142  eta_1 = m_base[1]->GetZ();
143  eta_2 = m_base[2]->GetZ();
144 
145  // calculate 2.0/((1-eta_1)(1-eta_2)) Out_dEta0
146 
147  NekDouble *dEta0 = &out_dEta0[0];
148  NekDouble fac;
149  for(int k=0; k< Q2; ++k)
150  {
151  for(int j=0; j<Q1; ++j,dEta0+=Q0)
152  {
153  Vmath::Smul(Q0,2.0/(1.0-eta_1[j]),dEta0,1,dEta0,1);
154  }
155  fac = 1.0/(1.0-eta_2[k]);
156  Vmath::Smul(Q0*Q1,fac,&out_dEta0[0]+k*Q0*Q1,1,&out_dEta0[0]+k*Q0*Q1,1);
157  }
158 
159  if (out_dxi0.size() > 0)
160  {
161  // out_dxi0 = 4.0/((1-eta_1)(1-eta_2)) Out_dEta0
162  Vmath::Smul(Qtot,2.0,out_dEta0,1,out_dxi0,1);
163  }
164 
165  if (Do_1||Do_2)
166  {
167  Array<OneD, NekDouble> Fac0(Q0);
168  Vmath::Sadd(Q0,1.0,eta_0,1,Fac0,1);
169 
170 
171  // calculate 2.0*(1+eta_0)/((1-eta_1)(1-eta_2)) Out_dEta0
172  for(int k = 0; k < Q1*Q2; ++k)
173  {
174  Vmath::Vmul(Q0,&Fac0[0],1,&out_dEta0[0]+k*Q0,1,&out_dEta0[0]+k*Q0,1);
175  }
176  // calculate 2/(1.0-eta_2) out_dEta1
177  for(int k = 0; k < Q2; ++k)
178  {
179  Vmath::Smul(Q0*Q1,2.0/(1.0-eta_2[k]),&out_dEta1[0]+k*Q0*Q1,1,
180  &out_dEta1[0]+k*Q0*Q1,1);
181  }
182 
183  if(Do_1)
184  {
185  // calculate out_dxi1 = 2.0(1+eta_0)/((1-eta_1)(1-eta_2)) Out_dEta0
186  // + 2/(1.0-eta_2) out_dEta1
187  Vmath::Vadd(Qtot,out_dEta0,1,out_dEta1,1,out_dxi1,1);
188  }
189 
190 
191  if(Do_2)
192  {
193  // calculate (1 + eta_1)/(1 -eta_2)*out_dEta1
194  NekDouble *dEta1 = &out_dEta1[0];
195  for(int k=0; k< Q2; ++k)
196  {
197  for(int j=0; j<Q1; ++j,dEta1+=Q0)
198  {
199  Vmath::Smul(Q0,(1.0+eta_1[j])/2.0,dEta1,1,dEta1,1);
200  }
201  }
202 
203  // calculate out_dxi2 =
204  // 2.0(1+eta_0)/((1-eta_1)(1-eta_2)) Out_dEta0 +
205  // (1 + eta_1)/(1 -eta_2)*out_dEta1 + out_dEta2
206  Vmath::Vadd(Qtot,out_dEta0,1,out_dEta1,1,out_dxi2,1);
207  Vmath::Vadd(Qtot,out_dEta2,1,out_dxi2 ,1,out_dxi2,1);
208 
209  }
210  }
211  }
212 
213  /**
214  * @param dir Direction in which to compute derivative.
215  * Valid values are 0, 1, 2.
216  * @param inarray Input array.
217  * @param outarray Output array.
218  */
219  void StdTetExp::v_PhysDeriv(
220  const int dir,
221  const Array<OneD, const NekDouble>& inarray,
222  Array<OneD, NekDouble>& outarray)
223  {
224  switch(dir)
225  {
226  case 0:
227  {
228  v_PhysDeriv(inarray, outarray, NullNekDouble1DArray,
229  NullNekDouble1DArray);
230  break;
231  }
232  case 1:
233  {
234  v_PhysDeriv(inarray, NullNekDouble1DArray, outarray,
235  NullNekDouble1DArray);
236  break;
237  }
238  case 2:
239  {
240  v_PhysDeriv(inarray, NullNekDouble1DArray,
241  NullNekDouble1DArray, outarray);
242  break;
243  }
244  default:
245  {
246  ASSERTL1(false, "input dir is out of range");
247  }
248  break;
249  }
250  }
251 
252  void StdTetExp::v_StdPhysDeriv(
253  const Array<OneD, const NekDouble>& inarray,
254  Array<OneD, NekDouble>& out_d0,
255  Array<OneD, NekDouble>& out_d1,
256  Array<OneD, NekDouble>& out_d2)
257  {
258  StdTetExp::v_PhysDeriv(inarray, out_d0, out_d1, out_d2);
259  }
260 
261  void StdTetExp::v_StdPhysDeriv(
262  const int dir,
263  const Array<OneD, const NekDouble>& inarray,
264  Array<OneD, NekDouble>& outarray)
265  {
266  StdTetExp::v_PhysDeriv(dir, inarray, outarray);
267  }
268 
269  //---------------------------------------
270  // Transforms
271  //---------------------------------------
272 
273  /**
274  * @note 'r' (base[2]) runs fastest in this element
275  *
276  * \f$ u^{\delta} (\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{m(pqr)} \hat
277  * u_{pqr} \phi_{pqr} (\xi_{1i}, \xi_{2j}, \xi_{3k})\f$
278  *
279  * Backward transformation is three dimensional tensorial expansion
280  * \f$ u (\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{p=0}^{Q_x} \psi_p^a
281  * (\xi_{1i}) \lbrace { \sum_{q=0}^{Q_y} \psi_{pq}^b (\xi_{2j})
282  * \lbrace { \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{pqr}^c (\xi_{3k})
283  * \rbrace} \rbrace}. \f$ And sumfactorizing step of the form is as:\\
284  *
285  * \f$ f_{pq} (\xi_{3k}) = \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{pqr}^c
286  * (\xi_{3k}),\\
287  *
288  * g_{p} (\xi_{2j}, \xi_{3k}) = \sum_{r=0}^{Q_y} \psi_{pq}^b
289  * (\xi_{2j}) f_{pq} (\xi_{3k}),\\
290  *
291  * u(\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{p=0}^{Q_x} \psi_{p}^a
292  * (\xi_{1i}) g_{p} (\xi_{2j}, \xi_{3k}). \f$
293  */
294  void StdTetExp::v_BwdTrans(
295  const Array<OneD, const NekDouble>& inarray,
296  Array<OneD, NekDouble>& outarray)
297  {
298  ASSERTL1((m_base[1]->GetBasisType() != LibUtilities::eOrtho_B) ||
299  (m_base[1]->GetBasisType() != LibUtilities::eModified_B),
300  "Basis[1] is not a general tensor type");
301 
302  ASSERTL1((m_base[2]->GetBasisType() != LibUtilities::eOrtho_C) ||
303  (m_base[2]->GetBasisType() != LibUtilities::eModified_C),
304  "Basis[2] is not a general tensor type");
305 
306  if(m_base[0]->Collocation() && m_base[1]->Collocation()
307  && m_base[2]->Collocation())
308  {
309  Vmath::Vcopy(m_base[0]->GetNumPoints()
310  * m_base[1]->GetNumPoints()
311  * m_base[2]->GetNumPoints(),
312  inarray, 1, outarray, 1);
313  }
314  else
315  {
316  StdTetExp::v_BwdTrans_SumFac(inarray,outarray);
317  }
318  }
319 
320 
321  /**
322  * Sum-factorisation implementation of the BwdTrans operation.
323  */
324  void StdTetExp::v_BwdTrans_SumFac(
325  const Array<OneD, const NekDouble>& inarray,
326  Array<OneD, NekDouble>& outarray)
327  {
328  int nquad1 = m_base[1]->GetNumPoints();
329  int nquad2 = m_base[2]->GetNumPoints();
330  int order0 = m_base[0]->GetNumModes();
331  int order1 = m_base[1]->GetNumModes();
332 
333  Array<OneD, NekDouble> wsp(nquad2*order0*(2*order1-order0+1)/2+
334  nquad2*nquad1*order0);
335 
336  BwdTrans_SumFacKernel(m_base[0]->GetBdata(),
337  m_base[1]->GetBdata(),
338  m_base[2]->GetBdata(),
339  inarray,outarray,wsp,
340  true,true,true);
341  }
342 
343 
344  /**
345  * @param base0 x-dirn basis matrix
346  * @param base1 y-dirn basis matrix
347  * @param base2 z-dirn basis matrix
348  * @param inarray Input vector of modes.
349  * @param outarray Output vector of physical space data.
350  * @param wsp Workspace of size Q_x*P_z*(P_y+Q_y)
351  * @param doCheckCollDir0 Check for collocation of basis.
352  * @param doCheckCollDir1 Check for collocation of basis.
353  * @param doCheckCollDir2 Check for collocation of basis.
354  * @todo Account for some directions being collocated. See
355  * StdQuadExp as an example.
356  */
357  void StdTetExp::v_BwdTrans_SumFacKernel(
358  const Array<OneD, const NekDouble>& base0,
359  const Array<OneD, const NekDouble>& base1,
360  const Array<OneD, const NekDouble>& base2,
361  const Array<OneD, const NekDouble>& inarray,
362  Array<OneD, NekDouble>& outarray,
363  Array<OneD, NekDouble>& wsp,
364  bool doCheckCollDir0,
365  bool doCheckCollDir1,
366  bool doCheckCollDir2)
367  {
368  boost::ignore_unused(doCheckCollDir0, doCheckCollDir1,
369  doCheckCollDir2);
370 
371  int nquad0 = m_base[0]->GetNumPoints();
372  int nquad1 = m_base[1]->GetNumPoints();
373  int nquad2 = m_base[2]->GetNumPoints();
374 
375  int order0 = m_base[0]->GetNumModes();
376  int order1 = m_base[1]->GetNumModes();
377  int order2 = m_base[2]->GetNumModes();
378 
379  Array<OneD, NekDouble > tmp = wsp;
380  Array<OneD, NekDouble > tmp1 = tmp + nquad2*order0*(2*order1-order0+1)/2;
381 
382  int i, j, mode, mode1, cnt;
383 
384  // Perform summation over '2' direction
385  mode = mode1 = cnt = 0;
386  for(i = 0; i < order0; ++i)
387  {
388  for(j = 0; j < order1-i; ++j, ++cnt)
389  {
390  Blas::Dgemv('N', nquad2, order2-i-j,
391  1.0, base2.get()+mode*nquad2, nquad2,
392  inarray.get()+mode1, 1,
393  0.0, tmp.get()+cnt*nquad2, 1);
394  mode += order2-i-j;
395  mode1 += order2-i-j;
396  }
397  //increment mode in case order1!=order2
398  for(j = order1-i; j < order2-i; ++j)
399  {
400  mode += order2-i-j;
401  }
402  }
403 
404  // fix for modified basis by adding split of top singular
405  // vertex mode - currently (1+c)/2 x (1-b)/2 x (1-a)/2
406  // component is evaluated
407  if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
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]+nquad2,1);
412 
413  // top singular vertex - (1+c)/2 x (1-b)/2 x (1+a)/2 component
414  Blas::Daxpy(nquad2,inarray[1],base2.get()+nquad2,1,
415  &tmp[0]+order1*nquad2,1);
416  }
417 
418  // Perform summation over '1' direction
419  mode = 0;
420  for(i = 0; i < order0; ++i)
421  {
422  Blas::Dgemm('N', 'T', nquad1, nquad2, order1-i,
423  1.0, base1.get()+mode*nquad1, nquad1,
424  tmp.get()+mode*nquad2, nquad2,
425  0.0, tmp1.get()+i*nquad1*nquad2, nquad1);
426  mode += order1-i;
427  }
428 
429  // fix for modified basis by adding additional split of
430  // top and base singular vertex modes as well as singular
431  // edge
432  if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
433  {
434  // use tmp to sort out singular vertices and
435  // singular edge components with (1+b)/2 (1+a)/2 form
436  for(i = 0; i < nquad2; ++i)
437  {
438  Blas::Daxpy(nquad1,tmp[nquad2+i], base1.get()+nquad1,1,
439  &tmp1[nquad1*nquad2]+i*nquad1,1);
440  }
441  }
442 
443  // Perform summation over '0' direction
444  Blas::Dgemm('N', 'T', nquad0, nquad1*nquad2, order0,
445  1.0, base0.get(), nquad0,
446  tmp1.get(), nquad1*nquad2,
447  0.0, outarray.get(), nquad0);
448  }
449 
450 
451  /**
452  * @param inarray array of physical quadrature points to be
453  * transformed.
454  * @param outarray updated array of expansion coefficients.
455  */
456  void StdTetExp::v_FwdTrans(const Array<OneD, const NekDouble>& inarray,
457  Array<OneD, NekDouble> &outarray)
458  { //int numMax = nmodes0;
459  v_IProductWRTBase(inarray,outarray);
460 
461  // get Mass matrix inverse
462  StdMatrixKey masskey(eInvMass,DetShapeType(),*this);
463  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
464 
465  // copy inarray in case inarray == outarray
466  DNekVec in (m_ncoeffs, outarray);
467  DNekVec out(m_ncoeffs, outarray, eWrapper);
468 
469  out = (*matsys)*in;
470  }
471 
472 
473  //---------------------------------------
474  // Inner product functions
475  //---------------------------------------
476 
477  /**
478  * \f$ \begin{array}{rcl} I_{pqr} = (\phi_{pqr}, u)_{\delta} & = &
479  * \sum_{i=0}^{nq_0} \sum_{j=0}^{nq_1} \sum_{k=0}^{nq_2} \psi_{p}^{a}
480  * (\eta_{1i}) \psi_{pq}^{b} (\eta_{2j}) \psi_{pqr}^{c} (\eta_{3k})
481  * w_i w_j w_k u(\eta_{1,i} \eta_{2,j} \eta_{3,k}) J_{i,j,k}\\ & = &
482  * \sum_{i=0}^{nq_0} \psi_p^a(\eta_{1,i}) \sum_{j=0}^{nq_1}
483  * \psi_{pq}^b(\eta_{2,j}) \sum_{k=0}^{nq_2} \psi_{pqr}^c
484  * u(\eta_{1i},\eta_{2j},\eta_{3k}) J_{i,j,k} \end{array} \f$ \n
485  *
486  * where
487  *
488  * \f$ \phi_{pqr} (\xi_1 , \xi_2 , \xi_3) = \psi_p^a (\eta_1)
489  * \psi_{pq}^b (\eta_2) \psi_{pqr}^c (\eta_3) \f$
490  *
491  * which can be implemented as \n \f$f_{pqr} (\xi_{3k}) =
492  * \sum_{k=0}^{nq_3} \psi_{pqr}^c u(\eta_{1i},\eta_{2j},\eta_{3k})
493  *
494  * J_{i,j,k} = {\bf B_3 U} \f$ \n
495  *
496  * \f$ g_{pq} (\xi_{3k}) = \sum_{j=0}^{nq_1} \psi_{pq}^b (\xi_{2j})
497  * f_{pqr} (\xi_{3k}) = {\bf B_2 F} \f$ \n
498  *
499  * \f$ (\phi_{pqr}, u)_{\delta} = \sum_{k=0}^{nq_0} \psi_{p}^a
500  * (\xi_{3k}) g_{pq} (\xi_{3k}) = {\bf B_1 G} \f$
501  *
502  * @param inarray Function evaluated at physical collocation
503  * points.
504  * @param outarray Inner product with respect to each basis
505  * function over the element.
506  */
507  void StdTetExp::v_IProductWRTBase(
508  const Array<OneD, const NekDouble>& inarray,
509  Array<OneD, NekDouble> & outarray)
510  {
511  ASSERTL1( (m_base[1]->GetBasisType() != LibUtilities::eOrtho_B) ||
512  (m_base[1]->GetBasisType() != LibUtilities::eModified_B),
513  "Basis[1] is not a general tensor type");
514 
515  ASSERTL1( (m_base[2]->GetBasisType() != LibUtilities::eOrtho_C) ||
516  (m_base[2]->GetBasisType() != LibUtilities::eModified_C),
517  "Basis[2] is not a general tensor type");
518 
519  if(m_base[0]->Collocation() && m_base[1]->Collocation())
520  {
521  MultiplyByQuadratureMetric(inarray,outarray);
522  }
523  else
524  {
525  StdTetExp::v_IProductWRTBase_SumFac(inarray,outarray);
526  }
527  }
528 
529 
530  void StdTetExp::v_IProductWRTBase_MatOp(
531  const Array<OneD, const NekDouble>& inarray,
532  Array<OneD, NekDouble>& outarray)
533  {
534  int nq = GetTotPoints();
535  StdMatrixKey iprodmatkey(eIProductWRTBase,DetShapeType(),*this);
536  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
537 
538  Blas::Dgemv('N',m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
539  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
540  }
541 
542 
543  /**
544  * @param inarray Function evaluated at physical collocation
545  * points.
546  * @param outarray Inner product with respect to each basis
547  * function over the element.
548  */
549  void StdTetExp::v_IProductWRTBase_SumFac(
550  const Array<OneD, const NekDouble>& inarray,
551  Array<OneD, NekDouble>& outarray,
552  bool multiplybyweights)
553  {
554  int nquad0 = m_base[0]->GetNumPoints();
555  int nquad1 = m_base[1]->GetNumPoints();
556  int nquad2 = m_base[2]->GetNumPoints();
557  int order0 = m_base[0]->GetNumModes();
558  int order1 = m_base[1]->GetNumModes();
559 
560  Array<OneD, NekDouble> wsp (nquad1*nquad2*order0 +
561  nquad2*order0*(2*order1-order0+1)/2);
562 
563  if(multiplybyweights)
564  {
565  Array<OneD, NekDouble> tmp (nquad0*nquad1*nquad2);
566  MultiplyByQuadratureMetric(inarray, tmp);
567 
568  StdTetExp::IProductWRTBase_SumFacKernel(
569  m_base[0]->GetBdata(),
570  m_base[1]->GetBdata(),
571  m_base[2]->GetBdata(),
572  tmp, outarray, wsp, true, true, true);
573  }
574  else
575  {
576  StdTetExp::IProductWRTBase_SumFacKernel(
577  m_base[0]->GetBdata(),
578  m_base[1]->GetBdata(),
579  m_base[2]->GetBdata(),
580  inarray, outarray, wsp, true, true, true);
581  }
582  }
583 
584 
585  void StdTetExp::v_IProductWRTBase_SumFacKernel(
586  const Array<OneD, const NekDouble>& base0,
587  const Array<OneD, const NekDouble>& base1,
588  const Array<OneD, const NekDouble>& base2,
589  const Array<OneD, const NekDouble>& inarray,
590  Array<OneD, NekDouble> &outarray,
591  Array<OneD, NekDouble> &wsp,
592  bool doCheckCollDir0,
593  bool doCheckCollDir1,
594  bool doCheckCollDir2)
595  {
596  boost::ignore_unused(doCheckCollDir0, doCheckCollDir1,
597  doCheckCollDir2);
598 
599  int nquad0 = m_base[0]->GetNumPoints();
600  int nquad1 = m_base[1]->GetNumPoints();
601  int nquad2 = m_base[2]->GetNumPoints();
602 
603  int order0 = m_base[0]->GetNumModes();
604  int order1 = m_base[1]->GetNumModes();
605  int order2 = m_base[2]->GetNumModes();
606 
607  Array<OneD, NekDouble > tmp1 = wsp;
608  Array<OneD, NekDouble > tmp2 = wsp + nquad1*nquad2*order0;
609 
610  int i,j, mode,mode1, cnt;
611 
612  // Inner product with respect to the '0' direction
613  Blas::Dgemm('T', 'N', nquad1*nquad2, order0, nquad0,
614  1.0, inarray.get(), nquad0,
615  base0.get(), nquad0,
616  0.0, tmp1.get(), nquad1*nquad2);
617 
618  // Inner product with respect to the '1' direction
619  for(mode=i=0; i < order0; ++i)
620  {
621  Blas::Dgemm('T', 'N', nquad2, order1-i, nquad1,
622  1.0, tmp1.get()+i*nquad1*nquad2, nquad1,
623  base1.get()+mode*nquad1, nquad1,
624  0.0, tmp2.get()+mode*nquad2, nquad2);
625  mode += order1-i;
626  }
627 
628  // fix for modified basis for base singular vertex
629  if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
630  {
631  //base singular vertex and singular edge (1+b)/2
632  //(1+a)/2 components (makes tmp[nquad2] entry into (1+b)/2)
633  Blas::Dgemv('T', nquad1, nquad2,
634  1.0, tmp1.get()+nquad1*nquad2, nquad1,
635  base1.get()+nquad1, 1,
636  1.0, tmp2.get()+nquad2, 1);
637  }
638 
639  // Inner product with respect to the '2' direction
640  mode = mode1 = cnt = 0;
641  for(i = 0; i < order0; ++i)
642  {
643  for(j = 0; j < order1-i; ++j, ++cnt)
644  {
645  Blas::Dgemv('T', nquad2, order2-i-j,
646  1.0, base2.get()+mode*nquad2, nquad2,
647  tmp2.get()+cnt*nquad2, 1,
648  0.0, outarray.get()+mode1, 1);
649  mode += order2-i-j;
650  mode1 += order2-i-j;
651  }
652  //increment mode in case order1!=order2
653  for(j = order1-i; j < order2-i; ++j)
654  {
655  mode += order2-i-j;
656  }
657  }
658 
659  // fix for modified basis for top singular vertex component
660  // Already have evaluated (1+c)/2 (1-b)/2 (1-a)/2
661  if(m_base[0]->GetBasisType() == LibUtilities::eModified_A)
662  {
663  // add in (1+c)/2 (1+b)/2 component
664  outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
665  &tmp2[nquad2],1);
666 
667  // add in (1+c)/2 (1-b)/2 (1+a)/2 component
668  outarray[1] += Blas::Ddot(nquad2,base2.get()+nquad2,1,
669  &tmp2[nquad2*order1],1);
670  }
671  }
672 
673 
674  void StdTetExp::v_IProductWRTDerivBase(
675  const int dir,
676  const Array<OneD, const NekDouble>& inarray,
677  Array<OneD, NekDouble>& outarray)
678  {
679  StdTetExp::v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
680  }
681 
682 
683  void StdTetExp::v_IProductWRTDerivBase_MatOp(
684  const int dir,
685  const Array<OneD, const NekDouble>& inarray,
686  Array<OneD, NekDouble>& outarray)
687  {
688  ASSERTL0((dir==0)||(dir==1)||(dir==2),"input dir is out of range");
689 
690  int nq = GetTotPoints();
691  MatrixType mtype = eIProductWRTDerivBase0;
692 
693  switch (dir)
694  {
695  case 0:
696  mtype = eIProductWRTDerivBase0;
697  break;
698  case 1:
699  mtype = eIProductWRTDerivBase1;
700  break;
701  case 2:
702  mtype = eIProductWRTDerivBase2;
703  break;
704  }
705 
706  StdMatrixKey iprodmatkey(mtype,DetShapeType(),*this);
707  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
708 
709  Blas::Dgemv('N',m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
710  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
711  }
712 
713 
714  /**
715  * @param inarray Function evaluated at physical collocation
716  * points.
717  * @param outarray Inner product with respect to each basis
718  * function over the element.
719  */
720  void StdTetExp::v_IProductWRTDerivBase_SumFac(
721  const int dir,
722  const Array<OneD, const NekDouble>& inarray,
723  Array<OneD, NekDouble>& outarray)
724  {
725  int i;
726  int nquad0 = m_base[0]->GetNumPoints();
727  int nquad1 = m_base[1]->GetNumPoints();
728  int nquad2 = m_base[2]->GetNumPoints();
729  int nqtot = nquad0*nquad1*nquad2;
730  int nmodes0 = m_base[0]->GetNumModes();
731  int nmodes1 = m_base[1]->GetNumModes();
732  int wspsize = nquad0 + nquad1 + nquad2 + max(nqtot,m_ncoeffs)
733  + nquad1*nquad2*nmodes0 + nquad2*nmodes0*(2*nmodes1-nmodes0+1)/2;
734 
735  Array<OneD, NekDouble> gfac0(wspsize);
736  Array<OneD, NekDouble> gfac1(gfac0 + nquad0);
737  Array<OneD, NekDouble> gfac2(gfac1 + nquad1);
738  Array<OneD, NekDouble> tmp0 (gfac2 + nquad2);
739  Array<OneD, NekDouble> wsp(tmp0 + max(nqtot,m_ncoeffs));
740 
741  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
742  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
743  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
744 
745  // set up geometric factor: (1+z0)/2
746  for(i = 0; i < nquad0; ++i)
747  {
748  gfac0[i] = 0.5*(1+z0[i]);
749  }
750 
751  // set up geometric factor: 2/(1-z1)
752  for(i = 0; i < nquad1; ++i)
753  {
754  gfac1[i] = 2.0/(1-z1[i]);
755  }
756 
757  // Set up geometric factor: 2/(1-z2)
758  for(i = 0; i < nquad2; ++i)
759  {
760  gfac2[i] = 2.0/(1-z2[i]);
761  }
762 
763  // Derivative in first direction is always scaled as follows
764  for(i = 0; i < nquad1*nquad2; ++i)
765  {
766  Vmath::Smul(nquad0,gfac1[i%nquad1],&inarray[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
767  }
768  for(i = 0; i < nquad2; ++i)
769  {
770  Vmath::Smul(nquad0*nquad1,gfac2[i],&tmp0[0]+i*nquad0*nquad1,1,&tmp0[0]+i*nquad0*nquad1,1);
771  }
772 
773  MultiplyByQuadratureMetric(tmp0,tmp0);
774 
775  switch(dir)
776  {
777  case 0:
778  {
779  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
780  m_base[1]->GetBdata(),
781  m_base[2]->GetBdata(),
782  tmp0,outarray,wsp,
783  false, true, true);
784  }
785  break;
786  case 1:
787  {
788  Array<OneD, NekDouble> tmp3(m_ncoeffs);
789 
790  for(i = 0; i < nquad1*nquad2; ++i)
791  {
792  Vmath::Vmul(nquad0,&gfac0[0],1,&tmp0[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
793  }
794 
795  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
796  m_base[1]->GetBdata(),
797  m_base[2]->GetBdata(),
798  tmp0,tmp3,wsp,
799  false, true, true);
800 
801  for(i = 0; i < nquad2; ++i)
802  {
803  Vmath::Smul(nquad0*nquad1,gfac2[i],&inarray[0]+i*nquad0*nquad1,1,&tmp0[0]+i*nquad0*nquad1,1);
804  }
805  MultiplyByQuadratureMetric(tmp0,tmp0);
806  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
807  m_base[1]->GetDbdata(),
808  m_base[2]->GetBdata(),
809  tmp0,outarray,wsp,
810  true, false, true);
811  Vmath::Vadd(m_ncoeffs,&tmp3[0],1,&outarray[0],1,&outarray[0],1);
812  }
813  break;
814  case 2:
815  {
816  Array<OneD, NekDouble> tmp3(m_ncoeffs);
817  Array<OneD, NekDouble> tmp4(m_ncoeffs);
818  for(i = 0; i < nquad1; ++i)
819  {
820  gfac1[i] = (1+z1[i])/2;
821  }
822 
823  for(i = 0; i < nquad1*nquad2; ++i)
824  {
825  Vmath::Vmul(nquad0,&gfac0[0],1,&tmp0[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
826  }
827  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
828  m_base[1]->GetBdata(),
829  m_base[2]->GetBdata(),
830  tmp0,tmp3,wsp,
831  false, true, true);
832 
833  for(i = 0; i < nquad2; ++i)
834  {
835  Vmath::Smul(nquad0*nquad1,gfac2[i],&inarray[0]+i*nquad0*nquad1,1,&tmp0[0]+i*nquad0*nquad1,1);
836  }
837  for(i = 0; i < nquad1*nquad2; ++i)
838  {
839  Vmath::Smul(nquad0,gfac1[i%nquad1],&tmp0[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
840  }
841  MultiplyByQuadratureMetric(tmp0,tmp0);
842  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
843  m_base[1]->GetDbdata(),
844  m_base[2]->GetBdata(),
845  tmp0,tmp4,wsp,
846  true, false, true);
847 
848  MultiplyByQuadratureMetric(inarray,tmp0);
849  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
850  m_base[1]->GetBdata(),
851  m_base[2]->GetDbdata(),
852  tmp0,outarray,wsp,
853  true, true, false);
854 
855  Vmath::Vadd(m_ncoeffs,&tmp3[0],1,&outarray[0],1,&outarray[0],1);
856  Vmath::Vadd(m_ncoeffs,&tmp4[0],1,&outarray[0],1,&outarray[0],1);
857  }
858  break;
859  default:
860  {
861  ASSERTL1(false, "input dir is out of range");
862  }
863  break;
864  }
865  }
866 
867 
868  //---------------------------------------
869  // Evaluation functions
870  //---------------------------------------
871 
872 
873  void StdTetExp::v_LocCoordToLocCollapsed(
874  const Array<OneD, const NekDouble>& xi,
875  Array<OneD, NekDouble>& eta)
876  {
877  NekDouble d2 = 1.0-xi[2];
878  NekDouble d12 = -xi[1]-xi[2];
879  if(fabs(d2) < NekConstants::kNekZeroTol)
880  {
881  if(d2>=0.)
882  {
883  d2 = NekConstants::kNekZeroTol;
884  }
885  else
886  {
887  d2 = -NekConstants::kNekZeroTol;
888  }
889  }
890  if(fabs(d12) < NekConstants::kNekZeroTol)
891  {
892  if(d12>=0.)
893  {
894  d12 = NekConstants::kNekZeroTol;
895  }
896  else
897  {
898  d12 = -NekConstants::kNekZeroTol;
899  }
900  }
901  eta[0] = 2.0*(1.0+xi[0])/d12 - 1.0;
902  eta[1] = 2.0*(1.0+xi[1])/d2 - 1.0;
903  eta[2] = xi[2];
904  }
905 
906  void StdTetExp::v_LocCollapsedToLocCoord(
907  const Array<OneD, const NekDouble>& eta,
908  Array<OneD, NekDouble>& xi)
909  {
910  xi[1] = (1.0 + eta[0]) * ( 1.0 - eta[2]) * 0.5 - 1.0;
911  xi[0] = (1.0 + eta[0]) * (-xi[1] - eta[2]) * 0.5 - 1.0;
912  xi[2] = eta[2];
913  }
914 
915  void StdTetExp::v_FillMode(
916  const int mode,
917  Array<OneD, NekDouble> &outarray)
918  {
919  Array<OneD, NekDouble> tmp(m_ncoeffs,0.0);
920  tmp[mode] = 1.0;
921  StdTetExp::v_BwdTrans(tmp, outarray);
922  }
923 
924  NekDouble StdTetExp::v_PhysEvaluateBasis(
925  const Array<OneD, const NekDouble>& coords,
926  int mode)
927  {
928  Array<OneD, NekDouble> coll(3);
929  LocCoordToLocCollapsed(coords, coll);
930 
931  const int nm1 = m_base[1]->GetNumModes();
932  const int nm2 = m_base[2]->GetNumModes();
933 
934  const int b = 2 * nm2 + 1;
935  const int mode0 = floor(0.5 * (b - sqrt(b * b - 8.0 * mode / nm1)));
936  const int tmp =
937  mode - nm1*(mode0 * (nm2-1) + 1 - (mode0 - 2)*(mode0 - 1) / 2);
938  const int mode1 = tmp / (nm2 - mode0);
939  const int mode2 = tmp % (nm2 - mode0);
940 
941  if (m_base[0]->GetBasisType() == LibUtilities::eModified_A)
942  {
943  // Handle the collapsed vertices and edges in the modified
944  // basis.
945  if (mode == 1)
946  {
947  // Collapsed top vertex
948  return StdExpansion::BaryEvaluateBasis<2>(coll[2], 1);
949  }
950  else if (mode0 == 0 && mode2 == 1)
951  {
952  return
953  StdExpansion::BaryEvaluateBasis<1>(coll[1], 0) *
954  StdExpansion::BaryEvaluateBasis<2>(coll[2], 1);
955  }
956  else if (mode0 == 1 && mode1 == 1 && mode2 == 0)
957  {
958  return
959  StdExpansion::BaryEvaluateBasis<0>(coll[0], 0) *
960  StdExpansion::BaryEvaluateBasis<1>(coll[1], 1);
961  }
962  }
963 
964  return
965  StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
966  StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
967  StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
968  }
969 
970  void StdTetExp::v_GetTraceNumModes(
971  const int fid,
972 
973  int &numModes0,
974  int &numModes1,
975  Orientation faceOrient)
976  {
977  boost::ignore_unused(faceOrient);
978 
979  int nummodes [3] = {m_base[0]->GetNumModes(),
980  m_base[1]->GetNumModes(),
981  m_base[2]->GetNumModes()};
982  switch(fid)
983  {
984  case 0:
985  {
986  numModes0 = nummodes[0];
987  numModes1 = nummodes[1];
988  }
989  break;
990  case 1:
991  {
992  numModes0 = nummodes[0];
993  numModes1 = nummodes[2];
994  }
995  break;
996  case 2:
997  case 3:
998  {
999  numModes0 = nummodes[1];
1000  numModes1 = nummodes[2];
1001  }
1002  break;
1003  }
1004  }
1005 
1006  //---------------------------
1007  // Helper functions
1008  //---------------------------
1009 
1010  int StdTetExp::v_GetNverts() const
1011  {
1012  return 4;
1013  }
1014 
1015  int StdTetExp::v_GetNedges() const
1016  {
1017  return 6;
1018  }
1019 
1020  int StdTetExp::v_GetNtraces() const
1021  {
1022  return 4;
1023  }
1024 
1025  LibUtilities::ShapeType StdTetExp::v_DetShapeType() const
1026  {
1027  return DetShapeType();
1028  }
1029 
1030  int StdTetExp::v_NumBndryCoeffs() const
1031  {
1032  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1033  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1034  "BasisType is not a boundary interior form");
1035  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
1036  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1037  "BasisType is not a boundary interior form");
1038  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1039  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1040  "BasisType is not a boundary interior form");
1041 
1042  int P = m_base[0]->GetNumModes();
1043  int Q = m_base[1]->GetNumModes();
1044  int R = m_base[2]->GetNumModes();
1045 
1046  return LibUtilities::StdTetData::
1047  getNumberOfBndCoefficients(P, Q, R);
1048  }
1049 
1050  int StdTetExp::v_NumDGBndryCoeffs() const
1051  {
1052  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1053  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1054  "BasisType is not a boundary interior form");
1055  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
1056  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1057  "BasisType is not a boundary interior form");
1058  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1059  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1060  "BasisType is not a boundary interior form");
1061 
1062  int P = m_base[0]->GetNumModes()-1;
1063  int Q = m_base[1]->GetNumModes()-1;
1064  int R = m_base[2]->GetNumModes()-1;
1065 
1066 
1067  return (Q+1) + P*(1 + 2*Q - P)/2 // base face
1068  + (R+1) + P*(1 + 2*R - P)/2 // front face
1069  + 2*(R+1) + Q*(1 + 2*R - Q); // back two faces
1070  }
1071 
1072  int StdTetExp::v_GetTraceNcoeffs(const int i) const
1073  {
1074  ASSERTL2((i >= 0) && (i <= 3), "face id is out of range");
1075  int nFaceCoeffs = 0;
1076  int nummodesA, nummodesB, P, Q;
1077  if (i == 0)
1078  {
1079  nummodesA = GetBasisNumModes(0);
1080  nummodesB = GetBasisNumModes(1);
1081  }
1082  else if ((i == 1) || (i == 2))
1083  {
1084  nummodesA = GetBasisNumModes(0);
1085  nummodesB = GetBasisNumModes(2);
1086  }
1087  else
1088  {
1089  nummodesA = GetBasisNumModes(1);
1090  nummodesB = GetBasisNumModes(2);
1091  }
1092  P = nummodesA - 1;
1093  Q = nummodesB - 1;
1094  nFaceCoeffs = Q+1 + (P*(1 + 2*Q - P))/2;
1095  return nFaceCoeffs;
1096  }
1097 
1098  int StdTetExp::v_GetTraceIntNcoeffs(const int i) const
1099  {
1100  ASSERTL2((i >= 0) && (i <= 3), "face id is out of range");
1101  int Pi = m_base[0]->GetNumModes() - 2;
1102  int Qi = m_base[1]->GetNumModes() - 2;
1103  int Ri = m_base[2]->GetNumModes() - 2;
1104 
1105  if((i == 0))
1106  {
1107  return Pi * (2*Qi - Pi - 1) / 2;
1108  }
1109  else if((i == 1))
1110  {
1111  return Pi * (2*Ri - Pi - 1) / 2;
1112  }
1113  else
1114  {
1115  return Qi * (2*Ri - Qi - 1) / 2;
1116  }
1117  }
1118 
1119  int StdTetExp::v_GetTotalTraceIntNcoeffs() const
1120  {
1121  int Pi = m_base[0]->GetNumModes() - 2;
1122  int Qi = m_base[1]->GetNumModes() - 2;
1123  int Ri = m_base[2]->GetNumModes() - 2;
1124 
1125  return Pi * (2*Qi - Pi - 1) / 2 +
1126  Pi * (2*Ri - Pi - 1) / 2 +
1127  Qi * (2*Ri - Qi - 1);
1128  }
1129 
1130  int StdTetExp::v_GetTraceNumPoints(const int i) const
1131  {
1132  ASSERTL2(i >= 0 && i <= 3, "face id is out of range");
1133 
1134  if (i == 0)
1135  {
1136  return m_base[0]->GetNumPoints()*
1137  m_base[1]->GetNumPoints();
1138  }
1139  else if (i == 1)
1140  {
1141  return m_base[0]->GetNumPoints()*
1142  m_base[2]->GetNumPoints();
1143  }
1144  else
1145  {
1146  return m_base[1]->GetNumPoints()*
1147  m_base[2]->GetNumPoints();
1148  }
1149  }
1150 
1151 
1152  int StdTetExp::v_GetEdgeNcoeffs(const int i) const
1153  {
1154  ASSERTL2((i >= 0) && (i <= 5), "edge id is out of range");
1155  int P = m_base[0]->GetNumModes();
1156  int Q = m_base[1]->GetNumModes();
1157  int R = m_base[2]->GetNumModes();
1158 
1159  if (i == 0)
1160  {
1161  return P;
1162  }
1163  else if (i == 1 || i == 2)
1164  {
1165  return Q;
1166  }
1167  else
1168  {
1169  return R;
1170  }
1171  }
1172 
1173  LibUtilities::PointsKey StdTetExp::v_GetTracePointsKey(
1174  const int i, const int j) const
1175  {
1176  ASSERTL2(i >= 0 && i <= 3, "face id is out of range");
1177  ASSERTL2(j == 0 || j == 1, "face direction is out of range");
1178 
1179  if (i == 0)
1180  {
1181  return m_base[j]->GetPointsKey();
1182  }
1183  else if (i == 1)
1184  {
1185  return m_base[2*j]->GetPointsKey();
1186  }
1187  else
1188  {
1189  return m_base[j+1]->GetPointsKey();
1190  }
1191  }
1192 
1193  int StdTetExp::v_CalcNumberOfCoefficients(
1194  const std::vector<unsigned int>& nummodes,
1195  int & modes_offset)
1196  {
1197  int nmodes = LibUtilities::StdTetData::getNumberOfCoefficients(
1198  nummodes[modes_offset],
1199  nummodes[modes_offset+1],
1200  nummodes[modes_offset+2]);
1201  modes_offset += 3;
1202 
1203  return nmodes;
1204  }
1205 
1206  const LibUtilities::BasisKey StdTetExp::v_GetTraceBasisKey(
1207  const int i, const int k) const
1208  {
1209  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
1210  ASSERTL2(k == 0 || k == 1, "face direction out of range");
1211 
1212  int dir = k;
1213  switch(i)
1214  {
1215  case 0:
1216  dir = k;
1217  break;
1218  case 1:
1219  dir = 2*k;
1220  break;
1221  case 2:
1222  case 3:
1223  dir = k+1;
1224  break;
1225  }
1226 
1227  return EvaluateTriFaceBasisKey(k,
1228  m_base[dir]->GetBasisType(),
1229  m_base[dir]->GetNumPoints(),
1230  m_base[dir]->GetNumModes());
1231 
1232  // Should not get here.
1233  return LibUtilities::NullBasisKey;
1234  }
1235 
1236 
1237  void StdTetExp::v_GetCoords(
1238  Array<OneD, NekDouble> &xi_x,
1239  Array<OneD, NekDouble> &xi_y,
1240  Array<OneD, NekDouble> &xi_z)
1241  {
1242  Array<OneD, const NekDouble> eta_x = m_base[0]->GetZ();
1243  Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
1244  Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
1245  int Qx = GetNumPoints(0);
1246  int Qy = GetNumPoints(1);
1247  int Qz = GetNumPoints(2);
1248 
1249  // Convert collapsed coordinates into cartesian coordinates: eta
1250  // --> xi
1251  for( int k = 0; k < Qz; ++k ) {
1252  for( int j = 0; j < Qy; ++j ) {
1253  for( int i = 0; i < Qx; ++i ) {
1254  int s = i + Qx*(j + Qy*k);
1255  xi_x[s] = (eta_x[i] + 1.0) * (1.0 - eta_y[j]) * (1.0 - eta_z[k]) / 4 - 1.0;
1256  xi_y[s] = (eta_y[j] + 1.0) * (1.0 - eta_z[k]) / 2 - 1.0;
1257  xi_z[s] = eta_z[k];
1258  }
1259  }
1260  }
1261  }
1262 
1263  bool StdTetExp::v_IsBoundaryInteriorExpansion()
1264  {
1265  return (m_base[0]->GetBasisType() == LibUtilities::eModified_A) &&
1266  (m_base[1]->GetBasisType() == LibUtilities::eModified_B) &&
1267  (m_base[2]->GetBasisType() == LibUtilities::eModified_C);
1268  }
1269 
1270 
1271  //--------------------------
1272  // Mappings
1273  //--------------------------
1274  int StdTetExp::v_GetVertexMap(const int localVertexId, bool useCoeffPacking)
1275  {
1276  ASSERTL0((GetBasisType(0)==LibUtilities::eModified_A)||
1277  (GetBasisType(1)==LibUtilities::eModified_B)||
1278  (GetBasisType(2)==LibUtilities::eModified_C),
1279  "Mapping not defined for this type of basis");
1280 
1281  int localDOF = 0;
1282  if(useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1283  {
1284  switch(localVertexId)
1285  {
1286  case 0:
1287  {
1288  localDOF = GetMode(0,0,0);
1289  break;
1290  }
1291  case 1:
1292  {
1293  localDOF = GetMode(0,0,1);
1294  break;
1295  }
1296  case 2:
1297  {
1298  localDOF = GetMode(0,1,0);
1299  break;
1300  }
1301  case 3:
1302  {
1303  localDOF = GetMode(1,0,0);
1304  break;
1305  }
1306  default:
1307  {
1308  ASSERTL0(false,"Vertex ID must be between 0 and 3");
1309  break;
1310  }
1311  }
1312  }
1313  else
1314  {
1315  switch(localVertexId)
1316  {
1317  case 0:
1318  {
1319  localDOF = GetMode(0,0,0);
1320  break;
1321  }
1322  case 1:
1323  {
1324  localDOF = GetMode(1,0,0);
1325  break;
1326  }
1327  case 2:
1328  {
1329  localDOF = GetMode(0,1,0);
1330  break;
1331  }
1332  case 3:
1333  {
1334  localDOF = GetMode(0,0,1);
1335  break;
1336  }
1337  default:
1338  {
1339  ASSERTL0(false,"Vertex ID must be between 0 and 3");
1340  break;
1341  }
1342  }
1343 
1344  }
1345 
1346  return localDOF;
1347  }
1348 
1349  /**
1350  * Maps interior modes of an edge to the elemental modes.
1351  */
1352 
1353  /**
1354  * List of all interior modes in the expansion.
1355  */
1356  void StdTetExp::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
1357  {
1358  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1359  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1360  "BasisType is not a boundary interior form");
1361  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
1362  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1363  "BasisType is not a boundary interior form");
1364  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1365  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1366  "BasisType is not a boundary interior form");
1367 
1368  int P = m_base[0]->GetNumModes();
1369  int Q = m_base[1]->GetNumModes();
1370  int R = m_base[2]->GetNumModes();
1371 
1372  int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1373 
1374  if(outarray.size() != nIntCoeffs)
1375  {
1376  outarray = Array<OneD, unsigned int>(nIntCoeffs);
1377  }
1378 
1379  int idx = 0;
1380  for (int i = 2; i < P-2; ++i)
1381  {
1382  for (int j = 1; j < Q-i-1; ++j)
1383  {
1384  for (int k = 1; k < R-i-j; ++k)
1385  {
1386  outarray[idx++] = GetMode(i,j,k);
1387  }
1388  }
1389  }
1390  }
1391 
1392  /**
1393  * List of all boundary modes in the the expansion.
1394  */
1395  void StdTetExp::v_GetBoundaryMap(Array<OneD, unsigned int>& outarray)
1396  {
1397  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1398  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1399  "BasisType is not a boundary interior form");
1400  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
1401  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1402  "BasisType is not a boundary interior form");
1403  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1404  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1405  "BasisType is not a boundary interior form");
1406 
1407  int P = m_base[0]->GetNumModes();
1408  int Q = m_base[1]->GetNumModes();
1409  int R = m_base[2]->GetNumModes();
1410 
1411  int i,j,k;
1412  int idx = 0;
1413 
1414  int nBnd = NumBndryCoeffs();
1415 
1416  if (outarray.size() != nBnd)
1417  {
1418  outarray = Array<OneD, unsigned int>(nBnd);
1419  }
1420 
1421  for (i = 0; i < P; ++i)
1422  {
1423  // First two Q-R planes are entirely boundary modes
1424  if (i < 2)
1425  {
1426  for (j = 0; j < Q-i; j++)
1427  {
1428  for (k = 0; k < R-i-j; ++k)
1429  {
1430  outarray[idx++] = GetMode(i,j,k);
1431  }
1432  }
1433  }
1434  // Remaining Q-R planes contain boundary modes on bottom and
1435  // left edge.
1436  else
1437  {
1438  for (k = 0; k < R-i; ++k)
1439  {
1440  outarray[idx++] = GetMode(i,0,k);
1441  }
1442  for (j = 1; j < Q-i; ++j)
1443  {
1444  outarray[idx++] = GetMode(i,j,0);
1445  }
1446  }
1447  }
1448  }
1449 
1450  /**
1451  * Maps Expansion2D modes of a 2D face to the corresponding expansion
1452  * modes.
1453  */
1454  void StdTetExp::v_GetTraceToElementMap(
1455  const int fid,
1456  Array<OneD, unsigned int> &maparray,
1457  Array<OneD, int> &signarray,
1458  const Orientation faceOrient,
1459  int P,
1460  int Q)
1461  {
1462  int nummodesA=0,nummodesB=0, i, j, k, idx;
1463 
1464  ASSERTL1(v_IsBoundaryInteriorExpansion(),
1465  "Method only implemented for Modified_A BasisType (x "
1466  "direction), Modified_B BasisType (y direction), and "
1467  "Modified_C BasisType(z direction)");
1468 
1469  int nFaceCoeffs = 0;
1470 
1471  switch(fid)
1472  {
1473  case 0:
1474  nummodesA = m_base[0]->GetNumModes();
1475  nummodesB = m_base[1]->GetNumModes();
1476  break;
1477  case 1:
1478  nummodesA = m_base[0]->GetNumModes();
1479  nummodesB = m_base[2]->GetNumModes();
1480  break;
1481  case 2:
1482  case 3:
1483  nummodesA = m_base[1]->GetNumModes();
1484  nummodesB = m_base[2]->GetNumModes();
1485  break;
1486  default:
1487  ASSERTL0(false,"fid must be between 0 and 3");
1488  }
1489 
1490  bool CheckForZeroedModes = false;
1491  if(P == -1)
1492  {
1493  P = nummodesA;
1494  Q = nummodesB;
1495  }
1496  else
1497  {
1498  CheckForZeroedModes = true;
1499  }
1500 
1501  nFaceCoeffs = P*(2*Q-P+1)/2;
1502 
1503  // Allocate the map array and sign array; set sign array to ones (+)
1504  if(maparray.size() != nFaceCoeffs)
1505  {
1506  maparray = Array<OneD, unsigned int>(nFaceCoeffs,1);
1507  }
1508 
1509  if(signarray.size() != nFaceCoeffs)
1510  {
1511  signarray = Array<OneD, int>(nFaceCoeffs,1);
1512  }
1513  else
1514  {
1515  fill(signarray.get(),signarray.get()+nFaceCoeffs, 1 );
1516  }
1517 
1518  switch (fid)
1519  {
1520  case 0:
1521  idx = 0;
1522  for (i = 0; i < P; ++i)
1523  {
1524  for (j = 0; j < Q-i; ++j)
1525  {
1526  if ((int)faceOrient == 7 && i > 1)
1527  {
1528  signarray[idx] = (i%2 ? -1 : 1);
1529  }
1530  maparray[idx++] = GetMode(i,j,0);
1531  }
1532  }
1533  break;
1534  case 1:
1535  idx = 0;
1536  for (i = 0; i < P; ++i)
1537  {
1538  for (k = 0; k < Q-i; ++k)
1539  {
1540  if ((int)faceOrient == 7 && i > 1)
1541  {
1542  signarray[idx] = (i%2 ? -1: 1);
1543  }
1544  maparray[idx++] = GetMode(i,0,k);
1545  }
1546  }
1547  break;
1548  case 2:
1549  idx = 0;
1550  for (j = 0; j < P-1; ++j)
1551  {
1552  for (k = 0; k < Q-1-j; ++k)
1553  {
1554  if ((int)faceOrient == 7 && j > 1)
1555  {
1556  signarray[idx] = ((j+1)%2 ? -1: 1);
1557  }
1558  maparray[idx++] = GetMode(1,j,k);
1559  // Incorporate modes from zeroth plane where needed.
1560  if (j == 0 && k == 0)
1561  {
1562  maparray[idx++] = GetMode(0,0,1);
1563  }
1564  if (j == 0 && k == Q-2)
1565  {
1566  for (int r = 0; r < Q-1; ++r)
1567  {
1568  maparray[idx++] = GetMode(0,1,r);
1569  }
1570  }
1571  }
1572  }
1573  break;
1574  case 3:
1575  idx = 0;
1576  for (j = 0; j < P; ++j)
1577  {
1578  for (k = 0; k < Q-j; ++k)
1579  {
1580  if ((int)faceOrient == 7 && j > 1)
1581  {
1582  signarray[idx] = (j%2 ? -1: 1);
1583  }
1584  maparray[idx++] = GetMode(0,j,k);
1585  }
1586  }
1587  break;
1588  default:
1589  ASSERTL0(false, "Element map not available.");
1590  }
1591 
1592  if ((int)faceOrient == 7)
1593  {
1594  swap(maparray[0], maparray[Q]);
1595 
1596  for (i = 1; i < Q-1; ++i)
1597  {
1598  swap(maparray[i+1], maparray[Q+i]);
1599  }
1600  }
1601 
1602  if(CheckForZeroedModes)
1603  {
1604  // zero signmap and set maparray to zero if elemental
1605  // modes are not as large as face modesl
1606  idx = 0;
1607  for (j = 0; j < nummodesA; ++j)
1608  {
1609  idx += nummodesB-j;
1610  for (k = nummodesB-j; k < Q-j; ++k)
1611  {
1612  signarray[idx] = 0.0;
1613  maparray[idx++] = maparray[0];
1614  }
1615  }
1616 
1617  for (j = nummodesA; j < P; ++j)
1618  {
1619  for (k = 0; k < Q-j; ++k)
1620  {
1621  signarray[idx] = 0.0;
1622  maparray[idx++] = maparray[0];
1623  }
1624  }
1625  }
1626  }
1627 
1628  /**
1629  * Maps interior modes of an edge to the elemental modes.
1630  */
1631  void StdTetExp::v_GetEdgeInteriorToElementMap(
1632  const int eid,
1633  Array<OneD, unsigned int> &maparray,
1634  Array<OneD, int> &signarray,
1635  const Orientation edgeOrient)
1636  {
1637  int i;
1638  const int P = m_base[0]->GetNumModes();
1639  const int Q = m_base[1]->GetNumModes();
1640  const int R = m_base[2]->GetNumModes();
1641 
1642  const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid)-2;
1643 
1644  if(maparray.size() != nEdgeIntCoeffs)
1645  {
1646  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1647  }
1648  else
1649  {
1650  fill( maparray.get(), maparray.get() + nEdgeIntCoeffs, 0);
1651  }
1652 
1653  if(signarray.size() != nEdgeIntCoeffs)
1654  {
1655  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1656  }
1657  else
1658  {
1659  fill( signarray.get() , signarray.get()+nEdgeIntCoeffs, 1 );
1660  }
1661 
1662  switch (eid)
1663  {
1664  case 0:
1665  for (i = 0; i < P-2; ++i)
1666  {
1667  maparray[i] = GetMode(i+2, 0, 0);
1668  }
1669  if(edgeOrient==eBackwards)
1670  {
1671  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1672  {
1673  signarray[i] = -1;
1674  }
1675  }
1676  break;
1677  case 1:
1678  for (i = 0; i < Q-2; ++i)
1679  {
1680  maparray[i] = GetMode(1, i+1, 0);
1681  }
1682  if(edgeOrient==eBackwards)
1683  {
1684  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1685  {
1686  signarray[i] = -1;
1687  }
1688  }
1689  break;
1690  case 2:
1691  for (i = 0; i < Q-2; ++i)
1692  {
1693  maparray[i] = GetMode(0, i+2, 0);
1694  }
1695  if(edgeOrient==eBackwards)
1696  {
1697  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1698  {
1699  signarray[i] = -1;
1700  }
1701  }
1702  break;
1703  case 3:
1704  for (i = 0; i < R-2; ++i)
1705  {
1706  maparray[i] = GetMode(0, 0, i+2);
1707  }
1708  if(edgeOrient==eBackwards)
1709  {
1710  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1711  {
1712  signarray[i] = -1;
1713  }
1714  }
1715  break;
1716  case 4:
1717  for (i = 0; i < R - 2; ++i)
1718  {
1719  maparray[i] = GetMode(1, 0, i+1);
1720  }
1721  if(edgeOrient==eBackwards)
1722  {
1723  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1724  {
1725  signarray[i] = -1;
1726  }
1727  }
1728  break;
1729  case 5:
1730  for (i = 0; i < R-2; ++i)
1731  {
1732  maparray[i] = GetMode(0, 1, i+1);
1733  }
1734  if(edgeOrient==eBackwards)
1735  {
1736  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1737  {
1738  signarray[i] = -1;
1739  }
1740  }
1741  break;
1742  default:
1743  ASSERTL0(false, "Edge not defined.");
1744  break;
1745  }
1746  }
1747 
1748  void StdTetExp::v_GetTraceInteriorToElementMap(
1749  const int fid,
1750  Array<OneD, unsigned int> &maparray,
1751  Array<OneD, int> &signarray,
1752  const Orientation faceOrient)
1753  {
1754  int i, j, idx, k;
1755  const int P = m_base[0]->GetNumModes();
1756  const int Q = m_base[1]->GetNumModes();
1757  const int R = m_base[2]->GetNumModes();
1758 
1759  const int nFaceIntCoeffs = v_GetTraceIntNcoeffs(fid);
1760 
1761  if(maparray.size() != nFaceIntCoeffs)
1762  {
1763  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1764  }
1765 
1766  if(signarray.size() != nFaceIntCoeffs)
1767  {
1768  signarray = Array<OneD, int>(nFaceIntCoeffs,1);
1769  }
1770  else
1771  {
1772  fill( signarray.get() , signarray.get()+nFaceIntCoeffs, 1 );
1773  }
1774 
1775  switch (fid)
1776  {
1777  case 0:
1778  idx = 0;
1779  for (i = 2; i < P-1; ++i)
1780  {
1781  for (j = 1; j < Q-i; ++j)
1782  {
1783  if ((int)faceOrient == 7)
1784  {
1785  signarray[idx] = (i%2 ? -1 : 1);
1786  }
1787  maparray[idx++] = GetMode(i,j,0);
1788  }
1789  }
1790  break;
1791  case 1:
1792  idx = 0;
1793  for (i = 2; i < P; ++i)
1794  {
1795  for (k = 1; k < R-i; ++k)
1796  {
1797  if ((int)faceOrient == 7)
1798  {
1799  signarray[idx] = (i%2 ? -1: 1);
1800  }
1801  maparray[idx++] = GetMode(i,0,k);
1802  }
1803  }
1804  break;
1805  case 2:
1806  idx = 0;
1807  for (j = 1; j < Q-2; ++j)
1808  {
1809  for (k = 1; k < R-1-j; ++k)
1810  {
1811  if ((int)faceOrient == 7)
1812  {
1813  signarray[idx] = ((j+1)%2 ? -1: 1);
1814  }
1815  maparray[idx++] = GetMode(1,j,k);
1816  }
1817  }
1818  break;
1819  case 3:
1820  idx = 0;
1821  for (j = 2; j < Q-1; ++j)
1822  {
1823  for (k = 1; k < R-j; ++k)
1824  {
1825  if ((int)faceOrient == 7)
1826  {
1827  signarray[idx] = (j%2 ? -1: 1);
1828  }
1829  maparray[idx++] = GetMode(0,j,k);
1830  }
1831  }
1832  break;
1833  default:
1834  ASSERTL0(false, "Face interior map not available.");
1835  break;
1836  }
1837  }
1838  //---------------------------------------
1839  // Wrapper functions
1840  //---------------------------------------
1841  DNekMatSharedPtr StdTetExp::v_GenMatrix(const StdMatrixKey &mkey)
1842  {
1843 
1844  MatrixType mtype = mkey.GetMatrixType();
1845 
1846  DNekMatSharedPtr Mat;
1847 
1848  switch(mtype)
1849  {
1850  case ePhysInterpToEquiSpaced:
1851  {
1852  int nq0 = m_base[0]->GetNumPoints();
1853  int nq1 = m_base[1]->GetNumPoints();
1854  int nq2 = m_base[2]->GetNumPoints();
1855  int nq;
1856 
1857  // take definition from key
1858  if(mkey.ConstFactorExists(eFactorConst))
1859  {
1860  nq = (int) mkey.GetConstFactor(eFactorConst);
1861  }
1862  else
1863  {
1864  nq = max(nq0,max(nq1,nq2));
1865  }
1866 
1867  int neq = LibUtilities::StdTetData::
1868  getNumberOfCoefficients(nq,nq,nq);
1869  Array<OneD, Array<OneD, NekDouble> > coords(neq);
1870  Array<OneD, NekDouble> coll(3);
1871  Array<OneD, DNekMatSharedPtr> I(3);
1872  Array<OneD, NekDouble> tmp(nq0);
1873 
1874  Mat = MemoryManager<DNekMat>::
1875  AllocateSharedPtr(neq, nq0 * nq1 * nq2);
1876  int cnt = 0;
1877 
1878  for(int i = 0; i < nq; ++i)
1879  {
1880  for(int j = 0; j < nq-i; ++j)
1881  {
1882  for(int k = 0; k < nq-i-j; ++k,++cnt)
1883  {
1884  coords[cnt] = Array<OneD, NekDouble>(3);
1885  coords[cnt][0] = -1.0 + 2*k/(NekDouble)(nq-1);
1886  coords[cnt][1] = -1.0 + 2*j/(NekDouble)(nq-1);
1887  coords[cnt][2] = -1.0 + 2*i/(NekDouble)(nq-1);
1888  }
1889  }
1890  }
1891 
1892  for(int i = 0; i < neq; ++i)
1893  {
1894  LocCoordToLocCollapsed(coords[i],coll);
1895 
1896  I[0] = m_base[0]->GetI(coll);
1897  I[1] = m_base[1]->GetI(coll+1);
1898  I[2] = m_base[2]->GetI(coll+2);
1899 
1900  // interpolate first coordinate direction
1901  NekDouble fac;
1902  for( int k = 0; k < nq2; ++k)
1903  {
1904  for (int j = 0; j < nq1; ++j)
1905  {
1906 
1907  fac = (I[1]->GetPtr())[j]*(I[2]->GetPtr())[k];
1908  Vmath::Smul(nq0,fac,I[0]->GetPtr(),1,tmp,1);
1909 
1910  Vmath::Vcopy(nq0, &tmp[0], 1,
1911  Mat->GetRawPtr()+
1912  k*nq0*nq1*neq+
1913  j*nq0*neq+i,neq);
1914  }
1915  }
1916  }
1917  }
1918  break;
1919  default:
1920  {
1921  Mat = StdExpansion::CreateGeneralMatrix(mkey);
1922  }
1923  break;
1924  }
1925 
1926  return Mat;
1927  }
1928 
1929  DNekMatSharedPtr StdTetExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
1930  {
1931  return v_GenMatrix(mkey);
1932  }
1933 
1934 
1935  //---------------------------------------
1936  // Private helper functions
1937  //---------------------------------------
1938 
1939  /**
1940  * @brief Compute the mode number in the expansion for a particular
1941  * tensorial combination.
1942  *
1943  * Modes are numbered with the r index travelling fastest, followed by
1944  * q and then p, and each q-r plane is of size
1945  * (Q+1)*(Q+2)/2+max(0,R-Q-p)*Q. For example, when P=2, Q=3 and R=4
1946  * the indexing inside each q-r plane (with r increasing upwards and q
1947  * to the right) is:
1948  *
1949  * p = 0: p = 1: p = 2:
1950  * ----------------------------------
1951  * 4
1952  * 3 8 17
1953  * 2 7 11 16 20 26
1954  * 1 6 10 13 15 19 22 25 28
1955  * 0 5 9 12 14 18 21 23 24 27 29
1956  *
1957  * Note that in this element, we must have that \f$ P \leq Q \leq
1958  * R\f$.
1959  */
1960  int StdTetExp::GetMode(const int I, const int J, const int K)
1961  {
1962  const int Q = m_base[1]->GetNumModes();
1963  const int R = m_base[2]->GetNumModes();
1964 
1965  int i,j,q_hat,k_hat;
1966  int cnt = 0;
1967 
1968  // Traverse to q-r plane number I
1969  for (i = 0; i < I; ++i)
1970  {
1971  // Size of triangle part
1972  q_hat = min(Q,R-i);
1973  // Size of rectangle part
1974  k_hat = max(R-Q-i,0);
1975  cnt += q_hat*(q_hat+1)/2 + k_hat*Q;
1976  }
1977 
1978  // Traverse to q column J
1979  q_hat = R-I;
1980  for (j = 0; j < J; ++j)
1981  {
1982  cnt += q_hat;
1983  q_hat--;
1984  }
1985 
1986  // Traverse up stacks to K
1987  cnt += K;
1988 
1989  return cnt;
1990  }
1991 
1992  void StdTetExp::v_MultiplyByStdQuadratureMetric(
1993  const Array<OneD, const NekDouble>& inarray,
1994  Array<OneD, NekDouble>& outarray)
1995  {
1996  int i, j;
1997 
1998  int nquad0 = m_base[0]->GetNumPoints();
1999  int nquad1 = m_base[1]->GetNumPoints();
2000  int nquad2 = m_base[2]->GetNumPoints();
2001 
2002  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
2003  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
2004  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
2005 
2006  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
2007  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
2008 
2009  // multiply by integration constants
2010  for(i = 0; i < nquad1*nquad2; ++i)
2011  {
2012  Vmath::Vmul(nquad0,(NekDouble*)&inarray[0]+i*nquad0,1,
2013  w0.get(),1, &outarray[0]+i*nquad0,1);
2014  }
2015 
2016  switch(m_base[1]->GetPointsType())
2017  {
2018  // (1,0) Jacobi Inner product.
2019  case LibUtilities::eGaussRadauMAlpha1Beta0:
2020  for(j = 0; j < nquad2; ++j)
2021  {
2022  for(i = 0; i < nquad1; ++i)
2023  {
2024  Blas::Dscal(nquad0,0.5*w1[i], &outarray[0]+i*nquad0+
2025  j*nquad0*nquad1,1);
2026  }
2027  }
2028  break;
2029 
2030  default:
2031  for(j = 0; j < nquad2; ++j)
2032  {
2033  for(i = 0; i < nquad1; ++i)
2034  {
2035  Blas::Dscal(nquad0,
2036  0.5*(1-z1[i])*w1[i],
2037  &outarray[0]+i*nquad0 + j*nquad0*nquad1,
2038  1 );
2039  }
2040  }
2041  break;
2042  }
2043 
2044  switch(m_base[2]->GetPointsType())
2045  {
2046  // (2,0) Jacobi inner product.
2047  case LibUtilities::eGaussRadauMAlpha2Beta0:
2048  for(i = 0; i < nquad2; ++i)
2049  {
2050  Blas::Dscal(nquad0*nquad1, 0.25*w2[i],
2051  &outarray[0]+i*nquad0*nquad1, 1);
2052  }
2053  break;
2054  // (1,0) Jacobi inner product.
2055  case LibUtilities::eGaussRadauMAlpha1Beta0:
2056  for(i = 0; i < nquad2; ++i)
2057  {
2058  Blas::Dscal(nquad0*nquad1, 0.25*(1-z2[i])*w2[i],
2059  &outarray[0]+i*nquad0*nquad1, 1);
2060  }
2061  break;
2062  default:
2063  for(i = 0; i < nquad2; ++i)
2064  {
2065  Blas::Dscal(nquad0*nquad1,0.25*(1-z2[i])*(1-z2[i])*w2[i],
2066  &outarray[0]+i*nquad0*nquad1,1);
2067  }
2068  break;
2069  }
2070  }
2071 
2072  void StdTetExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
2073  const StdMatrixKey &mkey)
2074  {
2075  //To do : 1) add a test to ensure 0 \leq SvvCutoff \leq 1.
2076  // 2) check if the transfer function needs an analytical
2077  // Fourier transform.
2078  // 3) if it doesn't : find a transfer function that renders
2079  // the if( cutoff_a ...) useless to reduce computational
2080  // cost.
2081  // 4) add SVVDiffCoef to both models!!
2082 
2083  int qa = m_base[0]->GetNumPoints();
2084  int qb = m_base[1]->GetNumPoints();
2085  int qc = m_base[2]->GetNumPoints();
2086  int nmodes_a = m_base[0]->GetNumModes();
2087  int nmodes_b = m_base[1]->GetNumModes();
2088  int nmodes_c = m_base[2]->GetNumModes();
2089 
2090  // Declare orthogonal basis.
2091  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
2092  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
2093  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
2094 
2095  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
2096  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B,nmodes_b,pb);
2097  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_C,nmodes_c,pc);
2098 
2099  StdTetExp OrthoExp(Ba,Bb,Bc);
2100 
2101 
2102  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2103  int i,j,k,cnt = 0;
2104 
2105  // project onto physical space.
2106  OrthoExp.FwdTrans(array,orthocoeffs);
2107 
2108  if(mkey.ConstFactorExists(eFactorSVVPowerKerDiffCoeff))
2109  {
2110  // Rodrigo's power kernel
2111  NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
2112  NekDouble SvvDiffCoeff =
2113  mkey.GetConstFactor(eFactorSVVPowerKerDiffCoeff)*
2114  mkey.GetConstFactor(eFactorSVVDiffCoeff);
2115 
2116  for(int i = 0; i < nmodes_a; ++i)
2117  {
2118  for(int j = 0; j < nmodes_b-j; ++j)
2119  {
2120  NekDouble fac1 = std::max(
2121  pow((1.0*i)/(nmodes_a-1),cutoff*nmodes_a),
2122  pow((1.0*j)/(nmodes_b-1),cutoff*nmodes_b));
2123 
2124  for(int k = 0; k < nmodes_c-i-j; ++k)
2125  {
2126  NekDouble fac = std::max(fac1,
2127  pow((1.0*k)/(nmodes_c-1),cutoff*nmodes_c));
2128 
2129  orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2130  cnt++;
2131  }
2132  }
2133  }
2134  }
2135  else if(mkey.ConstFactorExists(eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
2136  {
2137  NekDouble SvvDiffCoeff =
2138  mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff)*
2139  mkey.GetConstFactor(eFactorSVVDiffCoeff);
2140 
2141  int max_abc = max(nmodes_a-kSVVDGFiltermodesmin,
2142  nmodes_b-kSVVDGFiltermodesmin);
2143  max_abc = max(max_abc, nmodes_c-kSVVDGFiltermodesmin);
2144  // clamp max_abc
2145  max_abc = max(max_abc,0);
2146  max_abc = min(max_abc,kSVVDGFiltermodesmax-kSVVDGFiltermodesmin);
2147 
2148  for(int i = 0; i < nmodes_a; ++i)
2149  {
2150  for(int j = 0; j < nmodes_b-j; ++j)
2151  {
2152  int maxij = max(i,j);
2153 
2154  for(int k = 0; k < nmodes_c-i-j; ++k)
2155  {
2156  int maxijk = max(maxij,k);
2157  maxijk = min(maxijk,kSVVDGFiltermodesmax-1);
2158 
2159  orthocoeffs[cnt] *= SvvDiffCoeff *
2160  kSVVDGFilter[max_abc][maxijk];
2161  cnt++;
2162  }
2163  }
2164  }
2165  }
2166  else
2167  {
2168 
2169  //SVV filter paramaters (how much added diffusion
2170  //relative to physical one and fraction of modes from
2171  //which you start applying this added diffusion)
2172 
2173  NekDouble SvvDiffCoeff = mkey.GetConstFactor(StdRegions::eFactorSVVDiffCoeff);
2174  NekDouble SVVCutOff = mkey.GetConstFactor(StdRegions::eFactorSVVCutoffRatio);
2175 
2176  //Defining the cut of mode
2177  int cutoff_a = (int) (SVVCutOff*nmodes_a);
2178  int cutoff_b = (int) (SVVCutOff*nmodes_b);
2179  int cutoff_c = (int) (SVVCutOff*nmodes_c);
2180  int nmodes = min(min(nmodes_a,nmodes_b),nmodes_c);
2181  NekDouble cutoff = min(min(cutoff_a,cutoff_b),cutoff_c);
2182  NekDouble epsilon = 1;
2183 
2184 
2185  //------"New" Version August 22nd '13--------------------
2186  for(i = 0; i < nmodes_a; ++i)
2187  {
2188  for(j = 0; j < nmodes_b-i; ++j)
2189  {
2190  for(k = 0; k < nmodes_c-i-j; ++k)
2191  {
2192  if(i + j + k >= cutoff)
2193  {
2194  orthocoeffs[cnt] *= ((SvvDiffCoeff)*exp(-(i+j+k-nmodes)*(i+j+k-nmodes)/((NekDouble)((i+j+k-cutoff+epsilon)*(i+j+k-cutoff+epsilon)))));
2195  }
2196  else
2197  {
2198  orthocoeffs[cnt] *= 0.0;
2199  }
2200  cnt++;
2201  }
2202  }
2203  }
2204  }
2205 
2206  // backward transform to physical space
2207  OrthoExp.BwdTrans(orthocoeffs,array);
2208  }
2209 
2210 
2211  void StdTetExp::v_ReduceOrderCoeffs(
2212  int numMin,
2213  const Array<OneD, const NekDouble> &inarray,
2214  Array<OneD, NekDouble> &outarray)
2215  {
2216  int nquad0 = m_base[0]->GetNumPoints();
2217  int nquad1 = m_base[1]->GetNumPoints();
2218  int nquad2 = m_base[2]->GetNumPoints();
2219  int nqtot = nquad0 * nquad1 * nquad2;
2220  int nmodes0 = m_base[0]->GetNumModes();
2221  int nmodes1 = m_base[1]->GetNumModes();
2222  int nmodes2 = m_base[2]->GetNumModes();
2223  int numMax = nmodes0;
2224 
2225  Array<OneD, NekDouble> coeff (m_ncoeffs);
2226  Array<OneD, NekDouble> coeff_tmp1(m_ncoeffs, 0.0);
2227  Array<OneD, NekDouble> coeff_tmp2(m_ncoeffs, 0.0);
2228  Array<OneD, NekDouble> phys_tmp (nqtot, 0.0);
2229  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
2230 
2231  Vmath::Vcopy(m_ncoeffs,inarray,1,coeff_tmp2,1);
2232 
2233  const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
2234  const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
2235  const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
2236 
2237  LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A,
2238  nmodes0, Pkey0);
2239  LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_B,
2240  nmodes1, Pkey1);
2241  LibUtilities::BasisKey bortho2(LibUtilities::eOrtho_C,
2242  nmodes2, Pkey2);
2243 
2244  Vmath::Zero(m_ncoeffs, coeff_tmp2, 1);
2245 
2246  StdRegions::StdTetExpSharedPtr OrthoTetExp;
2247  OrthoTetExp = MemoryManager<StdRegions::StdTetExp>
2248  ::AllocateSharedPtr(bortho0, bortho1, bortho2);
2249 
2250  BwdTrans(inarray,phys_tmp);
2251  OrthoTetExp->FwdTrans(phys_tmp, coeff);
2252 
2253  Vmath::Zero(m_ncoeffs,outarray,1);
2254 
2255  // filtering
2256  int cnt = 0;
2257  for (int u = 0; u < numMin; ++u)
2258  {
2259  for (int i = 0; i < numMin-u; ++i)
2260  {
2261  Vmath::Vcopy(numMin - u - i, tmp = coeff + cnt, 1,
2262  tmp2 = coeff_tmp1 + cnt, 1);
2263  cnt += numMax - u - i;
2264  }
2265  for (int i = numMin; i < numMax-u; ++i)
2266  {
2267  cnt += numMax - u - i;
2268  }
2269  }
2270 
2271  OrthoTetExp->BwdTrans(coeff_tmp1,phys_tmp);
2272  FwdTrans(phys_tmp, outarray);
2273  }
2274 
2275 
2276  void StdTetExp::v_GetSimplexEquiSpacedConnectivity(
2277  Array<OneD, int> &conn,
2278  bool standard)
2279  {
2280  boost::ignore_unused(standard);
2281 
2282  int np0 = m_base[0]->GetNumPoints();
2283  int np1 = m_base[1]->GetNumPoints();
2284  int np2 = m_base[2]->GetNumPoints();
2285  int np = max(np0,max(np1,np2));
2286 
2287 
2288  conn = Array<OneD, int>(4*(np-1)*(np-1)*(np-1));
2289 
2290  int row = 0;
2291  int rowp1 = 0;
2292  int plane = 0;
2293  int row1 = 0;
2294  int row1p1 = 0;
2295  int planep1= 0;
2296  int cnt = 0;
2297  for(int i = 0; i < np-1; ++i)
2298  {
2299  planep1 += (np-i)*(np-i+1)/2;
2300  row = 0; // current plane row offset
2301  rowp1 = 0; // current plane row plus one offset
2302  row1 = 0; // next plane row offset
2303  row1p1 = 0; // nex plane row plus one offset
2304  for(int j = 0; j < np-i-1; ++j)
2305  {
2306  rowp1 += np-i-j;
2307  row1p1 += np-i-j-1;
2308  for(int k = 0; k < np-i-j-2; ++k)
2309  {
2310  conn[cnt++] = plane + row +k+1;
2311  conn[cnt++] = plane + row +k;
2312  conn[cnt++] = plane + rowp1 +k;
2313  conn[cnt++] = planep1 + row1 +k;
2314 
2315  conn[cnt++] = plane + row +k+1;
2316  conn[cnt++] = plane + rowp1 +k+1;
2317  conn[cnt++] = planep1 + row1 +k+1;
2318  conn[cnt++] = planep1 + row1 +k;
2319 
2320  conn[cnt++] = plane + rowp1 +k+1;
2321  conn[cnt++] = plane + row +k+1;
2322  conn[cnt++] = plane + rowp1 +k;
2323  conn[cnt++] = planep1 + row1 +k;
2324 
2325  conn[cnt++] = planep1 + row1 +k;
2326  conn[cnt++] = planep1 + row1p1+k;
2327  conn[cnt++] = plane + rowp1 +k;
2328  conn[cnt++] = plane + rowp1 +k+1;
2329 
2330  conn[cnt++] = planep1 + row1 +k;
2331  conn[cnt++] = planep1 + row1p1+k;
2332  conn[cnt++] = planep1 + row1 +k+1;
2333  conn[cnt++] = plane + rowp1 +k+1;
2334 
2335  if(k < np-i-j-3)
2336  {
2337  conn[cnt++] = plane + rowp1 +k+1;
2338  conn[cnt++] = planep1 + row1p1 +k+1;
2339  conn[cnt++] = planep1 + row1 +k+1;
2340  conn[cnt++] = planep1 + row1p1 +k;
2341  }
2342  }
2343 
2344  conn[cnt++] = plane + row + np-i-j-1;
2345  conn[cnt++] = plane + row + np-i-j-2;
2346  conn[cnt++] = plane + rowp1 + np-i-j-2;
2347  conn[cnt++] = planep1 + row1 + np-i-j-2;
2348 
2349  row += np-i-j;
2350  row1 += np-i-j-1;
2351  }
2352  plane += (np-i)*(np-i+1)/2;
2353  }
2354  }
2355 
2356  }//end namespace
2357 }//end namespace
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:274
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::LibUtilities::BasisKey::GetNumModes
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:86
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:60
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:171
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:114
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:128
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:69
Nektar::StdRegions::StdExpansion
The base class for all shapes.
Definition: StdExpansion.h:63
Nektar::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:124
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:134
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:158
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:617
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:982
Nektar::StdRegions::StdExpansion::CreateGeneralMatrix
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
Definition: StdExpansion.cpp:258
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:433
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:208
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:1745
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:221
Nektar::StdRegions::StdExpansion::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:733
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:171
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1162
Nektar::StdRegions::StdMatrixKey::GetMatrixType
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:81
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:121
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:130
Nektar::StdRegions::StdTetExp::v_BwdTrans
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:294
Nektar::StdRegions::StdTetExp::v_BwdTrans_SumFac
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:324
Nektar::StdRegions::StdTetExp::v_NumBndryCoeffs
virtual int v_NumBndryCoeffs() const
Definition: StdTetExp.cpp:1030
Nektar::StdRegions::StdTetExp::v_GetNedges
virtual int v_GetNedges() const
Definition: StdTetExp.cpp:1015
Nektar::StdRegions::StdTetExp::v_GetSimplexEquiSpacedConnectivity
virtual void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true)
Definition: StdTetExp.cpp:2276
Nektar::StdRegions::StdTetExp::v_GetBoundaryMap
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdTetExp.cpp:1395
Nektar::StdRegions::StdTetExp::v_GetTraceToElementMap
virtual void v_GetTraceToElementMap(const int fid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation faceOrient, int P, int Q)
Definition: StdTetExp.cpp:1454
Nektar::StdRegions::StdTetExp::v_ReduceOrderCoeffs
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:2211
Nektar::StdRegions::StdTetExp::v_IProductWRTBase_MatOp
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:530
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:252
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:585
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:1960
Nektar::StdRegions::StdTetExp::v_CalcNumberOfCoefficients
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdTetExp.cpp:1193
Nektar::StdRegions::StdTetExp::v_GetNtraces
virtual int v_GetNtraces() const
Definition: StdTetExp.cpp:1020
Nektar::StdRegions::StdTetExp::v_FwdTrans
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:456
Nektar::StdRegions::StdTetExp::v_PhysEvaluateBasis
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) final
Definition: StdTetExp.cpp:924
Nektar::StdRegions::StdTetExp::v_GetTraceNumPoints
virtual int v_GetTraceNumPoints(const int i) const
Definition: StdTetExp.cpp:1130
Nektar::StdRegions::StdTetExp::v_MultiplyByStdQuadratureMetric
virtual void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:1992
Nektar::StdRegions::StdTetExp::v_GetEdgeNcoeffs
virtual int v_GetEdgeNcoeffs(const int i) const
Definition: StdTetExp.cpp:1152
Nektar::StdRegions::StdTetExp::v_GetTotalTraceIntNcoeffs
virtual int v_GetTotalTraceIntNcoeffs() const
Definition: StdTetExp.cpp:1119
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:674
Nektar::StdRegions::StdTetExp::v_GetTraceNumModes
virtual void v_GetTraceNumModes(const int fid, int &numModes0, int &numModes1, Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
Definition: StdTetExp.cpp:970
Nektar::StdRegions::StdTetExp::v_GetTraceNcoeffs
virtual int v_GetTraceNcoeffs(const int i) const
Definition: StdTetExp.cpp:1072
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:549
Nektar::StdRegions::StdTetExp::v_FillMode
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:915
Nektar::StdRegions::StdTetExp::v_GetVertexMap
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdTetExp.cpp:1274
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:1237
Nektar::StdRegions::StdTetExp::v_IsBoundaryInteriorExpansion
virtual bool v_IsBoundaryInteriorExpansion()
Definition: StdTetExp.cpp:1263
Nektar::StdRegions::StdTetExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdTetExp.cpp:1356
Nektar::StdRegions::StdTetExp::v_LocCoordToLocCollapsed
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Definition: StdTetExp.cpp:873
Nektar::StdRegions::StdTetExp::v_GetTraceIntNcoeffs
virtual int v_GetTraceIntNcoeffs(const int i) const
Definition: StdTetExp.cpp:1098
Nektar::StdRegions::StdTetExp::v_IProductWRTBase
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:507
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:683
Nektar::StdRegions::StdTetExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Definition: StdTetExp.cpp:1841
Nektar::StdRegions::StdTetExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdTetExp.cpp:2072
Nektar::StdRegions::StdTetExp::v_GetNverts
virtual int v_GetNverts() const
Definition: StdTetExp.cpp:1010
Nektar::StdRegions::StdTetExp::v_GetTracePointsKey
virtual LibUtilities::PointsKey v_GetTracePointsKey(const int i, const int j) const
Definition: StdTetExp.cpp:1173
Nektar::StdRegions::StdTetExp::v_GetEdgeInteriorToElementMap
virtual void v_GetEdgeInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
Definition: StdTetExp.cpp:1631
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:107
Nektar::StdRegions::StdTetExp::v_GetTraceInteriorToElementMap
virtual void v_GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
Definition: StdTetExp.cpp:1748
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:357
Nektar::StdRegions::StdTetExp::v_LocCollapsedToLocCoord
virtual void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi)
Definition: StdTetExp.cpp:906
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:720
Nektar::StdRegions::StdTetExp::DetShapeType
LibUtilities::ShapeType DetShapeType() const
Definition: StdTetExp.h:69
Nektar::StdRegions::StdTetExp::v_GetTraceBasisKey
virtual const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int k) const
Definition: StdTetExp.cpp:1206
Nektar::StdRegions::StdTetExp::v_NumDGBndryCoeffs
virtual int v_NumDGBndryCoeffs() const
Definition: StdTetExp.cpp:1050
Nektar::StdRegions::StdTetExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdTetExp.cpp:1929
Nektar::StdRegions::StdTetExp::v_DetShapeType
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdTetExp.cpp:1025
Blas::Dgemv
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = A x where A[m x n].
Definition: Blas.hpp:265
Blas::Dscal
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:182
Blas::Ddot
static double Ddot(const int &n, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: output = .
Definition: Blas.hpp:197
Blas::Dgemm
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition: Blas.hpp:394
Blas::Daxpy
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
Definition: Blas.hpp:167
Nektar::LibUtilities::StdTetData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:217
Nektar::LibUtilities::StdTetData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:194
Nektar::LibUtilities::NullBasisKey
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
Nektar::LibUtilities::eGaussRadauMAlpha2Beta0
@ eGaussRadauMAlpha2Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:59
Nektar::LibUtilities::eGaussRadauMAlpha1Beta0
@ eGaussRadauMAlpha1Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
Nektar::LibUtilities::eModified_B
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::LibUtilities::eOrtho_A
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:45
Nektar::LibUtilities::eModified_C
@ eModified_C
Principle Modified Functions .
Definition: BasisType.h:50
Nektar::LibUtilities::eGLL_Lagrange
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:54
Nektar::LibUtilities::eOrtho_C
@ eOrtho_C
Principle Orthogonal Functions .
Definition: BasisType.h:47
Nektar::LibUtilities::eOrtho_B
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:46
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:51
Nektar::StdRegions::EvaluateTriFaceBasisKey
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Definition: StdExpansion3D.cpp:499
Nektar::StdRegions::StdTetExpSharedPtr
std::shared_ptr< StdTetExp > StdTetExpSharedPtr
Definition: StdTetExp.h:279
Nektar::StdRegions::kSVVDGFiltermodesmin
const int kSVVDGFiltermodesmin
Definition: StdRegions.hpp:388
Nektar::StdRegions::kSVVDGFiltermodesmax
const int kSVVDGFiltermodesmax
Definition: StdRegions.hpp:389
Nektar::StdRegions::kSVVDGFilter
const NekDouble kSVVDGFilter[9][11]
Definition: StdRegions.hpp:391
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:820
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
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:192
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:322
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:436
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:341
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
main.P
P
Definition: main.py:133
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267
Generated on Wed Nov 24 2021 10:33:29 for Nektar++ by doxygen 1.9.1