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