Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StdTetExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdTetExp.h
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Header field for tetrahedral routines built upon
33 // StdExpansion3D
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 
38 
39 #include <StdRegions/StdTetExp.h>
40 
41 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  {
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 
910  //---------------------------
911  // Helper functions
912  //---------------------------
913 
914  int StdTetExp::v_GetNverts() const
915  {
916  return 4;
917  }
918 
919  int StdTetExp::v_GetNedges() const
920  {
921  return 6;
922  }
923 
924  int StdTetExp::v_GetNfaces() const
925  {
926  return 4;
927  }
928 
929  LibUtilities::ShapeType StdTetExp::v_DetShapeType() const
930  {
931  return DetShapeType();
932  }
933 
934  int StdTetExp::v_NumBndryCoeffs() const
935  {
936  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
937  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
938  "BasisType is not a boundary interior form");
939  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
940  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
941  "BasisType is not a boundary interior form");
942  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
943  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
944  "BasisType is not a boundary interior form");
945 
946  int P = m_base[0]->GetNumModes();
947  int Q = m_base[1]->GetNumModes();
948  int R = m_base[2]->GetNumModes();
949 
950  return LibUtilities::StdTetData::
951  getNumberOfBndCoefficients(P, Q, R);
952  }
953 
954  int StdTetExp::v_NumDGBndryCoeffs() const
955  {
956  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
957  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
958  "BasisType is not a boundary interior form");
959  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
960  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
961  "BasisType is not a boundary interior form");
962  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
963  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
964  "BasisType is not a boundary interior form");
965 
966  int P = m_base[0]->GetNumModes()-1;
967  int Q = m_base[1]->GetNumModes()-1;
968  int R = m_base[2]->GetNumModes()-1;
969 
970 
971  return (Q+1) + P*(1 + 2*Q - P)/2 // base face
972  + (R+1) + P*(1 + 2*R - P)/2 // front face
973  + 2*(R+1) + Q*(1 + 2*R - Q); // back two faces
974  }
975 
976  int StdTetExp::v_GetEdgeNcoeffs(const int i) const
977  {
978  ASSERTL2((i >= 0) && (i <= 5), "edge id is out of range");
979  int P = m_base[0]->GetNumModes();
980  int Q = m_base[1]->GetNumModes();
981  int R = m_base[2]->GetNumModes();
982 
983  if (i == 0)
984  {
985  return P;
986  }
987  else if (i == 1 || i == 2)
988  {
989  return Q;
990  }
991  else
992  {
993  return R;
994  }
995  }
996 
997  int StdTetExp::v_GetTotalEdgeIntNcoeffs() const
998  {
999  int P = m_base[0]->GetNumModes()-2;
1000  int Q = m_base[1]->GetNumModes()-2;
1001  int R = m_base[2]->GetNumModes()-2;
1002 
1003  return P+Q+4*R;
1004  }
1005 
1006  int StdTetExp::v_GetFaceNcoeffs(const int i) const
1007  {
1008  ASSERTL2((i >= 0) && (i <= 3), "face id is out of range");
1009  int nFaceCoeffs = 0;
1010  int nummodesA, nummodesB, P, Q;
1011  if (i == 0)
1012  {
1013  nummodesA = GetBasisNumModes(0);
1014  nummodesB = GetBasisNumModes(1);
1015  }
1016  else if ((i == 1) || (i == 2))
1017  {
1018  nummodesA = GetBasisNumModes(0);
1019  nummodesB = GetBasisNumModes(2);
1020  }
1021  else
1022  {
1023  nummodesA = GetBasisNumModes(1);
1024  nummodesB = GetBasisNumModes(2);
1025  }
1026  P = nummodesA - 1;
1027  Q = nummodesB - 1;
1028  nFaceCoeffs = Q+1 + (P*(1 + 2*Q - P))/2;
1029  return nFaceCoeffs;
1030  }
1031 
1032  int StdTetExp::v_GetFaceIntNcoeffs(const int i) const
1033  {
1034  ASSERTL2((i >= 0) && (i <= 3), "face id is out of range");
1035  int Pi = m_base[0]->GetNumModes() - 2;
1036  int Qi = m_base[1]->GetNumModes() - 2;
1037  int Ri = m_base[2]->GetNumModes() - 2;
1038 
1039  if((i == 0))
1040  {
1041  return Pi * (2*Qi - Pi - 1) / 2;
1042  }
1043  else if((i == 1))
1044  {
1045  return Pi * (2*Ri - Pi - 1) / 2;
1046  }
1047  else
1048  {
1049  return Qi * (2*Ri - Qi - 1) / 2;
1050  }
1051  }
1052 
1053  int StdTetExp::v_GetTotalFaceIntNcoeffs() const
1054  {
1055  int Pi = m_base[0]->GetNumModes() - 2;
1056  int Qi = m_base[1]->GetNumModes() - 2;
1057  int Ri = m_base[2]->GetNumModes() - 2;
1058 
1059  return Pi * (2*Qi - Pi - 1) / 2 +
1060  Pi * (2*Ri - Pi - 1) / 2 +
1061  Qi * (2*Ri - Qi - 1);
1062  }
1063 
1064  int StdTetExp::v_GetFaceNumPoints(const int i) const
1065  {
1066  ASSERTL2(i >= 0 && i <= 3, "face id is out of range");
1067 
1068  if (i == 0)
1069  {
1070  return m_base[0]->GetNumPoints()*
1071  m_base[1]->GetNumPoints();
1072  }
1073  else if (i == 1)
1074  {
1075  return m_base[0]->GetNumPoints()*
1076  m_base[2]->GetNumPoints();
1077  }
1078  else
1079  {
1080  return m_base[1]->GetNumPoints()*
1081  m_base[2]->GetNumPoints();
1082  }
1083  }
1084 
1085  LibUtilities::PointsKey StdTetExp::v_GetFacePointsKey(
1086  const int i, const int j) const
1087  {
1088  ASSERTL2(i >= 0 && i <= 3, "face id is out of range");
1089  ASSERTL2(j == 0 || j == 1, "face direction is out of range");
1090 
1091  if (i == 0)
1092  {
1093  return m_base[j]->GetPointsKey();
1094  }
1095  else if (i == 1)
1096  {
1097  return m_base[2*j]->GetPointsKey();
1098  }
1099  else
1100  {
1101  return m_base[j+1]->GetPointsKey();
1102  }
1103  }
1104 
1105  int StdTetExp::v_CalcNumberOfCoefficients(
1106  const std::vector<unsigned int>& nummodes,
1107  int & modes_offset)
1108  {
1109  int nmodes = LibUtilities::StdTetData::getNumberOfCoefficients(
1110  nummodes[modes_offset],
1111  nummodes[modes_offset+1],
1112  nummodes[modes_offset+2]);
1113  modes_offset += 3;
1114 
1115  return nmodes;
1116  }
1117 
1118  const LibUtilities::BasisKey StdTetExp::v_DetFaceBasisKey(
1119  const int i, const int k) const
1120  {
1121  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
1122  ASSERTL2(k == 0 || k == 1, "face direction out of range");
1123 
1124  int dir = k;
1125  switch(i)
1126  {
1127  case 0:
1128  dir = k;
1129  break;
1130  case 1:
1131  dir = 2*k;
1132  break;
1133  case 2:
1134  case 3:
1135  dir = k+1;
1136  break;
1137  }
1138 
1139  return EvaluateTriFaceBasisKey(k,
1140  m_base[dir]->GetBasisType(),
1141  m_base[dir]->GetNumPoints(),
1142  m_base[dir]->GetNumModes());
1143 
1144  // Should not get here.
1145  return LibUtilities::NullBasisKey;
1146  }
1147 
1148  LibUtilities::BasisType StdTetExp::v_GetEdgeBasisType(const int i) const
1149  {
1150  ASSERTL2(i >= 0 && i <= 5, "edge id is out of range");
1151 
1152  if (i == 0)
1153  {
1154  return GetBasisType(0);
1155  }
1156  else if (i == 1 || i == 2)
1157  {
1158  return GetBasisType(1);
1159  }
1160  else
1161  {
1162  return GetBasisType(2);
1163  }
1164  }
1165 
1166  void StdTetExp::v_GetCoords(
1167  Array<OneD, NekDouble> &xi_x,
1168  Array<OneD, NekDouble> &xi_y,
1169  Array<OneD, NekDouble> &xi_z)
1170  {
1171  Array<OneD, const NekDouble> eta_x = m_base[0]->GetZ();
1172  Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
1173  Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
1174  int Qx = GetNumPoints(0);
1175  int Qy = GetNumPoints(1);
1176  int Qz = GetNumPoints(2);
1177 
1178  // Convert collapsed coordinates into cartesian coordinates: eta
1179  // --> xi
1180  for( int k = 0; k < Qz; ++k ) {
1181  for( int j = 0; j < Qy; ++j ) {
1182  for( int i = 0; i < Qx; ++i ) {
1183  int s = i + Qx*(j + Qy*k);
1184  xi_x[s] = (eta_x[i] + 1.0) * (1.0 - eta_y[j]) * (1.0 - eta_z[k]) / 4 - 1.0;
1185  xi_y[s] = (eta_y[j] + 1.0) * (1.0 - eta_z[k]) / 2 - 1.0;
1186  xi_z[s] = eta_z[k];
1187  }
1188  }
1189  }
1190  }
1191 
1192  bool StdTetExp::v_IsBoundaryInteriorExpansion()
1193  {
1194  return (m_base[0]->GetBasisType() == LibUtilities::eModified_A) &&
1195  (m_base[1]->GetBasisType() == LibUtilities::eModified_B) &&
1196  (m_base[2]->GetBasisType() == LibUtilities::eModified_C);
1197  }
1198 
1199 
1200  //--------------------------
1201  // Mappings
1202  //--------------------------
1203 
1204  /**
1205  * Maps Expansion2D modes of a 2D face to the corresponding expansion
1206  * modes.
1207  */
1208  void StdTetExp::v_GetFaceToElementMap(
1209  const int fid,
1210  const Orientation faceOrient,
1211  Array<OneD, unsigned int> &maparray,
1212  Array<OneD, int> &signarray,
1213  int P,
1214  int Q)
1215  {
1216  int nummodesA,nummodesB, i, j, k, idx;
1217 
1218  ASSERTL1(v_IsBoundaryInteriorExpansion(),
1219  "Method only implemented for Modified_A BasisType (x "
1220  "direction), Modified_B BasisType (y direction), and "
1221  "Modified_C BasisType(z direction)");
1222 
1223  int nFaceCoeffs = 0;
1224 
1225  switch(fid)
1226  {
1227  case 0:
1228  nummodesA = m_base[0]->GetNumModes();
1229  nummodesB = m_base[1]->GetNumModes();
1230  break;
1231  case 1:
1232  nummodesA = m_base[0]->GetNumModes();
1233  nummodesB = m_base[2]->GetNumModes();
1234  break;
1235  case 2:
1236  case 3:
1237  nummodesA = m_base[1]->GetNumModes();
1238  nummodesB = m_base[2]->GetNumModes();
1239  break;
1240  }
1241 
1242  bool CheckForZeroedModes = false;
1243  if(P == -1)
1244  {
1245  P = nummodesA;
1246  Q = nummodesB;
1247  }
1248  else
1249  {
1250  CheckForZeroedModes = true;
1251  }
1252 
1253  nFaceCoeffs = P*(2*Q-P+1)/2;
1254 
1255  // Allocate the map array and sign array; set sign array to ones (+)
1256  if(maparray.num_elements() != nFaceCoeffs)
1257  {
1258  maparray = Array<OneD, unsigned int>(nFaceCoeffs,1);
1259  }
1260 
1261  if(signarray.num_elements() != nFaceCoeffs)
1262  {
1263  signarray = Array<OneD, int>(nFaceCoeffs,1);
1264  }
1265  else
1266  {
1267  fill(signarray.get(),signarray.get()+nFaceCoeffs, 1 );
1268  }
1269 
1270  switch (fid)
1271  {
1272  case 0:
1273  idx = 0;
1274  for (i = 0; i < P; ++i)
1275  {
1276  for (j = 0; j < Q-i; ++j)
1277  {
1278  if ((int)faceOrient == 7 && i > 1)
1279  {
1280  signarray[idx] = (i%2 ? -1 : 1);
1281  }
1282  maparray[idx++] = GetMode(i,j,0);
1283  }
1284  }
1285  break;
1286  case 1:
1287  idx = 0;
1288  for (i = 0; i < P; ++i)
1289  {
1290  for (k = 0; k < Q-i; ++k)
1291  {
1292  if ((int)faceOrient == 7 && i > 1)
1293  {
1294  signarray[idx] = (i%2 ? -1: 1);
1295  }
1296  maparray[idx++] = GetMode(i,0,k);
1297  }
1298  }
1299  break;
1300  case 2:
1301  idx = 0;
1302  for (j = 0; j < P-1; ++j)
1303  {
1304  for (k = 0; k < Q-1-j; ++k)
1305  {
1306  if ((int)faceOrient == 7 && j > 1)
1307  {
1308  signarray[idx] = ((j+1)%2 ? -1: 1);
1309  }
1310  maparray[idx++] = GetMode(1,j,k);
1311  // Incorporate modes from zeroth plane where needed.
1312  if (j == 0 && k == 0)
1313  {
1314  maparray[idx++] = GetMode(0,0,1);
1315  }
1316  if (j == 0 && k == Q-2)
1317  {
1318  for (int r = 0; r < Q-1; ++r)
1319  {
1320  maparray[idx++] = GetMode(0,1,r);
1321  }
1322  }
1323  }
1324  }
1325  break;
1326  case 3:
1327  idx = 0;
1328  for (j = 0; j < P; ++j)
1329  {
1330  for (k = 0; k < Q-j; ++k)
1331  {
1332  if ((int)faceOrient == 7 && j > 1)
1333  {
1334  signarray[idx] = (j%2 ? -1: 1);
1335  }
1336  maparray[idx++] = GetMode(0,j,k);
1337  }
1338  }
1339  break;
1340  default:
1341  ASSERTL0(false, "Element map not available.");
1342  }
1343 
1344  if ((int)faceOrient == 7)
1345  {
1346  swap(maparray[0], maparray[Q]);
1347 
1348  for (i = 1; i < Q-1; ++i)
1349  {
1350  swap(maparray[i+1], maparray[Q+i]);
1351  }
1352  }
1353 
1354  if(CheckForZeroedModes)
1355  {
1356  // zero signmap and set maparray to zero if elemental
1357  // modes are not as large as face modesl
1358  idx = 0;
1359  for (j = 0; j < nummodesA; ++j)
1360  {
1361  idx += nummodesB-j;
1362  for (k = nummodesB-j; k < Q-j; ++k)
1363  {
1364  signarray[idx] = 0.0;
1365  maparray[idx++] = maparray[0];
1366  }
1367  }
1368 
1369  for (j = nummodesA; j < P; ++j)
1370  {
1371  for (k = 0; k < Q-j; ++k)
1372  {
1373  signarray[idx] = 0.0;
1374  maparray[idx++] = maparray[0];
1375  }
1376  }
1377  }
1378  }
1379 
1380  int StdTetExp::v_GetVertexMap(const int localVertexId, bool useCoeffPacking)
1381  {
1382  ASSERTL0((GetEdgeBasisType(localVertexId)==LibUtilities::eModified_A)||
1383  (GetEdgeBasisType(localVertexId)==LibUtilities::eModified_B)||
1384  (GetEdgeBasisType(localVertexId)==LibUtilities::eModified_C),
1385  "Mapping not defined for this type of basis");
1386 
1387  int localDOF = 0;
1388  if(useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1389  {
1390  switch(localVertexId)
1391  {
1392  case 0:
1393  {
1394  localDOF = GetMode(0,0,0);
1395  break;
1396  }
1397  case 1:
1398  {
1399  localDOF = GetMode(0,0,1);
1400  break;
1401  }
1402  case 2:
1403  {
1404  localDOF = GetMode(0,1,0);
1405  break;
1406  }
1407  case 3:
1408  {
1409  localDOF = GetMode(1,0,0);
1410  break;
1411  }
1412  default:
1413  {
1414  ASSERTL0(false,"Vertex ID must be between 0 and 3");
1415  break;
1416  }
1417  }
1418  }
1419  else
1420  {
1421  switch(localVertexId)
1422  {
1423  case 0:
1424  {
1425  localDOF = GetMode(0,0,0);
1426  break;
1427  }
1428  case 1:
1429  {
1430  localDOF = GetMode(1,0,0);
1431  break;
1432  }
1433  case 2:
1434  {
1435  localDOF = GetMode(0,1,0);
1436  break;
1437  }
1438  case 3:
1439  {
1440  localDOF = GetMode(0,0,1);
1441  break;
1442  }
1443  default:
1444  {
1445  ASSERTL0(false,"Vertex ID must be between 0 and 3");
1446  break;
1447  }
1448  }
1449 
1450  }
1451 
1452  return localDOF;
1453  }
1454 
1455  /**
1456  * Maps interior modes of an edge to the elemental modes.
1457  */
1458  void StdTetExp::v_GetEdgeInteriorMap(
1459  const int eid,
1460  const Orientation edgeOrient,
1461  Array<OneD, unsigned int> &maparray,
1462  Array<OneD, int> &signarray)
1463  {
1464  int i;
1465  const int P = m_base[0]->GetNumModes();
1466  const int Q = m_base[1]->GetNumModes();
1467  const int R = m_base[2]->GetNumModes();
1468 
1469  const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid)-2;
1470 
1471  if(maparray.num_elements() != nEdgeIntCoeffs)
1472  {
1473  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1474  }
1475  else
1476  {
1477  fill( maparray.get(), maparray.get() + nEdgeIntCoeffs, 0);
1478  }
1479 
1480  if(signarray.num_elements() != nEdgeIntCoeffs)
1481  {
1482  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1483  }
1484  else
1485  {
1486  fill( signarray.get() , signarray.get()+nEdgeIntCoeffs, 1 );
1487  }
1488 
1489  switch (eid)
1490  {
1491  case 0:
1492  for (i = 0; i < P-2; ++i)
1493  {
1494  maparray[i] = GetMode(i+2, 0, 0);
1495  }
1496  if(edgeOrient==eBackwards)
1497  {
1498  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1499  {
1500  signarray[i] = -1;
1501  }
1502  }
1503  break;
1504  case 1:
1505  for (i = 0; i < Q-2; ++i)
1506  {
1507  maparray[i] = GetMode(1, i+1, 0);
1508  }
1509  if(edgeOrient==eBackwards)
1510  {
1511  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1512  {
1513  signarray[i] = -1;
1514  }
1515  }
1516  break;
1517  case 2:
1518  for (i = 0; i < Q-2; ++i)
1519  {
1520  maparray[i] = GetMode(0, i+2, 0);
1521  }
1522  if(edgeOrient==eBackwards)
1523  {
1524  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1525  {
1526  signarray[i] = -1;
1527  }
1528  }
1529  break;
1530  case 3:
1531  for (i = 0; i < R-2; ++i)
1532  {
1533  maparray[i] = GetMode(0, 0, i+2);
1534  }
1535  if(edgeOrient==eBackwards)
1536  {
1537  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1538  {
1539  signarray[i] = -1;
1540  }
1541  }
1542  break;
1543  case 4:
1544  for (i = 0; i < R - 2; ++i)
1545  {
1546  maparray[i] = GetMode(1, 0, i+1);
1547  }
1548  if(edgeOrient==eBackwards)
1549  {
1550  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1551  {
1552  signarray[i] = -1;
1553  }
1554  }
1555  break;
1556  case 5:
1557  for (i = 0; i < R-2; ++i)
1558  {
1559  maparray[i] = GetMode(0, 1, i+1);
1560  }
1561  if(edgeOrient==eBackwards)
1562  {
1563  for(i = 1; i < nEdgeIntCoeffs; i+=2)
1564  {
1565  signarray[i] = -1;
1566  }
1567  }
1568  break;
1569  default:
1570  ASSERTL0(false, "Edge not defined.");
1571  break;
1572  }
1573  }
1574 
1575  void StdTetExp::v_GetFaceInteriorMap(
1576  const int fid,
1577  const Orientation faceOrient,
1578  Array<OneD, unsigned int> &maparray,
1579  Array<OneD, int> &signarray)
1580  {
1581  int i, j, idx, k;
1582  const int P = m_base[0]->GetNumModes();
1583  const int Q = m_base[1]->GetNumModes();
1584  const int R = m_base[2]->GetNumModes();
1585 
1586  const int nFaceIntCoeffs = v_GetFaceIntNcoeffs(fid);
1587 
1588  if(maparray.num_elements() != nFaceIntCoeffs)
1589  {
1590  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1591  }
1592 
1593  if(signarray.num_elements() != nFaceIntCoeffs)
1594  {
1595  signarray = Array<OneD, int>(nFaceIntCoeffs,1);
1596  }
1597  else
1598  {
1599  fill( signarray.get() , signarray.get()+nFaceIntCoeffs, 1 );
1600  }
1601 
1602  switch (fid)
1603  {
1604  case 0:
1605  idx = 0;
1606  for (i = 2; i < P-1; ++i)
1607  {
1608  for (j = 1; j < Q-i; ++j)
1609  {
1610  if ((int)faceOrient == 7)
1611  {
1612  signarray[idx] = (i%2 ? -1 : 1);
1613  }
1614  maparray[idx++] = GetMode(i,j,0);
1615  }
1616  }
1617  break;
1618  case 1:
1619  idx = 0;
1620  for (i = 2; i < P; ++i)
1621  {
1622  for (k = 1; k < R-i; ++k)
1623  {
1624  if ((int)faceOrient == 7)
1625  {
1626  signarray[idx] = (i%2 ? -1: 1);
1627  }
1628  maparray[idx++] = GetMode(i,0,k);
1629  }
1630  }
1631  break;
1632  case 2:
1633  idx = 0;
1634  for (j = 1; j < Q-2; ++j)
1635  {
1636  for (k = 1; k < R-1-j; ++k)
1637  {
1638  if ((int)faceOrient == 7)
1639  {
1640  signarray[idx] = ((j+1)%2 ? -1: 1);
1641  }
1642  maparray[idx++] = GetMode(1,j,k);
1643  }
1644  }
1645  break;
1646  case 3:
1647  idx = 0;
1648  for (j = 2; j < Q-1; ++j)
1649  {
1650  for (k = 1; k < R-j; ++k)
1651  {
1652  if ((int)faceOrient == 7)
1653  {
1654  signarray[idx] = (j%2 ? -1: 1);
1655  }
1656  maparray[idx++] = GetMode(0,j,k);
1657  }
1658  }
1659  break;
1660  default:
1661  ASSERTL0(false, "Face interior map not available.");
1662  break;
1663  }
1664  }
1665 
1666  /**
1667  * List of all interior modes in the expansion.
1668  */
1669  void StdTetExp::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
1670  {
1671  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1672  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1673  "BasisType is not a boundary interior form");
1674  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
1675  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1676  "BasisType is not a boundary interior form");
1677  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1678  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1679  "BasisType is not a boundary interior form");
1680 
1681  int P = m_base[0]->GetNumModes();
1682  int Q = m_base[1]->GetNumModes();
1683  int R = m_base[2]->GetNumModes();
1684 
1685  int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1686 
1687  if(outarray.num_elements() != nIntCoeffs)
1688  {
1689  outarray = Array<OneD, unsigned int>(nIntCoeffs);
1690  }
1691 
1692  int idx = 0;
1693  for (int i = 2; i < P-2; ++i)
1694  {
1695  for (int j = 1; j < Q-i-1; ++j)
1696  {
1697  for (int k = 1; k < R-i-j; ++k)
1698  {
1699  outarray[idx++] = GetMode(i,j,k);
1700  }
1701  }
1702  }
1703  }
1704 
1705  /**
1706  * List of all boundary modes in the the expansion.
1707  */
1708  void StdTetExp::v_GetBoundaryMap(Array<OneD, unsigned int>& outarray)
1709  {
1710  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1711  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1712  "BasisType is not a boundary interior form");
1713  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_B ||
1714  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1715  "BasisType is not a boundary interior form");
1716  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1717  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1718  "BasisType is not a boundary interior form");
1719 
1720  int P = m_base[0]->GetNumModes();
1721  int Q = m_base[1]->GetNumModes();
1722  int R = m_base[2]->GetNumModes();
1723 
1724  int i,j,k;
1725  int idx = 0;
1726 
1727  int nBnd = NumBndryCoeffs();
1728 
1729  if (outarray.num_elements() != nBnd)
1730  {
1731  outarray = Array<OneD, unsigned int>(nBnd);
1732  }
1733 
1734  for (i = 0; i < P; ++i)
1735  {
1736  // First two Q-R planes are entirely boundary modes
1737  if (i < 2)
1738  {
1739  for (j = 0; j < Q-i; j++)
1740  {
1741  for (k = 0; k < R-i-j; ++k)
1742  {
1743  outarray[idx++] = GetMode(i,j,k);
1744  }
1745  }
1746  }
1747  // Remaining Q-R planes contain boundary modes on bottom and
1748  // left edge.
1749  else
1750  {
1751  for (k = 0; k < R-i; ++k)
1752  {
1753  outarray[idx++] = GetMode(i,0,k);
1754  }
1755  for (j = 1; j < Q-i; ++j)
1756  {
1757  outarray[idx++] = GetMode(i,j,0);
1758  }
1759  }
1760  }
1761  }
1762 
1763 
1764  //---------------------------------------
1765  // Wrapper functions
1766  //---------------------------------------
1767  DNekMatSharedPtr StdTetExp::v_GenMatrix(const StdMatrixKey &mkey)
1768  {
1769 
1770  MatrixType mtype = mkey.GetMatrixType();
1771 
1772  DNekMatSharedPtr Mat;
1773 
1774  switch(mtype)
1775  {
1776  case ePhysInterpToEquiSpaced:
1777  {
1778  int nq0 = m_base[0]->GetNumPoints();
1779  int nq1 = m_base[1]->GetNumPoints();
1780  int nq2 = m_base[2]->GetNumPoints();
1781  int nq;
1782 
1783  // take definition from key
1784  if(mkey.ConstFactorExists(eFactorConst))
1785  {
1786  nq = (int) mkey.GetConstFactor(eFactorConst);
1787  }
1788  else
1789  {
1790  nq = max(nq0,max(nq1,nq2));
1791  }
1792 
1793  int neq = LibUtilities::StdTetData::
1794  getNumberOfCoefficients(nq,nq,nq);
1795  Array<OneD, Array<OneD, NekDouble> > coords(neq);
1796  Array<OneD, NekDouble> coll(3);
1797  Array<OneD, DNekMatSharedPtr> I(3);
1798  Array<OneD, NekDouble> tmp(nq0);
1799 
1800  Mat = MemoryManager<DNekMat>::
1801  AllocateSharedPtr(neq, nq0 * nq1 * nq2);
1802  int cnt = 0;
1803 
1804  for(int i = 0; i < nq; ++i)
1805  {
1806  for(int j = 0; j < nq-i; ++j)
1807  {
1808  for(int k = 0; k < nq-i-j; ++k,++cnt)
1809  {
1810  coords[cnt] = Array<OneD, NekDouble>(3);
1811  coords[cnt][0] = -1.0 + 2*k/(NekDouble)(nq-1);
1812  coords[cnt][1] = -1.0 + 2*j/(NekDouble)(nq-1);
1813  coords[cnt][2] = -1.0 + 2*i/(NekDouble)(nq-1);
1814  }
1815  }
1816  }
1817 
1818  for(int i = 0; i < neq; ++i)
1819  {
1820  LocCoordToLocCollapsed(coords[i],coll);
1821 
1822  I[0] = m_base[0]->GetI(coll);
1823  I[1] = m_base[1]->GetI(coll+1);
1824  I[2] = m_base[2]->GetI(coll+2);
1825 
1826  // interpolate first coordinate direction
1827  NekDouble fac;
1828  for( int k = 0; k < nq2; ++k)
1829  {
1830  for (int j = 0; j < nq1; ++j)
1831  {
1832 
1833  fac = (I[1]->GetPtr())[j]*(I[2]->GetPtr())[k];
1834  Vmath::Smul(nq0,fac,I[0]->GetPtr(),1,tmp,1);
1835 
1836  Vmath::Vcopy(nq0, &tmp[0], 1,
1837  Mat->GetRawPtr()+
1838  k*nq0*nq1*neq+
1839  j*nq0*neq+i,neq);
1840  }
1841  }
1842  }
1843  }
1844  break;
1845  default:
1846  {
1847  Mat = StdExpansion::CreateGeneralMatrix(mkey);
1848  }
1849  break;
1850  }
1851 
1852  return Mat;
1853  }
1854 
1855  DNekMatSharedPtr StdTetExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
1856  {
1857  return v_GenMatrix(mkey);
1858  }
1859 
1860 
1861  //---------------------------------------
1862  // Private helper functions
1863  //---------------------------------------
1864 
1865  /**
1866  * @brief Compute the mode number in the expansion for a particular
1867  * tensorial combination.
1868  *
1869  * Modes are numbered with the r index travelling fastest, followed by
1870  * q and then p, and each q-r plane is of size
1871  * (Q+1)*(Q+2)/2+max(0,R-Q-p)*Q. For example, when P=2, Q=3 and R=4
1872  * the indexing inside each q-r plane (with r increasing upwards and q
1873  * to the right) is:
1874  *
1875  * p = 0: p = 1: p = 2:
1876  * ----------------------------------
1877  * 4
1878  * 3 8 17
1879  * 2 7 11 16 20 26
1880  * 1 6 10 13 15 19 22 25 28
1881  * 0 5 9 12 14 18 21 23 24 27 29
1882  *
1883  * Note that in this element, we must have that \f$ P \leq Q \leq
1884  * R\f$.
1885  */
1886  int StdTetExp::GetMode(const int I, const int J, const int K)
1887  {
1888  const int Q = m_base[1]->GetNumModes();
1889  const int R = m_base[2]->GetNumModes();
1890 
1891  int i,j,q_hat,k_hat;
1892  int cnt = 0;
1893 
1894  // Traverse to q-r plane number I
1895  for (i = 0; i < I; ++i)
1896  {
1897  // Size of triangle part
1898  q_hat = min(Q,R-i);
1899  // Size of rectangle part
1900  k_hat = max(R-Q-i,0);
1901  cnt += q_hat*(q_hat+1)/2 + k_hat*Q;
1902  }
1903 
1904  // Traverse to q column J
1905  q_hat = R-I;
1906  for (j = 0; j < J; ++j)
1907  {
1908  cnt += q_hat;
1909  q_hat--;
1910  }
1911 
1912  // Traverse up stacks to K
1913  cnt += K;
1914 
1915  return cnt;
1916  }
1917 
1918  void StdTetExp::v_MultiplyByStdQuadratureMetric(
1919  const Array<OneD, const NekDouble>& inarray,
1920  Array<OneD, NekDouble>& outarray)
1921  {
1922  int i, j;
1923 
1924  int nquad0 = m_base[0]->GetNumPoints();
1925  int nquad1 = m_base[1]->GetNumPoints();
1926  int nquad2 = m_base[2]->GetNumPoints();
1927 
1928  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1929  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1930  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
1931 
1932  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
1933  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1934 
1935  // multiply by integration constants
1936  for(i = 0; i < nquad1*nquad2; ++i)
1937  {
1938  Vmath::Vmul(nquad0,(NekDouble*)&inarray[0]+i*nquad0,1,
1939  w0.get(),1, &outarray[0]+i*nquad0,1);
1940  }
1941 
1942  switch(m_base[1]->GetPointsType())
1943  {
1944  // (1,0) Jacobi Inner product.
1945  case LibUtilities::eGaussRadauMAlpha1Beta0:
1946  for(j = 0; j < nquad2; ++j)
1947  {
1948  for(i = 0; i < nquad1; ++i)
1949  {
1950  Blas::Dscal(nquad0,0.5*w1[i], &outarray[0]+i*nquad0+
1951  j*nquad0*nquad1,1);
1952  }
1953  }
1954  break;
1955 
1956  default:
1957  for(j = 0; j < nquad2; ++j)
1958  {
1959  for(i = 0; i < nquad1; ++i)
1960  {
1961  Blas::Dscal(nquad0,
1962  0.5*(1-z1[i])*w1[i],
1963  &outarray[0]+i*nquad0 + j*nquad0*nquad1,
1964  1 );
1965  }
1966  }
1967  break;
1968  }
1969 
1970  switch(m_base[2]->GetPointsType())
1971  {
1972  // (2,0) Jacobi inner product.
1973  case LibUtilities::eGaussRadauMAlpha2Beta0:
1974  for(i = 0; i < nquad2; ++i)
1975  {
1976  Blas::Dscal(nquad0*nquad1, 0.25*w2[i],
1977  &outarray[0]+i*nquad0*nquad1, 1);
1978  }
1979  break;
1980  // (1,0) Jacobi inner product.
1981  case LibUtilities::eGaussRadauMAlpha1Beta0:
1982  for(i = 0; i < nquad2; ++i)
1983  {
1984  Blas::Dscal(nquad0*nquad1, 0.25*(1-z2[i])*w2[i],
1985  &outarray[0]+i*nquad0*nquad1, 1);
1986  }
1987  break;
1988  default:
1989  for(i = 0; i < nquad2; ++i)
1990  {
1991  Blas::Dscal(nquad0*nquad1,0.25*(1-z2[i])*(1-z2[i])*w2[i],
1992  &outarray[0]+i*nquad0*nquad1,1);
1993  }
1994  break;
1995  }
1996  }
1997 
1998  void StdTetExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
1999  const StdMatrixKey &mkey)
2000  {
2001  //To do : 1) add a test to ensure 0 \leq SvvCutoff \leq 1.
2002  // 2) check if the transfer function needs an analytical
2003  // Fourier transform.
2004  // 3) if it doesn't : find a transfer function that renders
2005  // the if( cutoff_a ...) useless to reduce computational
2006  // cost.
2007  // 4) add SVVDiffCoef to both models!!
2008 
2009  int qa = m_base[0]->GetNumPoints();
2010  int qb = m_base[1]->GetNumPoints();
2011  int qc = m_base[2]->GetNumPoints();
2012  int nmodes_a = m_base[0]->GetNumModes();
2013  int nmodes_b = m_base[1]->GetNumModes();
2014  int nmodes_c = m_base[2]->GetNumModes();
2015 
2016  // Declare orthogonal basis.
2017  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
2018  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
2019  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
2020 
2021  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
2022  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B,nmodes_b,pb);
2023  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_C,nmodes_c,pc);
2024 
2025  StdTetExp OrthoExp(Ba,Bb,Bc);
2026 
2027 
2028  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2029  int i,j,k,cnt = 0;
2030 
2031  //SVV filter paramaters (how much added diffusion relative to physical one
2032  // and fraction of modes from which you start applying this added diffusion)
2033  //
2034  NekDouble SvvDiffCoeff = mkey.GetConstFactor(StdRegions::eFactorSVVDiffCoeff);
2035  NekDouble SVVCutOff = mkey.GetConstFactor(StdRegions::eFactorSVVCutoffRatio);
2036 
2037 
2038  //Defining the cut of mode
2039  int cutoff_a = (int) (SVVCutOff*nmodes_a);
2040  int cutoff_b = (int) (SVVCutOff*nmodes_b);
2041  int cutoff_c = (int) (SVVCutOff*nmodes_c);
2042  int nmodes = min(min(nmodes_a,nmodes_b),nmodes_c);
2043  NekDouble cutoff = min(min(cutoff_a,cutoff_b),cutoff_c);
2044  NekDouble epsilon = 1;
2045 
2046  // project onto physical space.
2047  OrthoExp.FwdTrans(array,orthocoeffs);
2048 
2049  //------"New" Version August 22nd '13--------------------
2050  for(i = 0; i < nmodes_a; ++i)
2051  {
2052  for(j = 0; j < nmodes_b-i; ++j)
2053  {
2054  for(k = 0; k < nmodes_c-i-j; ++k)
2055  {
2056  if(i + j + k >= cutoff)
2057  {
2058  orthocoeffs[cnt] *= ((SvvDiffCoeff)*exp(-(i+j+k-nmodes)*(i+j+k-nmodes)/((NekDouble)((i+j+k-cutoff+epsilon)*(i+j+k-cutoff+epsilon)))));
2059  }
2060  else
2061  {
2062  orthocoeffs[cnt] *= 0.0;
2063  }
2064  cnt++;
2065  }
2066  }
2067  }
2068  // backward transform to physical space
2069  OrthoExp.BwdTrans(orthocoeffs,array);
2070  }
2071 
2072 
2073  void StdTetExp::v_ReduceOrderCoeffs(
2074  int numMin,
2075  const Array<OneD, const NekDouble> &inarray,
2076  Array<OneD, NekDouble> &outarray)
2077  {
2078  int nquad0 = m_base[0]->GetNumPoints();
2079  int nquad1 = m_base[1]->GetNumPoints();
2080  int nquad2 = m_base[2]->GetNumPoints();
2081  int nqtot = nquad0 * nquad1 * nquad2;
2082  int nmodes0 = m_base[0]->GetNumModes();
2083  int nmodes1 = m_base[1]->GetNumModes();
2084  int nmodes2 = m_base[2]->GetNumModes();
2085  int numMax = nmodes0;
2086 
2087  Array<OneD, NekDouble> coeff (m_ncoeffs);
2088  Array<OneD, NekDouble> coeff_tmp1(m_ncoeffs, 0.0);
2089  Array<OneD, NekDouble> coeff_tmp2(m_ncoeffs, 0.0);
2090  Array<OneD, NekDouble> phys_tmp (nqtot, 0.0);
2091  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
2092 
2093  Vmath::Vcopy(m_ncoeffs,inarray,1,coeff_tmp2,1);
2094 
2095  const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
2096  const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
2097  const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
2098 
2099  LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A,
2100  nmodes0, Pkey0);
2101  LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_B,
2102  nmodes1, Pkey1);
2103  LibUtilities::BasisKey bortho2(LibUtilities::eOrtho_C,
2104  nmodes2, Pkey2);
2105 
2106  Vmath::Zero(m_ncoeffs, coeff_tmp2, 1);
2107 
2108  StdRegions::StdTetExpSharedPtr OrthoTetExp;
2109  OrthoTetExp = MemoryManager<StdRegions::StdTetExp>
2110  ::AllocateSharedPtr(bortho0, bortho1, bortho2);
2111 
2112  BwdTrans(inarray,phys_tmp);
2113  OrthoTetExp->FwdTrans(phys_tmp, coeff);
2114 
2115  Vmath::Zero(m_ncoeffs,outarray,1);
2116 
2117  // filtering
2118  int cnt = 0;
2119  for (int u = 0; u < numMin; ++u)
2120  {
2121  for (int i = 0; i < numMin-u; ++i)
2122  {
2123  Vmath::Vcopy(numMin - u - i, tmp = coeff + cnt, 1,
2124  tmp2 = coeff_tmp1 + cnt, 1);
2125  cnt += numMax - u - i;
2126  }
2127  for (int i = numMin; i < numMax-u; ++i)
2128  {
2129  cnt += numMax - u - i;
2130  }
2131  }
2132 
2133  OrthoTetExp->BwdTrans(coeff_tmp1,phys_tmp);
2134  FwdTrans(phys_tmp, outarray);
2135  }
2136 
2137 
2138  void StdTetExp::v_GetSimplexEquiSpacedConnectivity(
2139  Array<OneD, int> &conn,
2140  bool standard)
2141  {
2142  int np0 = m_base[0]->GetNumPoints();
2143  int np1 = m_base[1]->GetNumPoints();
2144  int np2 = m_base[2]->GetNumPoints();
2145  int np = max(np0,max(np1,np2));
2146 
2147 
2148  conn = Array<OneD, int>(4*(np-1)*(np-1)*(np-1));
2149 
2150  int row = 0;
2151  int rowp1 = 0;
2152  int plane = 0;
2153  int row1 = 0;
2154  int row1p1 = 0;
2155  int planep1= 0;
2156  int cnt = 0;
2157  for(int i = 0; i < np-1; ++i)
2158  {
2159  planep1 += (np-i)*(np-i+1)/2;
2160  row = 0; // current plane row offset
2161  rowp1 = 0; // current plane row plus one offset
2162  row1 = 0; // next plane row offset
2163  row1p1 = 0; // nex plane row plus one offset
2164  for(int j = 0; j < np-i-1; ++j)
2165  {
2166  rowp1 += np-i-j;
2167  row1p1 += np-i-j-1;
2168  for(int k = 0; k < np-i-j-2; ++k)
2169  {
2170  conn[cnt++] = plane + row +k+1;
2171  conn[cnt++] = plane + row +k;
2172  conn[cnt++] = plane + rowp1 +k;
2173  conn[cnt++] = planep1 + row1 +k;
2174 
2175  conn[cnt++] = plane + row +k+1;
2176  conn[cnt++] = plane + rowp1 +k+1;
2177  conn[cnt++] = planep1 + row1 +k+1;
2178  conn[cnt++] = planep1 + row1 +k;
2179 
2180  conn[cnt++] = plane + rowp1 +k+1;
2181  conn[cnt++] = plane + row +k+1;
2182  conn[cnt++] = plane + rowp1 +k;
2183  conn[cnt++] = planep1 + row1 +k;
2184 
2185  conn[cnt++] = planep1 + row1 +k;
2186  conn[cnt++] = planep1 + row1p1+k;
2187  conn[cnt++] = plane + rowp1 +k;
2188  conn[cnt++] = plane + rowp1 +k+1;
2189 
2190  conn[cnt++] = planep1 + row1 +k;
2191  conn[cnt++] = planep1 + row1p1+k;
2192  conn[cnt++] = planep1 + row1 +k+1;
2193  conn[cnt++] = plane + rowp1 +k+1;
2194 
2195  if(k < np-i-j-3)
2196  {
2197  conn[cnt++] = plane + rowp1 +k+1;
2198  conn[cnt++] = planep1 + row1p1 +k+1;
2199  conn[cnt++] = planep1 + row1 +k+1;
2200  conn[cnt++] = planep1 + row1p1 +k;
2201  }
2202  }
2203 
2204  conn[cnt++] = plane + row + np-i-j-1;
2205  conn[cnt++] = plane + row + np-i-j-2;
2206  conn[cnt++] = plane + rowp1 + np-i-j-2;
2207  conn[cnt++] = planep1 + row1 + np-i-j-2;
2208 
2209  row += np-i-j;
2210  row1 += np-i-j-1;
2211  }
2212  plane += (np-i)*(np-i+1)/2;
2213  }
2214  }
2215 
2216  }//end namespace
2217 }//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:1886
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::StdRegions::StdTetExp::v_NumBndryCoeffs
virtual int v_NumBndryCoeffs() const
Definition: StdTetExp.cpp:934
Nektar::StdRegions::StdTetExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Definition: StdTetExp.cpp:1767
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:2138
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:1148
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:784
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::StdRegions::StdExpansion::GetBasisNumModes
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
Definition: StdExpansion.h:178
Nektar::StdRegions::StdExpansion::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:942
Nektar::StdRegions::StdTetExp::v_GetEdgeNcoeffs
virtual int v_GetEdgeNcoeffs(const int i) const
Definition: StdTetExp.cpp:976
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:1166
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:1173
Nektar::StdRegions::StdTetExp::v_PhysDeriv
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_dx, Array< OneD, NekDouble > &out_dy, Array< OneD, NekDouble > &out_dz)
Calculate the derivative of the physical points.
Definition: StdTetExp.cpp:108
Nektar::StdRegions::StdTetExp::v_GetNverts
virtual int v_GetNverts() const
Definition: StdTetExp.cpp:914
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:1105
Nektar::DNekMatSharedPtr
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
Nektar::StdRegions::StdTetExp::v_NumDGBndryCoeffs
virtual int v_NumDGBndryCoeffs() const
Definition: StdTetExp.cpp:954
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:1458
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:700
Nektar::StdRegions::StdTetExp::v_GetFaceIntNcoeffs
virtual int v_GetFaceIntNcoeffs(const int i) const
Definition: StdTetExp.cpp:1032
Nektar::StdRegions::StdTetExp::DetShapeType
LibUtilities::ShapeType DetShapeType() const
Definition: StdTetExp.h:70
Nektar::LibUtilities::eGaussRadauMAlpha1Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
Nektar::StdRegions::StdTetExp::v_DetShapeType
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdTetExp.cpp:929
Nektar::StdRegions::StdTetExp::v_GetTotalFaceIntNcoeffs
virtual int v_GetTotalFaceIntNcoeffs() const
Definition: StdTetExp.cpp:1053
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:1085
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
Nektar::StdRegions::StdTetExp::v_GetTotalEdgeIntNcoeffs
virtual int v_GetTotalEdgeIntNcoeffs() const
Definition: StdTetExp.cpp:997
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:1708
Nektar::LibUtilities::eModified_B
Principle Modified Functions .
Definition: BasisType.h:50
Nektar::StdRegions::StdTetExp::v_IProductWRTDerivBase
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:669
Nektar::StdRegions::StdTetExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdTetExp.cpp:1855
Nektar::StdRegions::StdTetExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdTetExp.cpp:1669
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:1918
Nektar::StdRegions::StdExpansion::CreateGeneralMatrix
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
Definition: StdExpansion.cpp:322
Nektar::LibUtilities::eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:46
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:301
Nektar::StdRegions::StdTetExp::v_GetFaceToElementMap
virtual void v_GetFaceToElementMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int P=-1, int Q=-1)
Definition: StdTetExp.cpp:1208
Nektar::StdRegions::StdTetExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdTetExp.cpp:1998
Vmath::Ddot
T Ddot(int n, const Array< OneD, const T > &w, const int incw, const Array< OneD, const T > &x, const int incx, const Array< OneD, const int > &y, const int incy)
Definition: VmathArray.hpp:434
Nektar::StdRegions::StdTetExp::v_GetFaceInteriorMap
virtual void v_GetFaceInteriorMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
Definition: StdTetExp.cpp:1575
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:1192
Nektar::StdRegions::StdTetExp::v_GetNfaces
virtual int v_GetNfaces() const
Definition: StdTetExp.cpp:924
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:240
Nektar::StdRegions::StdTetExp::v_GetFaceNcoeffs
virtual int v_GetFaceNcoeffs(const int i) const
Definition: StdTetExp.cpp:1006
Nektar::LibUtilities::eGaussRadauMAlpha2Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
Nektar::StdRegions::StdTetExp::v_DetFaceBasisKey
virtual const LibUtilities::BasisKey v_DetFaceBasisKey(const int i, const int k) const
Definition: StdTetExp.cpp:1118
Nektar::StdRegions::StdTetExp::v_GetVertexMap
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdTetExp.cpp:1380
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:525
Nektar::LibUtilities::NullBasisKey
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
Nektar::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:131
Nektar::StdRegions::StdTetExp::v_GetFaceNumPoints
virtual int v_GetFaceNumPoints(const int i) const
Definition: StdTetExp.cpp:1064
Nektar::StdRegions::StdTetExpSharedPtr
boost::shared_ptr< StdTetExp > StdTetExpSharedPtr
Definition: StdTetExp.h:268
Nektar::LibUtilities::BasisKey::GetNumModes
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
Nektar::LibUtilities::eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:53
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Nektar::StdRegions::StdExpansion::GetPointsType
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:216
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:218
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1422
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Nektar::StdRegions::StdTetExp::v_GetNedges
virtual int v_GetNedges() const
Definition: StdTetExp.cpp:919
Nektar::StdRegions::StdTetExp::v_ReduceOrderCoeffs
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTetExp.cpp:2073
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:285
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
Nektar::LibUtilities::StdTetData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:209
StdTetExp.h
Nektar::StdRegions::StdExpansion::FwdTrans
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space...
Definition: StdExpansion.h:1905
Generated on Fri Dec 2 2016 13:10:52 for Nektar++ by   doxygen 1.8.8