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.num_elements() > 0)? true:false;
124  bool Do_1 = (out_dxi1.num_elements() > 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.num_elements() > 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  if( fabs(xi[2]-1.0) < NekConstants::kNekZeroTol)
878  {
879  // Very top point of the tetrahedron
880  eta[0] = -1.0;
881  eta[1] = -1.0;
882  eta[2] = xi[2];
883  }
884  else
885  {
886  if( fabs(xi[1]-1.0) < NekConstants::kNekZeroTol )
887  {
888  // Distant diagonal edge shared by all eta_x
889  // coordinate planes: the xi_y == -xi_z line
890  eta[0] = -1.0;
891  }
892  else if (fabs(xi[1] + xi[2]) < NekConstants::kNekZeroTol)
893  {
894  eta[0] = -1.0;
895  }
896  else
897  {
898  eta[0] = 2.0*(1.0+xi[0])/(-xi[1]-xi[2]) - 1.0;
899  }
900  eta[1] = 2.0*(1.0+xi[1])/(1.0-xi[2]) - 1.0;
901  eta[2] = xi[2];
902  }
903  }
904 
905  void StdTetExp::v_FillMode(
906  const int mode,
907  Array<OneD, NekDouble> &outarray)
908  {
909  Array<OneD, NekDouble> tmp(m_ncoeffs,0.0);
910  tmp[mode] = 1.0;
911  StdTetExp::v_BwdTrans(tmp, outarray);
912  }
913 
914  void StdTetExp::v_GetFaceNumModes(
915  const int fid,
916  const Orientation faceOrient,
917  int &numModes0,
918  int &numModes1)
919  {
920  boost::ignore_unused(faceOrient);
921 
922  int nummodes [3] = {m_base[0]->GetNumModes(),
923  m_base[1]->GetNumModes(),
924  m_base[2]->GetNumModes()};
925  switch(fid)
926  {
927  case 0:
928  {
929  numModes0 = nummodes[0];
930  numModes1 = nummodes[1];
931  }
932  break;
933  case 1:
934  {
935  numModes0 = nummodes[0];
936  numModes1 = nummodes[2];
937  }
938  break;
939  case 2:
940  case 3:
941  {
942  numModes0 = nummodes[1];
943  numModes1 = nummodes[2];
944  }
945  break;
946  }
947  }
948 
949  //---------------------------
950  // Helper functions
951  //---------------------------
952 
953  int StdTetExp::v_GetNverts() const
954  {
955  return 4;
956  }
957 
958  int StdTetExp::v_GetNedges() const
959  {
960  return 6;
961  }
962 
963  int StdTetExp::v_GetNfaces() const
964  {
965  return 4;
966  }
967 
968  LibUtilities::ShapeType StdTetExp::v_DetShapeType() const
969  {
970  return DetShapeType();
971  }
972 
973  int StdTetExp::v_NumBndryCoeffs() const
974  {
975  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
976  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
977  "BasisType is not a boundary interior form");
978  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
979  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
980  "BasisType is not a boundary interior form");
981  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
982  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
983  "BasisType is not a boundary interior form");
984 
985  int P = m_base[0]->GetNumModes();
986  int Q = m_base[1]->GetNumModes();
987  int R = m_base[2]->GetNumModes();
988 
989  return LibUtilities::StdTetData::
990  getNumberOfBndCoefficients(P, Q, R);
991  }
992 
993  int StdTetExp::v_NumDGBndryCoeffs() const
994  {
995  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
996  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
997  "BasisType is not a boundary interior form");
998  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
999  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1000  "BasisType is not a boundary interior form");
1001  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1002  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1003  "BasisType is not a boundary interior form");
1004 
1005  int P = m_base[0]->GetNumModes()-1;
1006  int Q = m_base[1]->GetNumModes()-1;
1007  int R = m_base[2]->GetNumModes()-1;
1008 
1009 
1010  return (Q+1) + P*(1 + 2*Q - P)/2 // base face
1011  + (R+1) + P*(1 + 2*R - P)/2 // front face
1012  + 2*(R+1) + Q*(1 + 2*R - Q); // back two faces
1013  }
1014 
1015  int StdTetExp::v_GetEdgeNcoeffs(const int i) const
1016  {
1017  ASSERTL2((i >= 0) && (i <= 5), "edge id is out of range");
1018  int P = m_base[0]->GetNumModes();
1019  int Q = m_base[1]->GetNumModes();
1020  int R = m_base[2]->GetNumModes();
1021 
1022  if (i == 0)
1023  {
1024  return P;
1025  }
1026  else if (i == 1 || i == 2)
1027  {
1028  return Q;
1029  }
1030  else
1031  {
1032  return R;
1033  }
1034  }
1035 
1036  int StdTetExp::v_GetTotalEdgeIntNcoeffs() const
1037  {
1038  int P = m_base[0]->GetNumModes()-2;
1039  int Q = m_base[1]->GetNumModes()-2;
1040  int R = m_base[2]->GetNumModes()-2;
1041 
1042  return P+Q+4*R;
1043  }
1044 
1045  int StdTetExp::v_GetFaceNcoeffs(const int i) const
1046  {
1047  ASSERTL2((i >= 0) && (i <= 3), "face id is out of range");
1048  int nFaceCoeffs = 0;
1049  int nummodesA, nummodesB, P, Q;
1050  if (i == 0)
1051  {
1052  nummodesA = GetBasisNumModes(0);
1053  nummodesB = GetBasisNumModes(1);
1054  }
1055  else if ((i == 1) || (i == 2))
1056  {
1057  nummodesA = GetBasisNumModes(0);
1058  nummodesB = GetBasisNumModes(2);
1059  }
1060  else
1061  {
1062  nummodesA = GetBasisNumModes(1);
1063  nummodesB = GetBasisNumModes(2);
1064  }
1065  P = nummodesA - 1;
1066  Q = nummodesB - 1;
1067  nFaceCoeffs = Q+1 + (P*(1 + 2*Q - P))/2;
1068  return nFaceCoeffs;
1069  }
1070 
1071  int StdTetExp::v_GetFaceIntNcoeffs(const int i) const
1072  {
1073  ASSERTL2((i >= 0) && (i <= 3), "face id is out of range");
1074  int Pi = m_base[0]->GetNumModes() - 2;
1075  int Qi = m_base[1]->GetNumModes() - 2;
1076  int Ri = m_base[2]->GetNumModes() - 2;
1077 
1078  if((i == 0))
1079  {
1080  return Pi * (2*Qi - Pi - 1) / 2;
1081  }
1082  else if((i == 1))
1083  {
1084  return Pi * (2*Ri - Pi - 1) / 2;
1085  }
1086  else
1087  {
1088  return Qi * (2*Ri - Qi - 1) / 2;
1089  }
1090  }
1091 
1092  int StdTetExp::v_GetTotalFaceIntNcoeffs() const
1093  {
1094  int Pi = m_base[0]->GetNumModes() - 2;
1095  int Qi = m_base[1]->GetNumModes() - 2;
1096  int Ri = m_base[2]->GetNumModes() - 2;
1097 
1098  return Pi * (2*Qi - Pi - 1) / 2 +
1099  Pi * (2*Ri - Pi - 1) / 2 +
1100  Qi * (2*Ri - Qi - 1);
1101  }
1102 
1103  int StdTetExp::v_GetFaceNumPoints(const int i) const
1104  {
1105  ASSERTL2(i >= 0 && i <= 3, "face id is out of range");
1106 
1107  if (i == 0)
1108  {
1109  return m_base[0]->GetNumPoints()*
1110  m_base[1]->GetNumPoints();
1111  }
1112  else if (i == 1)
1113  {
1114  return m_base[0]->GetNumPoints()*
1115  m_base[2]->GetNumPoints();
1116  }
1117  else
1118  {
1119  return m_base[1]->GetNumPoints()*
1120  m_base[2]->GetNumPoints();
1121  }
1122  }
1123 
1124  LibUtilities::PointsKey StdTetExp::v_GetFacePointsKey(
1125  const int i, const int j) const
1126  {
1127  ASSERTL2(i >= 0 && i <= 3, "face id is out of range");
1128  ASSERTL2(j == 0 || j == 1, "face direction is out of range");
1129 
1130  if (i == 0)
1131  {
1132  return m_base[j]->GetPointsKey();
1133  }
1134  else if (i == 1)
1135  {
1136  return m_base[2*j]->GetPointsKey();
1137  }
1138  else
1139  {
1140  return m_base[j+1]->GetPointsKey();
1141  }
1142  }
1143 
1144  int StdTetExp::v_CalcNumberOfCoefficients(
1145  const std::vector<unsigned int>& nummodes,
1146  int & modes_offset)
1147  {
1148  int nmodes = LibUtilities::StdTetData::getNumberOfCoefficients(
1149  nummodes[modes_offset],
1150  nummodes[modes_offset+1],
1151  nummodes[modes_offset+2]);
1152  modes_offset += 3;
1153 
1154  return nmodes;
1155  }
1156 
1157  const LibUtilities::BasisKey StdTetExp::v_DetFaceBasisKey(
1158  const int i, const int k) const
1159  {
1160  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
1161  ASSERTL2(k == 0 || k == 1, "face direction out of range");
1162 
1163  int dir = k;
1164  switch(i)
1165  {
1166  case 0:
1167  dir = k;
1168  break;
1169  case 1:
1170  dir = 2*k;
1171  break;
1172  case 2:
1173  case 3:
1174  dir = k+1;
1175  break;
1176  }
1177 
1178  return EvaluateTriFaceBasisKey(k,
1179  m_base[dir]->GetBasisType(),
1180  m_base[dir]->GetNumPoints(),
1181  m_base[dir]->GetNumModes());
1182 
1183  // Should not get here.
1184  return LibUtilities::NullBasisKey;
1185  }
1186 
1187  LibUtilities::BasisType StdTetExp::v_GetEdgeBasisType(const int i) const
1188  {
1189  ASSERTL2(i >= 0 && i <= 5, "edge id is out of range");
1190 
1191  if (i == 0)
1192  {
1193  return GetBasisType(0);
1194  }
1195  else if (i == 1 || i == 2)
1196  {
1197  return GetBasisType(1);
1198  }
1199  else
1200  {
1201  return GetBasisType(2);
1202  }
1203  }
1204 
1205  void StdTetExp::v_GetCoords(
1206  Array<OneD, NekDouble> &xi_x,
1207  Array<OneD, NekDouble> &xi_y,
1208  Array<OneD, NekDouble> &xi_z)
1209  {
1210  Array<OneD, const NekDouble> eta_x = m_base[0]->GetZ();
1211  Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
1212  Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
1213  int Qx = GetNumPoints(0);
1214  int Qy = GetNumPoints(1);
1215  int Qz = GetNumPoints(2);
1216 
1217  // Convert collapsed coordinates into cartesian coordinates: eta
1218  // --> xi
1219  for( int k = 0; k < Qz; ++k ) {
1220  for( int j = 0; j < Qy; ++j ) {
1221  for( int i = 0; i < Qx; ++i ) {
1222  int s = i + Qx*(j + Qy*k);
1223  xi_x[s] = (eta_x[i] + 1.0) * (1.0 - eta_y[j]) * (1.0 - eta_z[k]) / 4 - 1.0;
1224  xi_y[s] = (eta_y[j] + 1.0) * (1.0 - eta_z[k]) / 2 - 1.0;
1225  xi_z[s] = eta_z[k];
1226  }
1227  }
1228  }
1229  }
1230 
1231  bool StdTetExp::v_IsBoundaryInteriorExpansion()
1232  {
1233  return (m_base[0]->GetBasisType() == LibUtilities::eModified_A) &&
1234  (m_base[1]->GetBasisType() == LibUtilities::eModified_B) &&
1235  (m_base[2]->GetBasisType() == LibUtilities::eModified_C);
1236  }
1237 
1238 
1239  //--------------------------
1240  // Mappings
1241  //--------------------------
1242 
1243  void StdTetExp::v_GetEdgeToElementMap(
1244  const int eid,
1245  const Orientation edgeOrient,
1246  Array<OneD, unsigned int>& maparray,
1247  Array<OneD, int>& signarray,
1248  int P)
1249  {
1250  boost::ignore_unused(P);
1251 
1252  ASSERTL2(eid >= 0 && eid < 6, "Invalid edge");
1253  ASSERTL2(v_IsBoundaryInteriorExpansion(),
1254  "Method only implemented for Modified_A BasisType (x "
1255  "direction), Modified_B BasisType (y direction), and "
1256  "Modified_C BasisType(z direction)");
1257 
1258  int edgeVerts[6][2] = {{0,1},{1,2},{0,2},{0,3},{1,3},{2,3}};
1259  int nEdgeModes = v_GetEdgeNcoeffs(eid);
1260 
1261  if (maparray.num_elements() != nEdgeModes)
1262  {
1263  maparray = Array<OneD, unsigned int>(nEdgeModes);
1264  }
1265 
1266  if (signarray.num_elements() != nEdgeModes)
1267  {
1268  signarray = Array<OneD, int>(nEdgeModes, 1.0);
1269  }
1270 
1271  if (edgeOrient == StdRegions::eForwards)
1272  {
1273  maparray[0] = v_GetVertexMap(edgeVerts[eid][0]);
1274  maparray[1] = v_GetVertexMap(edgeVerts[eid][1]);
1275  }
1276  else if (edgeOrient == StdRegions::eBackwards)
1277  {
1278  maparray[0] = v_GetVertexMap(edgeVerts[eid][1]);
1279  maparray[1] = v_GetVertexMap(edgeVerts[eid][0]);
1280  }
1281  else
1282  {
1283  ASSERTL2(false, "Unsupported edge orientation");
1284  }
1285 
1286  Array<OneD, unsigned int> tmp1(nEdgeModes-2);
1287  Array<OneD, int> tmp2(nEdgeModes-2);
1288  v_GetEdgeInteriorMap(eid, edgeOrient, tmp1, tmp2);
1289 
1290  for (int i = 0; i < nEdgeModes-2; ++i)
1291  {
1292  maparray[i+2] = tmp1[i];
1293  signarray[i+2] = tmp2[i];
1294  }
1295  }
1296 
1297  /**
1298  * Maps Expansion2D modes of a 2D face to the corresponding expansion
1299  * modes.
1300  */
1301  void StdTetExp::v_GetFaceToElementMap(
1302  const int fid,
1303  const Orientation faceOrient,
1304  Array<OneD, unsigned int> &maparray,
1305  Array<OneD, int> &signarray,
1306  int P,
1307  int Q)
1308  {
1309  int nummodesA=0,nummodesB=0, i, j, k, idx;
1310 
1311  ASSERTL1(v_IsBoundaryInteriorExpansion(),
1312  "Method only implemented for Modified_A BasisType (x "
1313  "direction), Modified_B BasisType (y direction), and "
1314  "Modified_C BasisType(z direction)");
1315 
1316  int nFaceCoeffs = 0;
1317 
1318  switch(fid)
1319  {
1320  case 0:
1321  nummodesA = m_base[0]->GetNumModes();
1322  nummodesB = m_base[1]->GetNumModes();
1323  break;
1324  case 1:
1325  nummodesA = m_base[0]->GetNumModes();
1326  nummodesB = m_base[2]->GetNumModes();
1327  break;
1328  case 2:
1329  case 3:
1330  nummodesA = m_base[1]->GetNumModes();
1331  nummodesB = m_base[2]->GetNumModes();
1332  break;
1333  default:
1334  ASSERTL0(false,"fid must be between 0 and 3");
1335  }
1336 
1337  bool CheckForZeroedModes = false;
1338  if(P == -1)
1339  {
1340  P = nummodesA;
1341  Q = nummodesB;
1342  }
1343  else
1344  {
1345  CheckForZeroedModes = true;
1346  }
1347 
1348  nFaceCoeffs = P*(2*Q-P+1)/2;
1349 
1350  // Allocate the map array and sign array; set sign array to ones (+)
1351  if(maparray.num_elements() != nFaceCoeffs)
1352  {
1353  maparray = Array<OneD, unsigned int>(nFaceCoeffs,1);
1354  }
1355 
1356  if(signarray.num_elements() != nFaceCoeffs)
1357  {
1358  signarray = Array<OneD, int>(nFaceCoeffs,1);
1359  }
1360  else
1361  {
1362  fill(signarray.get(),signarray.get()+nFaceCoeffs, 1 );
1363  }
1364 
1365  switch (fid)
1366  {
1367  case 0:
1368  idx = 0;
1369  for (i = 0; i < P; ++i)
1370  {
1371  for (j = 0; j < Q-i; ++j)
1372  {
1373  if ((int)faceOrient == 7 && i > 1)
1374  {
1375  signarray[idx] = (i%2 ? -1 : 1);
1376  }
1377  maparray[idx++] = GetMode(i,j,0);
1378  }
1379  }
1380  break;
1381  case 1:
1382  idx = 0;
1383  for (i = 0; i < P; ++i)
1384  {
1385  for (k = 0; k < Q-i; ++k)
1386  {
1387  if ((int)faceOrient == 7 && i > 1)
1388  {
1389  signarray[idx] = (i%2 ? -1: 1);
1390  }
1391  maparray[idx++] = GetMode(i,0,k);
1392  }
1393  }
1394  break;
1395  case 2:
1396  idx = 0;
1397  for (j = 0; j < P-1; ++j)
1398  {
1399  for (k = 0; k < Q-1-j; ++k)
1400  {
1401  if ((int)faceOrient == 7 && j > 1)
1402  {
1403  signarray[idx] = ((j+1)%2 ? -1: 1);
1404  }
1405  maparray[idx++] = GetMode(1,j,k);
1406  // Incorporate modes from zeroth plane where needed.
1407  if (j == 0 && k == 0)
1408  {
1409  maparray[idx++] = GetMode(0,0,1);
1410  }
1411  if (j == 0 && k == Q-2)
1412  {
1413  for (int r = 0; r < Q-1; ++r)
1414  {
1415  maparray[idx++] = GetMode(0,1,r);
1416  }
1417  }
1418  }
1419  }
1420  break;
1421  case 3:
1422  idx = 0;
1423  for (j = 0; j < P; ++j)
1424  {
1425  for (k = 0; k < Q-j; ++k)
1426  {
1427  if ((int)faceOrient == 7 && j > 1)
1428  {
1429  signarray[idx] = (j%2 ? -1: 1);
1430  }
1431  maparray[idx++] = GetMode(0,j,k);
1432  }
1433  }
1434  break;
1435  default:
1436  ASSERTL0(false, "Element map not available.");
1437  }
1438 
1439  if ((int)faceOrient == 7)
1440  {
1441  swap(maparray[0], maparray[Q]);
1442 
1443  for (i = 1; i < Q-1; ++i)
1444  {
1445  swap(maparray[i+1], maparray[Q+i]);
1446  }
1447  }
1448 
1449  if(CheckForZeroedModes)
1450  {
1451  // zero signmap and set maparray to zero if elemental
1452  // modes are not as large as face modesl
1453  idx = 0;
1454  for (j = 0; j < nummodesA; ++j)
1455  {
1456  idx += nummodesB-j;
1457  for (k = nummodesB-j; k < Q-j; ++k)
1458  {
1459  signarray[idx] = 0.0;
1460  maparray[idx++] = maparray[0];
1461  }
1462  }
1463 
1464  for (j = nummodesA; j < P; ++j)
1465  {
1466  for (k = 0; k < Q-j; ++k)
1467  {
1468  signarray[idx] = 0.0;
1469  maparray[idx++] = maparray[0];
1470  }
1471  }
1472  }
1473  }
1474 
1475  int StdTetExp::v_GetVertexMap(const int localVertexId, bool useCoeffPacking)
1476  {
1477  ASSERTL0((GetEdgeBasisType(localVertexId)==LibUtilities::eModified_A)||
1478  (GetEdgeBasisType(localVertexId)==LibUtilities::eModified_B)||
1479  (GetEdgeBasisType(localVertexId)==LibUtilities::eModified_C),
1480  "Mapping not defined for this type of basis");
1481 
1482  int localDOF = 0;
1483  if(useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1484  {
1485  switch(localVertexId)
1486  {
1487  case 0:
1488  {
1489  localDOF = GetMode(0,0,0);
1490  break;
1491  }
1492  case 1:
1493  {
1494  localDOF = GetMode(0,0,1);
1495  break;
1496  }
1497  case 2:
1498  {
1499  localDOF = GetMode(0,1,0);
1500  break;
1501  }
1502  case 3:
1503  {
1504  localDOF = GetMode(1,0,0);
1505  break;
1506  }
1507  default:
1508  {
1509  ASSERTL0(false,"Vertex ID must be between 0 and 3");
1510  break;
1511  }
1512  }
1513  }
1514  else
1515  {
1516  switch(localVertexId)
1517  {
1518  case 0:
1519  {
1520  localDOF = GetMode(0,0,0);
1521  break;
1522  }
1523  case 1:
1524  {
1525  localDOF = GetMode(1,0,0);
1526  break;
1527  }
1528  case 2:
1529  {
1530  localDOF = GetMode(0,1,0);
1531  break;
1532  }
1533  case 3:
1534  {
1535  localDOF = GetMode(0,0,1);
1536  break;
1537  }
1538  default:
1539  {
1540  ASSERTL0(false,"Vertex ID must be between 0 and 3");
1541  break;
1542  }
1543  }
1544 
1545  }
1546 
1547  return localDOF;
1548  }
1549 
1550  /**
1551  * Maps interior modes of an edge to the elemental modes.
1552  */
1553  void StdTetExp::v_GetEdgeInteriorMap(
1554  const int eid,
1555  const Orientation edgeOrient,
1556  Array<OneD, unsigned int> &maparray,
1557  Array<OneD, int> &signarray)
1558  {
1559  int i;
1560  const int P = m_base[0]->GetNumModes();
1561  const int Q = m_base[1]->GetNumModes();
1562  const int R = m_base[2]->GetNumModes();
1563 
1564  const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid)-2;
1565 
1566  if(maparray.num_elements() != nEdgeIntCoeffs)
1567  {
1568  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1569  }
1570  else
1571  {
1572  fill( maparray.get(), maparray.get() + nEdgeIntCoeffs, 0);
1573  }
1574 
1575  if(signarray.num_elements() != nEdgeIntCoeffs)
1576  {
1577  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1578  }
1579  else
1580  {
1581  fill( signarray.get() , signarray.get()+nEdgeIntCoeffs, 1 );
1582  }
1583 
1584  switch (eid)
1585  {
1586  case 0:
1587  for (i = 0; i < P-2; ++i)
1588  {
1589  maparray[i] = GetMode(i+2, 0, 0);
1590  }
1591  if(edgeOrient==eBackwards)
1592  {
1593  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1594  {
1595  signarray[i] = -1;
1596  }
1597  }
1598  break;
1599  case 1:
1600  for (i = 0; i < Q-2; ++i)
1601  {
1602  maparray[i] = GetMode(1, i+1, 0);
1603  }
1604  if(edgeOrient==eBackwards)
1605  {
1606  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1607  {
1608  signarray[i] = -1;
1609  }
1610  }
1611  break;
1612  case 2:
1613  for (i = 0; i < Q-2; ++i)
1614  {
1615  maparray[i] = GetMode(0, i+2, 0);
1616  }
1617  if(edgeOrient==eBackwards)
1618  {
1619  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1620  {
1621  signarray[i] = -1;
1622  }
1623  }
1624  break;
1625  case 3:
1626  for (i = 0; i < R-2; ++i)
1627  {
1628  maparray[i] = GetMode(0, 0, i+2);
1629  }
1630  if(edgeOrient==eBackwards)
1631  {
1632  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1633  {
1634  signarray[i] = -1;
1635  }
1636  }
1637  break;
1638  case 4:
1639  for (i = 0; i < R - 2; ++i)
1640  {
1641  maparray[i] = GetMode(1, 0, i+1);
1642  }
1643  if(edgeOrient==eBackwards)
1644  {
1645  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1646  {
1647  signarray[i] = -1;
1648  }
1649  }
1650  break;
1651  case 5:
1652  for (i = 0; i < R-2; ++i)
1653  {
1654  maparray[i] = GetMode(0, 1, i+1);
1655  }
1656  if(edgeOrient==eBackwards)
1657  {
1658  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1659  {
1660  signarray[i] = -1;
1661  }
1662  }
1663  break;
1664  default:
1665  ASSERTL0(false, "Edge not defined.");
1666  break;
1667  }
1668  }
1669 
1670  void StdTetExp::v_GetFaceInteriorMap(
1671  const int fid,
1672  const Orientation faceOrient,
1673  Array<OneD, unsigned int> &maparray,
1674  Array<OneD, int> &signarray)
1675  {
1676  int i, j, idx, k;
1677  const int P = m_base[0]->GetNumModes();
1678  const int Q = m_base[1]->GetNumModes();
1679  const int R = m_base[2]->GetNumModes();
1680 
1681  const int nFaceIntCoeffs = v_GetFaceIntNcoeffs(fid);
1682 
1683  if(maparray.num_elements() != nFaceIntCoeffs)
1684  {
1685  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1686  }
1687 
1688  if(signarray.num_elements() != nFaceIntCoeffs)
1689  {
1690  signarray = Array<OneD, int>(nFaceIntCoeffs,1);
1691  }
1692  else
1693  {
1694  fill( signarray.get() , signarray.get()+nFaceIntCoeffs, 1 );
1695  }
1696 
1697  switch (fid)
1698  {
1699  case 0:
1700  idx = 0;
1701  for (i = 2; i < P-1; ++i)
1702  {
1703  for (j = 1; j < Q-i; ++j)
1704  {
1705  if ((int)faceOrient == 7)
1706  {
1707  signarray[idx] = (i%2 ? -1 : 1);
1708  }
1709  maparray[idx++] = GetMode(i,j,0);
1710  }
1711  }
1712  break;
1713  case 1:
1714  idx = 0;
1715  for (i = 2; i < P; ++i)
1716  {
1717  for (k = 1; k < R-i; ++k)
1718  {
1719  if ((int)faceOrient == 7)
1720  {
1721  signarray[idx] = (i%2 ? -1: 1);
1722  }
1723  maparray[idx++] = GetMode(i,0,k);
1724  }
1725  }
1726  break;
1727  case 2:
1728  idx = 0;
1729  for (j = 1; j < Q-2; ++j)
1730  {
1731  for (k = 1; k < R-1-j; ++k)
1732  {
1733  if ((int)faceOrient == 7)
1734  {
1735  signarray[idx] = ((j+1)%2 ? -1: 1);
1736  }
1737  maparray[idx++] = GetMode(1,j,k);
1738  }
1739  }
1740  break;
1741  case 3:
1742  idx = 0;
1743  for (j = 2; j < Q-1; ++j)
1744  {
1745  for (k = 1; k < R-j; ++k)
1746  {
1747  if ((int)faceOrient == 7)
1748  {
1749  signarray[idx] = (j%2 ? -1: 1);
1750  }
1751  maparray[idx++] = GetMode(0,j,k);
1752  }
1753  }
1754  break;
1755  default:
1756  ASSERTL0(false, "Face interior map not available.");
1757  break;
1758  }
1759  }
1760 
1761  /**
1762  * List of all interior modes in the expansion.
1763  */
1764  void StdTetExp::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
1765  {
1766  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1767  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1768  "BasisType is not a boundary interior form");
1769  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
1770  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1771  "BasisType is not a boundary interior form");
1772  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1773  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1774  "BasisType is not a boundary interior form");
1775 
1776  int P = m_base[0]->GetNumModes();
1777  int Q = m_base[1]->GetNumModes();
1778  int R = m_base[2]->GetNumModes();
1779 
1780  int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1781 
1782  if(outarray.num_elements() != nIntCoeffs)
1783  {
1784  outarray = Array<OneD, unsigned int>(nIntCoeffs);
1785  }
1786 
1787  int idx = 0;
1788  for (int i = 2; i < P-2; ++i)
1789  {
1790  for (int j = 1; j < Q-i-1; ++j)
1791  {
1792  for (int k = 1; k < R-i-j; ++k)
1793  {
1794  outarray[idx++] = GetMode(i,j,k);
1795  }
1796  }
1797  }
1798  }
1799 
1800  /**
1801  * List of all boundary modes in the the expansion.
1802  */
1803  void StdTetExp::v_GetBoundaryMap(Array<OneD, unsigned int>& outarray)
1804  {
1805  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1806  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1807  "BasisType is not a boundary interior form");
1808  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
1809  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1810  "BasisType is not a boundary interior form");
1811  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1812  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1813  "BasisType is not a boundary interior form");
1814 
1815  int P = m_base[0]->GetNumModes();
1816  int Q = m_base[1]->GetNumModes();
1817  int R = m_base[2]->GetNumModes();
1818 
1819  int i,j,k;
1820  int idx = 0;
1821 
1822  int nBnd = NumBndryCoeffs();
1823 
1824  if (outarray.num_elements() != nBnd)
1825  {
1826  outarray = Array<OneD, unsigned int>(nBnd);
1827  }
1828 
1829  for (i = 0; i < P; ++i)
1830  {
1831  // First two Q-R planes are entirely boundary modes
1832  if (i < 2)
1833  {
1834  for (j = 0; j < Q-i; j++)
1835  {
1836  for (k = 0; k < R-i-j; ++k)
1837  {
1838  outarray[idx++] = GetMode(i,j,k);
1839  }
1840  }
1841  }
1842  // Remaining Q-R planes contain boundary modes on bottom and
1843  // left edge.
1844  else
1845  {
1846  for (k = 0; k < R-i; ++k)
1847  {
1848  outarray[idx++] = GetMode(i,0,k);
1849  }
1850  for (j = 1; j < Q-i; ++j)
1851  {
1852  outarray[idx++] = GetMode(i,j,0);
1853  }
1854  }
1855  }
1856  }
1857 
1858 
1859  //---------------------------------------
1860  // Wrapper functions
1861  //---------------------------------------
1862  DNekMatSharedPtr StdTetExp::v_GenMatrix(const StdMatrixKey &mkey)
1863  {
1864 
1865  MatrixType mtype = mkey.GetMatrixType();
1866 
1867  DNekMatSharedPtr Mat;
1868 
1869  switch(mtype)
1870  {
1871  case ePhysInterpToEquiSpaced:
1872  {
1873  int nq0 = m_base[0]->GetNumPoints();
1874  int nq1 = m_base[1]->GetNumPoints();
1875  int nq2 = m_base[2]->GetNumPoints();
1876  int nq;
1877 
1878  // take definition from key
1879  if(mkey.ConstFactorExists(eFactorConst))
1880  {
1881  nq = (int) mkey.GetConstFactor(eFactorConst);
1882  }
1883  else
1884  {
1885  nq = max(nq0,max(nq1,nq2));
1886  }
1887 
1888  int neq = LibUtilities::StdTetData::
1889  getNumberOfCoefficients(nq,nq,nq);
1890  Array<OneD, Array<OneD, NekDouble> > coords(neq);
1891  Array<OneD, NekDouble> coll(3);
1892  Array<OneD, DNekMatSharedPtr> I(3);
1893  Array<OneD, NekDouble> tmp(nq0);
1894 
1895  Mat = MemoryManager<DNekMat>::
1896  AllocateSharedPtr(neq, nq0 * nq1 * nq2);
1897  int cnt = 0;
1898 
1899  for(int i = 0; i < nq; ++i)
1900  {
1901  for(int j = 0; j < nq-i; ++j)
1902  {
1903  for(int k = 0; k < nq-i-j; ++k,++cnt)
1904  {
1905  coords[cnt] = Array<OneD, NekDouble>(3);
1906  coords[cnt][0] = -1.0 + 2*k/(NekDouble)(nq-1);
1907  coords[cnt][1] = -1.0 + 2*j/(NekDouble)(nq-1);
1908  coords[cnt][2] = -1.0 + 2*i/(NekDouble)(nq-1);
1909  }
1910  }
1911  }
1912 
1913  for(int i = 0; i < neq; ++i)
1914  {
1915  LocCoordToLocCollapsed(coords[i],coll);
1916 
1917  I[0] = m_base[0]->GetI(coll);
1918  I[1] = m_base[1]->GetI(coll+1);
1919  I[2] = m_base[2]->GetI(coll+2);
1920 
1921  // interpolate first coordinate direction
1922  NekDouble fac;
1923  for( int k = 0; k < nq2; ++k)
1924  {
1925  for (int j = 0; j < nq1; ++j)
1926  {
1927 
1928  fac = (I[1]->GetPtr())[j]*(I[2]->GetPtr())[k];
1929  Vmath::Smul(nq0,fac,I[0]->GetPtr(),1,tmp,1);
1930 
1931  Vmath::Vcopy(nq0, &tmp[0], 1,
1932  Mat->GetRawPtr()+
1933  k*nq0*nq1*neq+
1934  j*nq0*neq+i,neq);
1935  }
1936  }
1937  }
1938  }
1939  break;
1940  default:
1941  {
1942  Mat = StdExpansion::CreateGeneralMatrix(mkey);
1943  }
1944  break;
1945  }
1946 
1947  return Mat;
1948  }
1949 
1950  DNekMatSharedPtr StdTetExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
1951  {
1952  return v_GenMatrix(mkey);
1953  }
1954 
1955 
1956  //---------------------------------------
1957  // Private helper functions
1958  //---------------------------------------
1959 
1960  /**
1961  * @brief Compute the mode number in the expansion for a particular
1962  * tensorial combination.
1963  *
1964  * Modes are numbered with the r index travelling fastest, followed by
1965  * q and then p, and each q-r plane is of size
1966  * (Q+1)*(Q+2)/2+max(0,R-Q-p)*Q. For example, when P=2, Q=3 and R=4
1967  * the indexing inside each q-r plane (with r increasing upwards and q
1968  * to the right) is:
1969  *
1970  * p = 0: p = 1: p = 2:
1971  * ----------------------------------
1972  * 4
1973  * 3 8 17
1974  * 2 7 11 16 20 26
1975  * 1 6 10 13 15 19 22 25 28
1976  * 0 5 9 12 14 18 21 23 24 27 29
1977  *
1978  * Note that in this element, we must have that \f$ P \leq Q \leq
1979  * R\f$.
1980  */
1981  int StdTetExp::GetMode(const int I, const int J, const int K)
1982  {
1983  const int Q = m_base[1]->GetNumModes();
1984  const int R = m_base[2]->GetNumModes();
1985 
1986  int i,j,q_hat,k_hat;
1987  int cnt = 0;
1988 
1989  // Traverse to q-r plane number I
1990  for (i = 0; i < I; ++i)
1991  {
1992  // Size of triangle part
1993  q_hat = min(Q,R-i);
1994  // Size of rectangle part
1995  k_hat = max(R-Q-i,0);
1996  cnt += q_hat*(q_hat+1)/2 + k_hat*Q;
1997  }
1998 
1999  // Traverse to q column J
2000  q_hat = R-I;
2001  for (j = 0; j < J; ++j)
2002  {
2003  cnt += q_hat;
2004  q_hat--;
2005  }
2006 
2007  // Traverse up stacks to K
2008  cnt += K;
2009 
2010  return cnt;
2011  }
2012 
2013  void StdTetExp::v_MultiplyByStdQuadratureMetric(
2014  const Array<OneD, const NekDouble>& inarray,
2015  Array<OneD, NekDouble>& outarray)
2016  {
2017  int i, j;
2018 
2019  int nquad0 = m_base[0]->GetNumPoints();
2020  int nquad1 = m_base[1]->GetNumPoints();
2021  int nquad2 = m_base[2]->GetNumPoints();
2022 
2023  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
2024  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
2025  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
2026 
2027  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
2028  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
2029 
2030  // multiply by integration constants
2031  for(i = 0; i < nquad1*nquad2; ++i)
2032  {
2033  Vmath::Vmul(nquad0,(NekDouble*)&inarray[0]+i*nquad0,1,
2034  w0.get(),1, &outarray[0]+i*nquad0,1);
2035  }
2036 
2037  switch(m_base[1]->GetPointsType())
2038  {
2039  // (1,0) Jacobi Inner product.
2040  case LibUtilities::eGaussRadauMAlpha1Beta0:
2041  for(j = 0; j < nquad2; ++j)
2042  {
2043  for(i = 0; i < nquad1; ++i)
2044  {
2045  Blas::Dscal(nquad0,0.5*w1[i], &outarray[0]+i*nquad0+
2046  j*nquad0*nquad1,1);
2047  }
2048  }
2049  break;
2050 
2051  default:
2052  for(j = 0; j < nquad2; ++j)
2053  {
2054  for(i = 0; i < nquad1; ++i)
2055  {
2056  Blas::Dscal(nquad0,
2057  0.5*(1-z1[i])*w1[i],
2058  &outarray[0]+i*nquad0 + j*nquad0*nquad1,
2059  1 );
2060  }
2061  }
2062  break;
2063  }
2064 
2065  switch(m_base[2]->GetPointsType())
2066  {
2067  // (2,0) Jacobi inner product.
2068  case LibUtilities::eGaussRadauMAlpha2Beta0:
2069  for(i = 0; i < nquad2; ++i)
2070  {
2071  Blas::Dscal(nquad0*nquad1, 0.25*w2[i],
2072  &outarray[0]+i*nquad0*nquad1, 1);
2073  }
2074  break;
2075  // (1,0) Jacobi inner product.
2076  case LibUtilities::eGaussRadauMAlpha1Beta0:
2077  for(i = 0; i < nquad2; ++i)
2078  {
2079  Blas::Dscal(nquad0*nquad1, 0.25*(1-z2[i])*w2[i],
2080  &outarray[0]+i*nquad0*nquad1, 1);
2081  }
2082  break;
2083  default:
2084  for(i = 0; i < nquad2; ++i)
2085  {
2086  Blas::Dscal(nquad0*nquad1,0.25*(1-z2[i])*(1-z2[i])*w2[i],
2087  &outarray[0]+i*nquad0*nquad1,1);
2088  }
2089  break;
2090  }
2091  }
2092 
2093  void StdTetExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
2094  const StdMatrixKey &mkey)
2095  {
2096  //To do : 1) add a test to ensure 0 \leq SvvCutoff \leq 1.
2097  // 2) check if the transfer function needs an analytical
2098  // Fourier transform.
2099  // 3) if it doesn't : find a transfer function that renders
2100  // the if( cutoff_a ...) useless to reduce computational
2101  // cost.
2102  // 4) add SVVDiffCoef to both models!!
2103 
2104  int qa = m_base[0]->GetNumPoints();
2105  int qb = m_base[1]->GetNumPoints();
2106  int qc = m_base[2]->GetNumPoints();
2107  int nmodes_a = m_base[0]->GetNumModes();
2108  int nmodes_b = m_base[1]->GetNumModes();
2109  int nmodes_c = m_base[2]->GetNumModes();
2110 
2111  // Declare orthogonal basis.
2112  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
2113  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
2114  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
2115 
2116  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
2117  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B,nmodes_b,pb);
2118  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_C,nmodes_c,pc);
2119 
2120  StdTetExp OrthoExp(Ba,Bb,Bc);
2121 
2122 
2123  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2124  int i,j,k,cnt = 0;
2125 
2126  // project onto physical space.
2127  OrthoExp.FwdTrans(array,orthocoeffs);
2128 
2129  if(mkey.ConstFactorExists(eFactorSVVPowerKerDiffCoeff))
2130  {
2131  // Rodrigo's power kernel
2132  NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
2133  NekDouble SvvDiffCoeff =
2134  mkey.GetConstFactor(eFactorSVVPowerKerDiffCoeff)*
2135  mkey.GetConstFactor(eFactorSVVDiffCoeff);
2136 
2137  for(int i = 0; i < nmodes_a; ++i)
2138  {
2139  for(int j = 0; j < nmodes_b-j; ++j)
2140  {
2141  NekDouble fac1 = std::max(
2142  pow((1.0*i)/(nmodes_a-1),cutoff*nmodes_a),
2143  pow((1.0*j)/(nmodes_b-1),cutoff*nmodes_b));
2144 
2145  for(int k = 0; k < nmodes_c-i-j; ++k)
2146  {
2147  NekDouble fac = std::max(fac1,
2148  pow((1.0*k)/(nmodes_c-1),cutoff*nmodes_c));
2149 
2150  orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2151  cnt++;
2152  }
2153  }
2154  }
2155  }
2156  else if(mkey.ConstFactorExists(eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
2157  {
2158  NekDouble SvvDiffCoeff =
2159  mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff)*
2160  mkey.GetConstFactor(eFactorSVVDiffCoeff);
2161 
2162  int max_abc = max(nmodes_a-kSVVDGFiltermodesmin,
2163  nmodes_b-kSVVDGFiltermodesmin);
2164  max_abc = max(max_abc, nmodes_c-kSVVDGFiltermodesmin);
2165  // clamp max_abc
2166  max_abc = max(max_abc,0);
2167  max_abc = min(max_abc,kSVVDGFiltermodesmax-kSVVDGFiltermodesmin);
2168 
2169  for(int i = 0; i < nmodes_a; ++i)
2170  {
2171  for(int j = 0; j < nmodes_b-j; ++j)
2172  {
2173  int maxij = max(i,j);
2174 
2175  for(int k = 0; k < nmodes_c-i-j; ++k)
2176  {
2177  int maxijk = max(maxij,k);
2178  maxijk = min(maxijk,kSVVDGFiltermodesmax-1);
2179 
2180  orthocoeffs[cnt] *= SvvDiffCoeff *
2181  kSVVDGFilter[max_abc][maxijk];
2182  cnt++;
2183  }
2184  }
2185  }
2186  }
2187  else
2188  {
2189 
2190  //SVV filter paramaters (how much added diffusion
2191  //relative to physical one and fraction of modes from
2192  //which you start applying this added diffusion)
2193 
2194  NekDouble SvvDiffCoeff = mkey.GetConstFactor(StdRegions::eFactorSVVDiffCoeff);
2195  NekDouble SVVCutOff = mkey.GetConstFactor(StdRegions::eFactorSVVCutoffRatio);
2196 
2197  //Defining the cut of mode
2198  int cutoff_a = (int) (SVVCutOff*nmodes_a);
2199  int cutoff_b = (int) (SVVCutOff*nmodes_b);
2200  int cutoff_c = (int) (SVVCutOff*nmodes_c);
2201  int nmodes = min(min(nmodes_a,nmodes_b),nmodes_c);
2202  NekDouble cutoff = min(min(cutoff_a,cutoff_b),cutoff_c);
2203  NekDouble epsilon = 1;
2204 
2205 
2206  //------"New" Version August 22nd '13--------------------
2207  for(i = 0; i < nmodes_a; ++i)
2208  {
2209  for(j = 0; j < nmodes_b-i; ++j)
2210  {
2211  for(k = 0; k < nmodes_c-i-j; ++k)
2212  {
2213  if(i + j + k >= cutoff)
2214  {
2215  orthocoeffs[cnt] *= ((SvvDiffCoeff)*exp(-(i+j+k-nmodes)*(i+j+k-nmodes)/((NekDouble)((i+j+k-cutoff+epsilon)*(i+j+k-cutoff+epsilon)))));
2216  }
2217  else
2218  {
2219  orthocoeffs[cnt] *= 0.0;
2220  }
2221  cnt++;
2222  }
2223  }
2224  }
2225  }
2226 
2227  // backward transform to physical space
2228  OrthoExp.BwdTrans(orthocoeffs,array);
2229  }
2230 
2231 
2232  void StdTetExp::v_ReduceOrderCoeffs(
2233  int numMin,
2234  const Array<OneD, const NekDouble> &inarray,
2235  Array<OneD, NekDouble> &outarray)
2236  {
2237  int nquad0 = m_base[0]->GetNumPoints();
2238  int nquad1 = m_base[1]->GetNumPoints();
2239  int nquad2 = m_base[2]->GetNumPoints();
2240  int nqtot = nquad0 * nquad1 * nquad2;
2241  int nmodes0 = m_base[0]->GetNumModes();
2242  int nmodes1 = m_base[1]->GetNumModes();
2243  int nmodes2 = m_base[2]->GetNumModes();
2244  int numMax = nmodes0;
2245 
2246  Array<OneD, NekDouble> coeff (m_ncoeffs);
2247  Array<OneD, NekDouble> coeff_tmp1(m_ncoeffs, 0.0);
2248  Array<OneD, NekDouble> coeff_tmp2(m_ncoeffs, 0.0);
2249  Array<OneD, NekDouble> phys_tmp (nqtot, 0.0);
2250  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
2251 
2252  Vmath::Vcopy(m_ncoeffs,inarray,1,coeff_tmp2,1);
2253 
2254  const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
2255  const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
2256  const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
2257 
2258  LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A,
2259  nmodes0, Pkey0);
2260  LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_B,
2261  nmodes1, Pkey1);
2262  LibUtilities::BasisKey bortho2(LibUtilities::eOrtho_C,
2263  nmodes2, Pkey2);
2264 
2265  Vmath::Zero(m_ncoeffs, coeff_tmp2, 1);
2266 
2267  StdRegions::StdTetExpSharedPtr OrthoTetExp;
2268  OrthoTetExp = MemoryManager<StdRegions::StdTetExp>
2269  ::AllocateSharedPtr(bortho0, bortho1, bortho2);
2270 
2271  BwdTrans(inarray,phys_tmp);
2272  OrthoTetExp->FwdTrans(phys_tmp, coeff);
2273 
2274  Vmath::Zero(m_ncoeffs,outarray,1);
2275 
2276  // filtering
2277  int cnt = 0;
2278  for (int u = 0; u < numMin; ++u)
2279  {
2280  for (int i = 0; i < numMin-u; ++i)
2281  {
2282  Vmath::Vcopy(numMin - u - i, tmp = coeff + cnt, 1,
2283  tmp2 = coeff_tmp1 + cnt, 1);
2284  cnt += numMax - u - i;
2285  }
2286  for (int i = numMin; i < numMax-u; ++i)
2287  {
2288  cnt += numMax - u - i;
2289  }
2290  }
2291 
2292  OrthoTetExp->BwdTrans(coeff_tmp1,phys_tmp);
2293  FwdTrans(phys_tmp, outarray);
2294  }
2295 
2296 
2297  void StdTetExp::v_GetSimplexEquiSpacedConnectivity(
2298  Array<OneD, int> &conn,
2299  bool standard)
2300  {
2301  boost::ignore_unused(standard);
2302 
2303  int np0 = m_base[0]->GetNumPoints();
2304  int np1 = m_base[1]->GetNumPoints();
2305  int np2 = m_base[2]->GetNumPoints();
2306  int np = max(np0,max(np1,np2));
2307 
2308 
2309  conn = Array<OneD, int>(4*(np-1)*(np-1)*(np-1));
2310 
2311  int row = 0;
2312  int rowp1 = 0;
2313  int plane = 0;
2314  int row1 = 0;
2315  int row1p1 = 0;
2316  int planep1= 0;
2317  int cnt = 0;
2318  for(int i = 0; i < np-1; ++i)
2319  {
2320  planep1 += (np-i)*(np-i+1)/2;
2321  row = 0; // current plane row offset
2322  rowp1 = 0; // current plane row plus one offset
2323  row1 = 0; // next plane row offset
2324  row1p1 = 0; // nex plane row plus one offset
2325  for(int j = 0; j < np-i-1; ++j)
2326  {
2327  rowp1 += np-i-j;
2328  row1p1 += np-i-j-1;
2329  for(int k = 0; k < np-i-j-2; ++k)
2330  {
2331  conn[cnt++] = plane + row +k+1;
2332  conn[cnt++] = plane + row +k;
2333  conn[cnt++] = plane + rowp1 +k;
2334  conn[cnt++] = planep1 + row1 +k;
2335 
2336  conn[cnt++] = plane + row +k+1;
2337  conn[cnt++] = plane + rowp1 +k+1;
2338  conn[cnt++] = planep1 + row1 +k+1;
2339  conn[cnt++] = planep1 + row1 +k;
2340 
2341  conn[cnt++] = plane + rowp1 +k+1;
2342  conn[cnt++] = plane + row +k+1;
2343  conn[cnt++] = plane + rowp1 +k;
2344  conn[cnt++] = planep1 + row1 +k;
2345 
2346  conn[cnt++] = planep1 + row1 +k;
2347  conn[cnt++] = planep1 + row1p1+k;
2348  conn[cnt++] = plane + rowp1 +k;
2349  conn[cnt++] = plane + rowp1 +k+1;
2350 
2351  conn[cnt++] = planep1 + row1 +k;
2352  conn[cnt++] = planep1 + row1p1+k;
2353  conn[cnt++] = planep1 + row1 +k+1;
2354  conn[cnt++] = plane + rowp1 +k+1;
2355 
2356  if(k < np-i-j-3)
2357  {
2358  conn[cnt++] = plane + rowp1 +k+1;
2359  conn[cnt++] = planep1 + row1p1 +k+1;
2360  conn[cnt++] = planep1 + row1 +k+1;
2361  conn[cnt++] = planep1 + row1p1 +k;
2362  }
2363  }
2364 
2365  conn[cnt++] = plane + row + np-i-j-1;
2366  conn[cnt++] = plane + row + np-i-j-2;
2367  conn[cnt++] = plane + rowp1 + np-i-j-2;
2368  conn[cnt++] = planep1 + row1 + np-i-j-2;
2369 
2370  row += np-i-j;
2371  row1 += np-i-j-1;
2372  }
2373  plane += (np-i)*(np-i+1)/2;
2374  }
2375  }
2376 
2377  }//end namespace
2378 }//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:1981
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::StdRegions::StdTetExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Definition: StdTetExp.cpp:1862
Nektar::StdRegions::StdTetExp::v_GetFacePointsKey
virtual LibUtilities::PointsKey v_GetFacePointsKey(const int i, const int j) const
Definition: StdTetExp.cpp:1124
Nektar::StdRegions::StdTetExp::v_GetEdgeNcoeffs
virtual int v_GetEdgeNcoeffs(const int i) const
Definition: StdTetExp.cpp:1015
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::LibUtilities::eModified_C
Principle Modified Functions .
Definition: BasisType.h:50
Nektar::StdRegions::kSVVDGFiltermodesmax
const int kSVVDGFiltermodesmax
Definition: StdRegions.hpp:385
Nektar::StdRegions::StdTetExp::v_GetSimplexEquiSpacedConnectivity
virtual void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true)
Definition: StdTetExp.cpp:2297
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::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:787
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:228
Nektar::StdRegions::StdExpansion::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:945
Nektar::StdRegions::StdTetExp::v_GetNedges
virtual int v_GetNedges() const
Definition: StdTetExp.cpp:958
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:177
Nektar::StdRegions::StdTetExp::v_GetNverts
virtual int v_GetNverts() const
Definition: StdTetExp.cpp:953
Nektar::StdRegions::StdTetExp::v_GetTotalEdgeIntNcoeffs
virtual int v_GetTotalEdgeIntNcoeffs() const
Definition: StdTetExp.cpp:1036
Nektar::StdRegions::StdTetExp::v_GetNfaces
virtual int v_GetNfaces() const
Definition: StdTetExp.cpp:963
Nektar::StdRegions::StdTetExp::v_DetShapeType
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdTetExp.cpp:968
Nektar::StdRegions::EvaluateTriFaceBasisKey
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Definition: StdExpansion3D.cpp:409
Nektar::LibUtilities::eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
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_GetCoords
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_x, Array< OneD, NekDouble > &coords_y, Array< OneD, NekDouble > &coords_z)
Definition: StdTetExp.cpp:1205
std
STL namespace.
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::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:1184
Nektar::StdRegions::StdTetExp::v_DetFaceBasisKey
virtual const LibUtilities::BasisKey v_DetFaceBasisKey(const int i, const int k) const
Definition: StdTetExp.cpp:1157
Nektar::StdRegions::StdTetExp::DetShapeType
LibUtilities::ShapeType DetShapeType() const
Definition: StdTetExp.h:69
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::kSVVDGFiltermodesmin
const int kSVVDGFiltermodesmin
Definition: StdRegions.hpp:384
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
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::StdMatrixKey::GetMatrixType
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:81
Nektar::LibUtilities::StdTetData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:194
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 A[m x n], B[n x k], C[m x k].
Definition: Blas.hpp:213
Nektar::StdRegions::StdTetExp::v_CalcNumberOfCoefficients
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdTetExp.cpp:1144
class_topology.P
P
Definition: class_topology.py:290
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:1553
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::v_FwdTrans
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:456
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:714
Nektar::StdRegions::StdTetExpSharedPtr
std::shared_ptr< StdTetExp > StdTetExpSharedPtr
Definition: StdTetExp.h:279
Nektar::LibUtilities::eGaussRadauMAlpha1Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
Nektar::LibUtilities::eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:46
Nektar::StdRegions::StdTetExp::v_FillMode
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:905
Nektar::StdRegions::StdTetExp::v_NumBndryCoeffs
virtual int v_NumBndryCoeffs() const
Definition: StdTetExp.cpp:973
Nektar::StdRegions::StdTetExp::v_GetFaceNumPoints
virtual int v_GetFaceNumPoints(const int i) const
Definition: StdTetExp.cpp:1103
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:68
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::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:50
Nektar::LibUtilities::BasisKey::GetNumModes
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:86
Nektar::StdRegions::StdTetExp::v_GetTotalFaceIntNcoeffs
virtual int v_GetTotalFaceIntNcoeffs() const
Definition: StdTetExp.cpp:1092
Nektar::StdRegions::StdTetExp::v_GetEdgeToElementMap
virtual void v_GetEdgeToElementMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int P)
Definition: StdTetExp.cpp:1243
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:216
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_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::v_NumDGBndryCoeffs
virtual int v_NumDGBndryCoeffs() const
Definition: StdTetExp.cpp:993
Nektar::StdRegions::StdTetExp::v_GetBoundaryMap
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdTetExp.cpp:1803
Nektar::LibUtilities::eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::StdRegions::StdTetExp::v_GetFaceNumModes
virtual void v_GetFaceNumModes(const int fid, const Orientation faceOrient, int &numModes0, int &numModes1)
Definition: StdTetExp.cpp:914
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_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdTetExp.cpp:1950
Nektar::StdRegions::StdTetExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdTetExp.cpp:1764
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:161
Nektar::LibUtilities::eOrtho_C
Principle Orthogonal Functions .
Definition: BasisType.h:47
Nektar::StdRegions::StdTetExp::v_MultiplyByStdQuadratureMetric
virtual void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:2013
Nektar::StdRegions::StdExpansion::CreateGeneralMatrix
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
Definition: StdExpansion.cpp:324
Nektar::LibUtilities::eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:45
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:59
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
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:168
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:318
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:1301
Nektar::StdRegions::StdTetExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdTetExp.cpp:2093
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:1670
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:121
Nektar::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:130
Nektar::StdRegions::kSVVDGFilter
const NekDouble kSVVDGFilter[9][11]
Definition: StdRegions.hpp:387
Nektar::StdRegions::StdTetExp::v_IsBoundaryInteriorExpansion
virtual bool v_IsBoundaryInteriorExpansion()
Definition: StdTetExp.cpp:1231
Nektar::StdRegions::StdTetExp::v_GetEdgeBasisType
virtual LibUtilities::BasisType v_GetEdgeBasisType(const int i) const
Definition: StdTetExp.cpp:1187
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_IProductWRTBase_MatOp
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:530
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:125
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:140
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:215
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:130
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:274
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:164
Nektar::LibUtilities::eGaussRadauMAlpha2Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:59
Nektar::StdRegions::StdTetExp::v_GetVertexMap
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdTetExp.cpp:1475
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::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:530
Nektar::StdRegions::StdTetExp::v_GetFaceIntNcoeffs
virtual int v_GetFaceIntNcoeffs(const int i) const
Definition: StdTetExp.cpp:1071
Nektar::LibUtilities::NullBasisKey
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:140
Nektar::LibUtilities::eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:54
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:376
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1429
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:110
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Nektar::StdRegions::StdTetExp::v_ReduceOrderCoeffs
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:2232
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:412
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:49
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_GetFaceNcoeffs
virtual int v_GetFaceNcoeffs(const int i) const
Definition: StdTetExp.cpp:1045
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:302
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:186
Nektar::LibUtilities::StdTetData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:217
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:1957
Generated on Wed Jan 20 2021 03:33:42 for Nektar++ by   doxygen 1.8.13