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