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