Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StdPyrExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdPyrExp.cpp
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: pyramadic routines built upon StdExpansion3D
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <StdRegions/StdPyrExp.h>
37 #include <LibUtilities/Foundations/ManagerAccess.h>
38 #include <iomanip>
39 
40 namespace Nektar
41 {
42  namespace StdRegions
43  {
44  StdPyrExp::StdPyrExp() // Deafult construct of standard expansion directly called.
45  {
46  }
47 
48  StdPyrExp::StdPyrExp(const LibUtilities::BasisKey &Ba,
49  const LibUtilities::BasisKey &Bb,
50  const LibUtilities::BasisKey &Bc)
51  : StdExpansion (LibUtilities::StdPyrData::getNumberOfCoefficients(
52  Ba.GetNumModes(),
53  Bb.GetNumModes(),
54  Bc.GetNumModes()),
55  3, Ba, Bb, Bc),
56  StdExpansion3D(LibUtilities::StdPyrData::getNumberOfCoefficients(
57  Ba.GetNumModes(),
58  Bb.GetNumModes(),
59  Bc.GetNumModes()),
60  Ba, Bb, Bc)
61  {
62  if (Ba.GetNumModes() > Bc.GetNumModes())
63  {
64  ASSERTL0(false, "order in 'a' direction is higher "
65  "than order in 'c' direction");
66  }
67  if (Bb.GetNumModes() > Bc.GetNumModes())
68  {
69  ASSERTL0(false, "order in 'b' direction is higher "
70  "than order in 'c' direction");
71  }
72 
73  // Set up mode mapping which takes 0\leq i\leq N_coeffs -> (p,q,r)
74  // of the 3D tensor product
75  const int P = Ba.GetNumModes() - 1;
76  const int Q = Bb.GetNumModes() - 1;
77  const int R = Bc.GetNumModes() - 1;
78  int cnt = 0;
79 
80  // Vertices
81  m_map[Mode(0, 0, 0, 0)] = cnt++;
82  m_map[Mode(1, 0, 0, 0)] = cnt++;
83  m_map[Mode(1, 1, 0, 0)] = cnt++;
84  m_map[Mode(0, 1, 0, 0)] = cnt++;
85  m_map[Mode(0, 0, 1, 1)] = cnt++;
86 
87  // Edge 0
88  for (int i = 2; i <= P; ++i)
89  {
90  m_map[Mode(i, 0, 0, GetTetMode(i, 0, 0))] = cnt++;
91  }
92 
93  // Edge 1
94  for (int i = 2; i <= Q; ++i)
95  {
96  m_map[Mode(1, i, 0, GetTetMode(0, i, 0))] = cnt++;
97  }
98 
99  // Edge 2
100  for (int i = 2; i <= P; ++i)
101  {
102  m_map[Mode(i, 1, 0, GetTetMode(i, 0, 0))] = cnt++;
103  }
104 
105  // Edge 3
106  for (int i = 2; i <= Q; ++i)
107  {
108  m_map[Mode(0, i, 0, GetTetMode(0, i, 0))] = cnt++;
109  }
110 
111  // Edge 4
112  for (int i = 2; i <= R; ++i)
113  {
114  m_map[Mode(0, 0, i, i)] = cnt++;
115  }
116 
117  // Edge 5
118  for (int i = 2; i <= R; ++i)
119  {
120  m_map[Mode(1, 0, i, i)] = cnt++;
121  }
122 
123  // Edge 6
124  for (int i = 2; i <= R; ++i)
125  {
126  m_map[Mode(1, 1, i, i)] = cnt++;
127  }
128 
129  // Edge 7
130  for (int i = 2; i <= R; ++i)
131  {
132  m_map[Mode(0, 1, i, i)] = cnt++;
133  }
134 
135  // Face 0 - TODO check this
136  for (int j = 2; j <= Q; ++j)
137  {
138  for (int i = 2; i <= P; ++i)
139  {
140  m_map[Mode(i, j, 0, GetTetMode((i-2+j-2) % (Q-1) + 2, 0, 0))] = cnt++;
141  }
142  }
143 
144  // Face 1
145  for (int i = 2; i <= P; ++i)
146  {
147  for (int j = 1; j <= R-i; ++j)
148  {
149  m_map[Mode(i, 0, j, GetTetMode(i, 0, j))] = cnt++;
150  }
151  }
152 
153  // Face 2
154  for (int i = 2; i <= Q; ++i)
155  {
156  for (int j = 1; j <= R-i; ++j)
157  {
158  m_map[Mode(1, i, j, GetTetMode(0, i, j))] = cnt++;
159  }
160  }
161 
162  // Face 3
163  for (int i = 2; i <= P; ++i)
164  {
165  for (int j = 1; j <= R-i; ++j)
166  {
167  m_map[Mode(i, 1, j, GetTetMode(i, 0, j))] = cnt++;
168  }
169  }
170 
171  // Face 4
172  for (int i = 2; i <= Q; ++i)
173  {
174  for (int j = 1; j <= R-i; ++j)
175  {
176  m_map[Mode(0, i, j, GetTetMode(0, i, j))] = cnt++;
177  }
178  }
179 
180  // Interior (tetrahedral modes)
181  for (int i = 2; i <= P+1; ++i)
182  {
183  for (int j = 1; j <= Q-i+1; ++j)
184  {
185  for (int k = 1; k <= R-i-j+1; ++k)
186  {
187  // need to go to j+1-th mode in the 'b' direction to
188  // select correct modified_a mode
189  m_map[Mode(i, j+1, k, GetTetMode(i-1, j, k))] = cnt++;
190  }
191  }
192  }
193 
194  ASSERTL0(m_map.size() == m_ncoeffs,
195  "Duplicate coefficient entries in map");
196 
197  map<Mode, unsigned int, cmpop>::iterator it;
198  for (it = m_map.begin(); it != m_map.end(); ++it)
199  {
200  const int p = it->first.get<0>();
201  const int q = it->first.get<1>();
202  const int r = it->first.get<2>();
203  const int rp = it->first.get<3>();
204  if (m_idxMap.find(p) == m_idxMap.end())
205  {
206  m_idxMap[p] = map<int, map<int, pair<int, int> > >();
207  }
208 
209  if (m_idxMap[p].find(q) == m_idxMap[p].end())
210  {
211  m_idxMap[p][q] = map<int, pair<int, int> >();
212  }
213 
214  if (m_idxMap[p][q].find(r) == m_idxMap[p][q].end())
215  {
216  m_idxMap[p][q][r] = pair<int, int>(it->second, rp);
217  }
218  }
219  }
220 
221  StdPyrExp::StdPyrExp(const StdPyrExp &T)
222  : StdExpansion (T),
223  StdExpansion3D(T)
224  {
225  }
226 
227 
228  // Destructor
229  StdPyrExp::~StdPyrExp()
230  {
231  }
232 
233 
234  //---------------------------------------
235  // Differentiation/integration Methods
236  //---------------------------------------
237 
238  /**
239  * \brief Calculate the derivative of the physical points
240  *
241  * The derivative is evaluated at the nodal physical points.
242  * Derivatives with respect to the local Cartesian coordinates.
243  *
244  * \f$\begin{Bmatrix} \frac {\partial} {\partial \xi_1} \\ \frac
245  * {\partial} {\partial \xi_2} \\ \frac {\partial} {\partial \xi_3}
246  * \end{Bmatrix} = \begin{Bmatrix} \frac 2 {(1-\eta_3)} \frac \partial
247  * {\partial \bar \eta_1} \\ \frac {\partial} {\partial \xi_2} \ \
248  * \frac {(1 + \bar \eta_1)} {(1 - \eta_3)} \frac \partial {\partial
249  * \bar \eta_1} + \frac {\partial} {\partial \eta_3} \end{Bmatrix}\f$
250  */
251  void StdPyrExp::v_PhysDeriv(
252  const Array<OneD, const NekDouble> &u_physical,
253  Array<OneD, NekDouble> &out_dxi1,
254  Array<OneD, NekDouble> &out_dxi2,
255  Array<OneD, NekDouble> &out_dxi3)
256  {
257  // PhysDerivative implementation based on Spen's book page 152.
258  int Qx = m_base[0]->GetNumPoints();
259  int Qy = m_base[1]->GetNumPoints();
260  int Qz = m_base[2]->GetNumPoints();
261 
262  Array<OneD, NekDouble> dEta_bar1(Qx*Qy*Qz,0.0);
263  Array<OneD, NekDouble> dXi2 (Qx*Qy*Qz,0.0);
264  Array<OneD, NekDouble> dEta3 (Qx*Qy*Qz,0.0);
265  PhysTensorDeriv(u_physical, dEta_bar1, dXi2, dEta3);
266 
267  Array<OneD, const NekDouble> eta_x, eta_y, eta_z;
268  eta_x = m_base[0]->GetZ();
269  eta_y = m_base[1]->GetZ();
270  eta_z = m_base[2]->GetZ();
271 
272  int i, j, k, n;
273 
274  for (k = 0, n = 0; k < Qz; ++k)
275  {
276  for (j = 0; j < Qy; ++j)
277  {
278  for (i = 0; i < Qx; ++i, ++n)
279  {
280  if (out_dxi1.num_elements() > 0)
281  out_dxi1[n] = 2.0/(1.0 - eta_z[k]) * dEta_bar1[n];
282  if (out_dxi2.num_elements() > 0)
283  out_dxi2[n] = 2.0/(1.0 - eta_z[k]) * dXi2[n];
284  if (out_dxi3.num_elements() > 0)
285  out_dxi3[n] = (1.0+eta_x[i])/(1.0-eta_z[k])*dEta_bar1[n] +
286  (1.0+eta_y[j])/(1.0-eta_z[k])*dXi2[n] + dEta3[n];
287  }
288  }
289  }
290  }
291 
292  void StdPyrExp::v_PhysDeriv(const int dir,
293  const Array<OneD, const NekDouble>& inarray,
294  Array<OneD, NekDouble>& outarray)
295  {
296  switch(dir)
297  {
298  case 0:
299  {
300  v_PhysDeriv(inarray, outarray, NullNekDouble1DArray,
301  NullNekDouble1DArray);
302  break;
303  }
304 
305  case 1:
306  {
307  v_PhysDeriv(inarray, NullNekDouble1DArray, outarray,
308  NullNekDouble1DArray);
309  break;
310  }
311 
312  case 2:
313  {
314  v_PhysDeriv(inarray, NullNekDouble1DArray,
315  NullNekDouble1DArray, outarray);
316  break;
317  }
318 
319  default:
320  {
321  ASSERTL1(false,"input dir is out of range");
322  }
323  break;
324  }
325  }
326 
327  void StdPyrExp::v_StdPhysDeriv(const Array<OneD, const NekDouble> &inarray,
328  Array<OneD, NekDouble> &out_d0,
329  Array<OneD, NekDouble> &out_d1,
330  Array<OneD, NekDouble> &out_d2)
331  {
332  StdPyrExp::v_PhysDeriv(inarray, out_d0, out_d1, out_d2);
333  }
334 
335  void StdPyrExp::v_StdPhysDeriv(const int dir,
336  const Array<OneD, const NekDouble>& inarray,
337  Array<OneD, NekDouble>& outarray)
338  {
339  StdPyrExp::v_PhysDeriv(dir, inarray, outarray);
340  }
341 
342  //---------------------------------------
343  // Transforms
344  //---------------------------------------
345 
346  /**
347  * \brief Backward transformation is evaluated at the quadrature
348  * points.
349  *
350  * \f$ u^{\delta} (\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{m(pqr)} \hat
351  * u_{pqr} \phi_{pqr} (\xi_{1i}, \xi_{2j}, \xi_{3k})\f$
352  *
353  * Backward transformation is three dimensional tensorial expansion
354  *
355  * \f$ u (\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{p=0}^{Q_x} \psi_p^a
356  * (\xi_{1i}) \lbrace { \sum_{q=0}^{Q_y} \psi_{q}^a (\xi_{2j})
357  * \lbrace { \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{pqr}^c (\xi_{3k})
358  * \rbrace} \rbrace}. \f$
359  *
360  * And sumfactorizing step of the form is as:\ \ \f$ f_{pqr}
361  * (\xi_{3k}) = \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{pqr}^c
362  * (\xi_{3k}),\\ g_{p} (\xi_{2j}, \xi_{3k}) = \sum_{r=0}^{Q_y}
363  * \psi_{p}^a (\xi_{2j}) f_{pqr} (\xi_{3k}),\\ u(\xi_{1i}, \xi_{2j},
364  * \xi_{3k}) = \sum_{p=0}^{Q_x} \psi_{p}^a (\xi_{1i}) g_{p}
365  * (\xi_{2j}, \xi_{3k}). \f$
366  **/
367  void StdPyrExp::v_BwdTrans(const Array<OneD, const NekDouble> &inarray,
368  Array<OneD, NekDouble> &outarray)
369  {
370  if (m_base[0]->Collocation() &&
371  m_base[1]->Collocation() &&
372  m_base[2]->Collocation())
373  {
374  Vmath::Vcopy(m_base[0]->GetNumPoints()*
375  m_base[1]->GetNumPoints()*
376  m_base[2]->GetNumPoints(),
377  inarray, 1, outarray, 1);
378  }
379  else
380  {
381  StdPyrExp::v_BwdTrans_SumFac(inarray,outarray);
382  }
383  }
384 
385  /**
386  * Sum-factorisation implementation of the BwdTrans operation.
387  */
388  void StdPyrExp::v_BwdTrans_SumFac(
389  const Array<OneD, const NekDouble>& inarray,
390  Array<OneD, NekDouble>& outarray)
391  {
392  Array<OneD, NekDouble> wsp;
393  v_BwdTrans_SumFacKernel(m_base[0]->GetBdata(),
394  m_base[1]->GetBdata(),
395  m_base[2]->GetBdata(),
396  inarray,outarray,wsp,
397  true,true,true);
398  }
399 
400  void StdPyrExp::v_BwdTrans_SumFacKernel(
401  const Array<OneD, const NekDouble> &base0,
402  const Array<OneD, const NekDouble> &base1,
403  const Array<OneD, const NekDouble> &base2,
404  const Array<OneD, const NekDouble> &inarray,
405  Array<OneD, NekDouble> &outarray,
406  Array<OneD, NekDouble> &wsp,
407  bool doCheckCollDir0,
408  bool doCheckCollDir1,
409  bool doCheckCollDir2)
410  {
411  const int Qx = m_base[0]->GetNumPoints();
412  const int Qy = m_base[1]->GetNumPoints();
413  const int Qz = m_base[2]->GetNumPoints();
414 
415  const NekDouble *bx = base0.get();
416  const NekDouble *by = base1.get();
417  const NekDouble *bz = base2.get();
418 
419  // Need to count coeffs for storage...
420  map<int, map<int, map<int, pair<int, int> > > >::iterator it_p;
421  map<int, map<int, pair<int, int> > > ::iterator it_q;
422  map<int, pair<int, int> > ::iterator it_r;
423 
424  int pqcnt = 0;
425  for (it_p = m_idxMap.begin(); it_p != m_idxMap.end(); ++it_p)
426  {
427  for (it_q = it_p->second.begin(); it_q != it_p->second.end(); ++it_q)
428  {
429  pqcnt++;
430  }
431  }
432 
433  Array<OneD, NekDouble> fpq(pqcnt);
434  Array<OneD, NekDouble> fp (m_base[0]->GetNumModes());
435  int i ,j, k, s = 0, cnt = 0, cnt2 = 0;
436 
437  for (k = 0; k < Qz; ++k)
438  {
439  NekDouble bz1 = bz[k+Qz];
440 
441  cnt = 0;
442  for (it_p = m_idxMap.begin(); it_p != m_idxMap.end(); ++it_p)
443  {
444  for (it_q = it_p->second.begin(); it_q != it_p->second.end(); ++it_q)
445  {
446  NekDouble sum = 0.0;
447  for (it_r = it_q->second.begin(); it_r != it_q->second.end(); ++it_r)
448  {
449  sum += inarray[it_r->second.first] * bz[k + Qz*it_r->second.second];
450  }
451  fpq[cnt++] = sum;
452  }
453  }
454 
455  for (j = 0; j < Qy; ++j)
456  {
457  NekDouble by0 = bz1*by[j];
458  NekDouble by1 = bz1*by[j+Qy];
459 
460  cnt = cnt2 = 0;
461  for (it_p = m_idxMap.begin(); it_p != m_idxMap.end(); ++it_p)
462  {
463  NekDouble sum = 0.0;
464  for (it_q = it_p->second.begin(); it_q != it_p->second.end(); ++it_q)
465  {
466  sum += by[j + Qy*it_q->first] * fpq[cnt++];
467  }
468  fp[cnt2++] = sum;
469  }
470 
471  for (i = 0; i < Qx; ++i, ++s)
472  {
473  cnt2 = 0;
474  NekDouble sum = 0.0;
475  for (it_p = m_idxMap.begin(); it_p != m_idxMap.end(); ++it_p)
476  {
477  sum += bx[i + Qx*it_p->first] * fp[cnt2++];
478  }
479  sum += inarray[4]*(by1*(bx[i] + bx[i+Qx]) + by0*bx[i+Qx]);
480  outarray[s] = sum;
481  }
482  }
483  }
484  }
485 
486  /** \brief Forward transform from physical quadrature space
487  stored in \a inarray and evaluate the expansion coefficients and
488  store in \a outarray
489 
490  Inputs:\n
491 
492  - \a inarray: array of physical quadrature points to be transformed
493 
494  Outputs:\n
495 
496  - \a outarray: updated array of expansion coefficients.
497 
498  */
499  void StdPyrExp::v_FwdTrans(const Array<OneD, const NekDouble> &inarray,
500  Array<OneD, NekDouble> &outarray)
501  {
502  v_IProductWRTBase(inarray,outarray);
503 
504  // get Mass matrix inverse
505  StdMatrixKey masskey(eInvMass,DetShapeType(),*this);
506  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
507 
508  // copy inarray in case inarray == outarray
509  DNekVec in (m_ncoeffs, outarray);
510  DNekVec out(m_ncoeffs, outarray, eWrapper);
511 
512  out = (*matsys)*in;
513  }
514 
515 
516  //---------------------------------------
517  // Inner product functions
518  //---------------------------------------
519 
520  /** \brief Inner product of \a inarray over region with respect to the
521  expansion basis m_base[0]->GetBdata(),m_base[1]->GetBdata(), m_base[2]->GetBdata() and return in \a outarray
522 
523  Wrapper call to StdPyrExp::IProductWRTBase
524 
525  Input:\n
526 
527  - \a inarray: array of function evaluated at the physical collocation points
528 
529  Output:\n
530 
531  - \a outarray: array of inner product with respect to each basis over region
532 
533  */
534  void StdPyrExp::v_IProductWRTBase(
535  const Array<OneD, const NekDouble> &inarray,
536  Array<OneD, NekDouble> &outarray)
537  {
538  if (m_base[0]->Collocation() &&
539  m_base[1]->Collocation() &&
540  m_base[2]->Collocation())
541  {
542  MultiplyByStdQuadratureMetric(inarray, outarray);
543  }
544  else
545  {
546  StdPyrExp::v_IProductWRTBase_SumFac(inarray,outarray);
547  }
548  }
549 
550  void StdPyrExp::v_IProductWRTBase_SumFac(
551  const Array<OneD, const NekDouble>& inarray,
552  Array<OneD, NekDouble>& outarray,
553  bool multiplybyweights)
554  {
555  Array<OneD, NekDouble> wsp;
556 
557  if(multiplybyweights)
558  {
559  Array<OneD, NekDouble> tmp(inarray.num_elements());
560 
561  v_MultiplyByStdQuadratureMetric(inarray, tmp);
562 
563  v_IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
564  m_base[1]->GetBdata(),
565  m_base[2]->GetBdata(),
566  tmp,outarray,wsp,
567  true,true,true);
568  }
569  else
570  {
571  v_IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
572  m_base[1]->GetBdata(),
573  m_base[2]->GetBdata(),
574  inarray,outarray,wsp,
575  true,true,true);
576  }
577  }
578 
579  void StdPyrExp::v_IProductWRTBase_SumFacKernel(
580  const Array<OneD, const NekDouble> &base0,
581  const Array<OneD, const NekDouble> &base1,
582  const Array<OneD, const NekDouble> &base2,
583  const Array<OneD, const NekDouble> &inarray,
584  Array<OneD, NekDouble> &outarray,
585  Array<OneD, NekDouble> &wsp,
586  bool doCheckCollDir0,
587  bool doCheckCollDir1,
588  bool doCheckCollDir2)
589  {
590  int i, j, k, s;
591  int Qx = m_base[0]->GetNumPoints();
592  int Qy = m_base[1]->GetNumPoints();
593  int Qz = m_base[2]->GetNumPoints();
594 
595  const NekDouble *bx = base0.get();
596  const NekDouble *by = base1.get();
597  const NekDouble *bz = base2.get();
598 
599  map<int, map<int, map<int, pair<int, int> > > >::iterator it_p;
600  map<int, map<int, pair<int, int> > > ::iterator it_q;
601  map<int, pair<int, int> > ::iterator it_r;
602 
603  Array<OneD, NekDouble> f (Qy*Qz);
604  Array<OneD, NekDouble> fb(Qz);
605 
606  for (it_p = m_idxMap.begin(); it_p != m_idxMap.end(); ++it_p)
607  {
608  const int p = it_p->first;
609  s = 0;
610  for (k = 0; k < Qz; ++k)
611  {
612  for (j = 0; j < Qy; ++j)
613  {
614  NekDouble sum = 0.0;
615  for (i = 0; i < Qx; ++i, ++s)
616  {
617  sum += bx[i + Qx*p]*inarray[s];
618  }
619  f[j+Qy*k] = sum;
620  }
621  }
622 
623  for (it_q = it_p->second.begin(); it_q != it_p->second.end(); ++it_q)
624  {
625  const int q = it_q->first;
626 
627  for (k = 0; k < Qz; ++k)
628  {
629  NekDouble sum = 0.0;
630  for (j = 0; j < Qy; ++j)
631  {
632  sum += by[j + Qy*q]*f[j+Qy*k];
633  }
634  fb[k] = sum;
635  }
636 
637  for (it_r = it_q->second.begin(); it_r != it_q->second.end(); ++it_r)
638  {
639  const int rpqr = it_r->second.second;
640  NekDouble sum = 0.0;
641  for (k = 0; k < Qz; ++k)
642  {
643  sum += bz[k + Qz*rpqr]*fb[k];
644  }
645 
646  outarray[it_r->second.first] = sum;
647  }
648  }
649  }
650 
651  // Correct for top mode
652  s = 0;
653  for (k = 0; k < Qz; ++k)
654  {
655  for (j = 0; j < Qy; ++j)
656  {
657  for (i = 0; i < Qx; ++i, ++s)
658  {
659  outarray[4] += inarray[s] * bz[k+Qz]*(
660  bx[i+Qx]*by[j+Qy] +
661  bx[i+Qx]*by[j ] +
662  bx[i ]*by[j+Qy]);
663  }
664  }
665  }
666  }
667 
668  void StdPyrExp::v_IProductWRTDerivBase(
669  const int dir,
670  const Array<OneD, const NekDouble>& inarray,
671  Array<OneD, NekDouble>& outarray)
672  {
673  StdPyrExp::v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
674  }
675 
676  /**
677  * @param inarray Function evaluated at physical collocation
678  * points.
679  * @param outarray Inner product with respect to each basis
680  * function over the element.
681  */
682  void StdPyrExp::v_IProductWRTDerivBase_SumFac(
683  const int dir,
684  const Array<OneD, const NekDouble>& inarray,
685  Array<OneD, NekDouble>& outarray)
686  {
687  int i;
688  int nquad0 = m_base[0]->GetNumPoints();
689  int nquad1 = m_base[1]->GetNumPoints();
690  int nquad2 = m_base[2]->GetNumPoints();
691  int nqtot = nquad0*nquad1*nquad2;
692 
693  Array<OneD, NekDouble> gfac0(nquad0);
694  Array<OneD, NekDouble> gfac1(nquad1);
695  Array<OneD, NekDouble> gfac2(nquad2);
696  Array<OneD, NekDouble> tmp0 (nqtot);
697  Array<OneD, NekDouble> wsp;
698 
699  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
700  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
701  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
702 
703  // set up geometric factor: (1+z0)/2
704  for(i = 0; i < nquad0; ++i)
705  {
706  gfac0[i] = 0.5*(1+z0[i]);
707  }
708 
709  // set up geometric factor: (1+z1)/2
710  for(i = 0; i < nquad1; ++i)
711  {
712  gfac1[i] = 0.5*(1+z1[i]);
713  }
714 
715  // Set up geometric factor: 2/(1-z2)
716  for(i = 0; i < nquad2; ++i)
717  {
718  gfac2[i] = 2.0/(1-z2[i]);
719  }
720 
721  // Derivative in first/second direction is always scaled as follows
722  const int nq01 = nquad0*nquad1;
723  for(i = 0; i < nquad2; ++i)
724  {
725  Vmath::Smul(nq01, gfac2[i],
726  &inarray[0] + i*nq01, 1,
727  &tmp0 [0] + i*nq01, 1);
728  }
729 
730  MultiplyByStdQuadratureMetric(tmp0, tmp0);
731 
732  switch(dir)
733  {
734  case 0:
735  {
736  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
737  m_base[1]->GetBdata (),
738  m_base[2]->GetBdata (),
739  tmp0, outarray, wsp,
740  false, true, true);
741  break;
742  }
743  case 1:
744  {
745  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata (),
746  m_base[1]->GetDbdata(),
747  m_base[2]->GetBdata (),
748  tmp0, outarray, wsp,
749  true, false, true);
750  break;
751  }
752  case 2:
753  {
754  Array<OneD, NekDouble> tmp3(m_ncoeffs);
755  Array<OneD, NekDouble> tmp4(m_ncoeffs);
756 
757  // Scale eta_1 derivative by gfac0
758  for(i = 0; i < nquad1*nquad2; ++i)
759  {
760  Vmath::Vmul(nquad0, tmp0 .get() + i*nquad0, 1,
761  gfac0.get(), 1,
762  tmp0 .get() + i*nquad0, 1);
763  }
764  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
765  m_base[1]->GetBdata(),
766  m_base[2]->GetBdata(),
767  tmp0, tmp3, wsp,
768  false, true, true);
769 
770  // Scale eta_2 derivative by gfac1*gfac2
771  for(i = 0; i < nquad2; ++i)
772  {
773  Vmath::Smul(nq01, gfac2[i],
774  &inarray[0] + i*nq01, 1,
775  &tmp0 [0] + i*nq01, 1);
776  }
777  for(i = 0; i < nquad1*nquad2; ++i)
778  {
779  Vmath::Smul(nquad0, gfac1[i%nquad1],
780  &tmp0[0] + i*nquad0, 1,
781  &tmp0[0] + i*nquad0, 1);
782  }
783 
784  MultiplyByStdQuadratureMetric(tmp0, tmp0);
785  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
786  m_base[1]->GetDbdata(),
787  m_base[2]->GetBdata(),
788  tmp0, tmp4, wsp,
789  true, false, true);
790 
791  MultiplyByStdQuadratureMetric(inarray,tmp0);
792  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
793  m_base[1]->GetBdata(),
794  m_base[2]->GetDbdata(),
795  tmp0,outarray,wsp,
796  true, true, false);
797 
798  Vmath::Vadd(m_ncoeffs,&tmp3[0],1,&outarray[0],1,&outarray[0],1);
799  Vmath::Vadd(m_ncoeffs,&tmp4[0],1,&outarray[0],1,&outarray[0],1);
800  break;
801  }
802  default:
803  {
804  ASSERTL1(false, "input dir is out of range");
805  break;
806  }
807  }
808  }
809 
810  //---------------------------------------
811  // Evaluation functions
812  //---------------------------------------
813 
814  void StdPyrExp::v_LocCoordToLocCollapsed(
815  const Array<OneD, const NekDouble>& xi,
816  Array<OneD, NekDouble>& eta)
817  {
818  if (fabs(xi[2]-1.0) < NekConstants::kNekZeroTol)
819  {
820  // Very top point of the pyramid
821  eta[0] = -1.0;
822  eta[1] = -1.0;
823  eta[2] = xi[2];
824  }
825  else
826  {
827  // Below the line-singularity -- Common case
828  eta[2] = xi[2]; // eta_z = xi_z
829  eta[1] = 2.0*(1.0 + xi[1])/(1.0 - xi[2]) - 1.0;
830  eta[0] = 2.0*(1.0 + xi[0])/(1.0 - xi[2]) - 1.0;
831  }
832  }
833 
834  void StdPyrExp::v_GetCoords(Array<OneD, NekDouble> &xi_x,
835  Array<OneD, NekDouble> &xi_y,
836  Array<OneD, NekDouble> &xi_z)
837  {
838  Array<OneD, const NekDouble> etaBar_x = m_base[0]->GetZ();
839  Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
840  Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
841  int Qx = GetNumPoints(0);
842  int Qy = GetNumPoints(1);
843  int Qz = GetNumPoints(2);
844 
845  // Convert collapsed coordinates into cartesian coordinates: eta --> xi
846  for (int k = 0; k < Qz; ++k )
847  {
848  for (int j = 0; j < Qy; ++j)
849  {
850  for (int i = 0; i < Qx; ++i)
851  {
852  int s = i + Qx*(j + Qy*k);
853 
854  xi_z[s] = eta_z[k];
855  xi_y[s] = (1.0 + eta_y[j]) * (1.0 - eta_z[k]) / 2.0 - 1.0;
856  xi_x[s] = (1.0 + etaBar_x[i]) * (1.0 - eta_z[k]) / 2.0 - 1.0;
857  }
858  }
859  }
860  }
861 
862  void StdPyrExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
863  {
864  Array<OneD, NekDouble> tmp(m_ncoeffs, 0.0);
865  tmp[mode] = 1.0;
866  v_BwdTrans(tmp, outarray);
867  }
868 
869 
870  //---------------------------------------
871  // Helper functions
872  //---------------------------------------
873 
874  int StdPyrExp::v_GetNverts() const
875  {
876  return 5;
877  }
878 
879  int StdPyrExp::v_GetNedges() const
880  {
881  return 8;
882  }
883 
884  int StdPyrExp::v_GetNfaces() const
885  {
886  return 5;
887  }
888 
889  LibUtilities::ShapeType StdPyrExp::v_DetShapeType() const
890  {
891  return LibUtilities::ePyramid;
892  }
893 
894  int StdPyrExp::v_NumBndryCoeffs() const
895  {
896  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
897  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
898  "BasisType is not a boundary interior form");
899  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
900  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
901  "BasisType is not a boundary interior form");
902  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
903  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
904  "BasisType is not a boundary interior form");
905 
906  int P = m_base[0]->GetNumModes();
907  int Q = m_base[1]->GetNumModes();
908  int R = m_base[2]->GetNumModes();
909 
910  return LibUtilities::StdPyrData::
911  getNumberOfBndCoefficients(P, Q, R);
912  }
913 
914  int StdPyrExp::v_GetEdgeNcoeffs(const int i) const
915  {
916  ASSERTL2(i >= 0 && i <= 7, "edge id is out of range");
917 
918  if (i == 0 || i == 2)
919  {
920  return GetBasisNumModes(0);
921  }
922  else if (i == 1 || i == 3)
923  {
924  return GetBasisNumModes(1);
925  }
926  else
927  {
928  return GetBasisNumModes(2);
929  }
930  }
931 
932  int StdPyrExp::v_GetFaceNcoeffs(const int i) const
933  {
934  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
935 
936  if (i == 0)
937  {
938  return GetBasisNumModes(0)*GetBasisNumModes(1);
939  }
940  else if (i == 1 || i == 3)
941  {
942  int P = GetBasisNumModes(0)-1, Q = GetBasisNumModes(2)-1;
943  return Q+1 + (P*(1 + 2*Q - P))/2;
944  }
945  else
946  {
947  int P = GetBasisNumModes(1)-1, Q = GetBasisNumModes(2)-1;
948  return Q+1 + (P*(1 + 2*Q - P))/2;
949  }
950  }
951 
952  int StdPyrExp::v_GetFaceIntNcoeffs(const int i) const
953  {
954  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
955 
956  int P = m_base[0]->GetNumModes()-1;
957  int Q = m_base[1]->GetNumModes()-1;
958  int R = m_base[2]->GetNumModes()-1;
959 
960  if (i == 0)
961  {
962  return (P-1)*(Q-1);
963  }
964  else if (i == 1 || i == 3)
965  {
966  return (P-1) * (2*(R-1) - (P-1) - 1) / 2;
967  }
968  else
969  {
970  return (Q-1) * (2*(R-1) - (Q-1) - 1) / 2;
971  }
972  }
973 
974  int StdPyrExp::v_GetFaceNumPoints(const int i) const
975  {
976  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
977 
978  if (i == 0)
979  {
980  return m_base[0]->GetNumPoints()*
981  m_base[1]->GetNumPoints();
982  }
983  else if (i == 1 || i == 3)
984  {
985  return m_base[0]->GetNumPoints()*
986  m_base[2]->GetNumPoints();
987  }
988  else
989  {
990  return m_base[1]->GetNumPoints()*
991  m_base[2]->GetNumPoints();
992  }
993  }
994 
995 
996  const LibUtilities::BasisKey StdPyrExp::v_DetFaceBasisKey(
997  const int i, const int k) const
998  {
999  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
1000  ASSERTL2(k >= 0 && k <= 1, "basis key id is out of range");
1001 
1002  switch(i)
1003  {
1004  case 0:
1005  {
1006  return EvaluateQuadFaceBasisKey(k,
1007  m_base[k]->GetBasisType(),
1008  m_base[k]->GetNumPoints(),
1009  m_base[k]->GetNumModes());
1010 
1011  }
1012  case 1:
1013  case 3:
1014  {
1015  return EvaluateTriFaceBasisKey(k,
1016  m_base[2*k]->GetBasisType(),
1017  m_base[2*k]->GetNumPoints(),
1018  m_base[2*k]->GetNumModes());
1019  }
1020  case 2:
1021  case 4:
1022  {
1023  return EvaluateTriFaceBasisKey(k,
1024  m_base[k+1]->GetBasisType(),
1025  m_base[k+1]->GetNumPoints(),
1026  m_base[k+1]->GetNumModes());
1027  }
1028  }
1029 
1030  // Should never get here.
1031  return LibUtilities::NullBasisKey;
1032  }
1033 
1034  int StdPyrExp::v_CalcNumberOfCoefficients(
1035  const std::vector<unsigned int> &nummodes,
1036  int &modes_offset)
1037  {
1038  int nmodes = LibUtilities::StdPyrData::getNumberOfCoefficients(
1039  nummodes[modes_offset],
1040  nummodes[modes_offset+1],
1041  nummodes[modes_offset+2]);
1042 
1043  modes_offset += 3;
1044  return nmodes;
1045  }
1046 
1047  LibUtilities::BasisType StdPyrExp::v_GetEdgeBasisType(const int i) const
1048  {
1049  ASSERTL2(i >= 0 && i <= 7, "edge id is out of range");
1050  if (i == 0 || i == 2)
1051  {
1052  return GetBasisType(0);
1053  }
1054  else if (i == 1 || i == 3)
1055  {
1056  return GetBasisType(1);
1057  }
1058  else
1059  {
1060  return GetBasisType(2);
1061  }
1062  }
1063 
1064 
1065  //---------------------------------------
1066  // Mappings
1067  //---------------------------------------
1068 
1069  void StdPyrExp::v_GetFaceToElementMap(
1070  const int fid,
1071  const Orientation faceOrient,
1072  Array<OneD, unsigned int> &maparray,
1073  Array<OneD, int> &signarray,
1074  int P,
1075  int Q)
1076  {
1077  ASSERTL1(GetEdgeBasisType(0) == GetEdgeBasisType(1),
1078  "Method only implemented if BasisType is identical"
1079  "in x and y directions");
1080  ASSERTL1(GetEdgeBasisType(0) == LibUtilities::eModified_A &&
1081  GetEdgeBasisType(4) == LibUtilities::eModified_C,
1082  "Method only implemented for Modified_A BasisType"
1083  "(x and y direction) and Modified_C BasisType (z "
1084  "direction)");
1085 
1086  int i, j, k, p, q, r, nFaceCoeffs;
1087  int nummodesA, nummodesB;
1088 
1089  int order0 = m_base[0]->GetNumModes();
1090  int order1 = m_base[1]->GetNumModes();
1091  int order2 = m_base[2]->GetNumModes();
1092 
1093  switch (fid)
1094  {
1095  case 0:
1096  nummodesA = order0;
1097  nummodesB = order1;
1098  break;
1099  case 1:
1100  case 3:
1101  nummodesA = order0;
1102  nummodesB = order2;
1103  break;
1104  case 2:
1105  case 4:
1106  nummodesA = order1;
1107  nummodesB = order2;
1108  break;
1109  }
1110 
1111  bool CheckForZeroedModes = false;
1112 
1113  if (P == -1)
1114  {
1115  P = nummodesA;
1116  Q = nummodesB;
1117  nFaceCoeffs = GetFaceNcoeffs(fid);
1118  }
1119  else if (fid > 0)
1120  {
1121  nFaceCoeffs = P*(2*Q-P+1)/2;
1122  CheckForZeroedModes = true;
1123  }
1124  else
1125  {
1126  nFaceCoeffs = P*Q;
1127  CheckForZeroedModes = true;
1128  }
1129 
1130  // Allocate the map array and sign array; set sign array to ones (+)
1131  if (maparray.num_elements() != nFaceCoeffs)
1132  {
1133  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1134  }
1135 
1136  if (signarray.num_elements() != nFaceCoeffs)
1137  {
1138  signarray = Array<OneD, int>(nFaceCoeffs,1);
1139  }
1140  else
1141  {
1142  fill(signarray.get(), signarray.get() + nFaceCoeffs, 1);
1143  }
1144 
1145  // Set up an array indexing for quads, since the ordering may need
1146  // to be transposed.
1147  Array<OneD, int> arrayindx(nFaceCoeffs,-1);
1148 
1149  if (fid == 0)
1150  {
1151  for (i = 0; i < Q; i++)
1152  {
1153  for (j = 0; j < P; j++)
1154  {
1155  if (faceOrient < 9)
1156  {
1157  arrayindx[i*P+j] = i*P+j;
1158  }
1159  else
1160  {
1161  arrayindx[i*P+j] = j*Q+i;
1162  }
1163  }
1164  }
1165  }
1166 
1167  // Set up ordering inside each 2D face. Also for triangular faces,
1168  // populate signarray.
1169  int cnt = 0, cnt2;
1170  switch (fid)
1171  {
1172  case 0: // Bottom quad
1173 
1174  // Fill in vertices
1175  maparray[arrayindx[0]] = 0;
1176  maparray[arrayindx[1]] = 1;
1177  maparray[arrayindx[P+1]] = 2;
1178  maparray[arrayindx[P]] = 3;
1179 
1180  // Edge 0
1181  cnt = 5;
1182  for (p = 2; p < P; ++p)
1183  {
1184  maparray[arrayindx[p]] = p-2 + cnt;
1185  }
1186 
1187  // Edge 1
1188  cnt += P-2;
1189  for (q = 2; q < Q; ++q)
1190  {
1191  maparray[arrayindx[q*P+1]] = q-2 + cnt;
1192  }
1193 
1194  // Edge 2
1195  cnt += Q-2;
1196  for (p = 2; p < P; ++p)
1197  {
1198  maparray[arrayindx[P+p]] = p-2 + cnt;
1199  }
1200 
1201  // Edge 3
1202  cnt += P-2;
1203  for (q = 2; q < Q; ++q)
1204  {
1205  maparray[arrayindx[q*P]] = q-2 + cnt;
1206  }
1207 
1208  // Interior
1209  cnt += Q-2 + 4*(P-2);
1210  for (q = 2; q < Q; ++q)
1211  {
1212  for (p = 2; p < P; ++p)
1213  {
1214  maparray[arrayindx[q*P+p]] = cnt + (q-2)*P+(p-2);
1215  }
1216  }
1217 
1218 
1219  break;
1220 
1221  case 1: // Left triangle
1222  // Vertices
1223  maparray[0] = 0;
1224  maparray[1] = 4;
1225  maparray[Q] = 1;
1226 
1227  // Edge 0 (pyramid edge 0)
1228  cnt = 5;
1229  q = 2*Q-1;
1230  for (p = 2; p < P; q += Q-p, ++p)
1231  {
1232  maparray[q] = cnt++;
1233  if ((int)faceOrient == 7)
1234  {
1235  signarray[q] = p % 2 ? -1 : 1;
1236  }
1237  }
1238 
1239  // Edge 1 (pyramid edge 5)
1240  cnt = 5 + 2*(order0-2) + 2*(order1-2) + (order2-2);
1241  for (q = 2; q < Q; ++q)
1242  {
1243  maparray[q] = cnt++;
1244  if ((int)faceOrient == 7)
1245  {
1246  signarray[q] = q % 2 ? -1 : 1;
1247  }
1248  }
1249 
1250  // Edge 2 (pyramid edge 4)
1251  cnt = 5 + 2*(order0-2) + 2*(order1-2);
1252  for (q = 2; q < Q; ++q)
1253  {
1254  maparray[Q+q-1] = cnt++;
1255  if ((int)faceOrient == 7)
1256  {
1257  signarray[Q+q-1] = q % 2 ? -1 : 1;
1258  }
1259  }
1260 
1261  // Interior
1262  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1263  + v_GetFaceIntNcoeffs(0);
1264  cnt2 = 2*Q + 1;
1265  for (p = 2; p < P; ++p)
1266  {
1267  for (r = 2; r < Q-p; ++r)
1268  {
1269  maparray[cnt2] = cnt++;
1270  if ((int)faceOrient == 7 && p > 1)
1271  {
1272  signarray[cnt2++] = p % 2 ? -1 : 1;
1273  }
1274  }
1275  cnt2++;
1276  }
1277  break;
1278 
1279  case 2:
1280  // Vertices
1281  maparray[0] = 1;
1282  maparray[1] = 4;
1283  maparray[Q] = 2;
1284 
1285  // Edge 0 (pyramid edge 1)
1286  cnt = 5 + (order0-2);
1287  q = 2*Q-1;
1288  for (p = 2; p < P; q += Q-p, ++p)
1289  {
1290  maparray[q] = cnt++;
1291  if ((int)faceOrient == 7)
1292  {
1293  signarray[q] = p % 2 ? -1 : 1;
1294  }
1295  }
1296 
1297  // Edge 1 (pyramid edge 6)
1298  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 2*(order2-2);
1299  for (q = 2; q < Q; ++q)
1300  {
1301  maparray[q] = cnt++;
1302  if ((int)faceOrient == 7)
1303  {
1304  signarray[q] = q % 2 ? -1 : 1;
1305  }
1306  }
1307 
1308  // Edge 2 (pyramid edge 5)
1309  cnt = 5 + 2*(order0-2) + 2*(order1-2) + (order2-2);
1310  for (q = 2; q < Q; ++q)
1311  {
1312  maparray[Q+q-1] = cnt++;
1313  if ((int)faceOrient == 7)
1314  {
1315  signarray[Q+q-1] = q % 2 ? -1 : 1;
1316  }
1317  }
1318 
1319  // Interior
1320  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1321  + v_GetFaceIntNcoeffs(0) + v_GetFaceIntNcoeffs(1);
1322  cnt2 = 2*Q + 1;
1323  for (p = 2; p < P; ++p)
1324  {
1325  for (r = 2; r < Q-p; ++r)
1326  {
1327  maparray[cnt2] = cnt++;
1328  if ((int)faceOrient == 7 && p > 1)
1329  {
1330  signarray[cnt2++] = p % 2 ? -1 : 1;
1331  }
1332  }
1333  cnt2++;
1334  }
1335  break;
1336 
1337  case 3: // Right triangle
1338  // Vertices
1339  maparray[0] = 3;
1340  maparray[1] = 4;
1341  maparray[Q] = 2;
1342 
1343  // Edge 0 (pyramid edge 2)
1344  cnt = 5 + (order0-2) + (order1-2);
1345  q = 2*Q-1;
1346  for (p = 2; p < P; q += Q-p, ++p)
1347  {
1348  maparray[q] = cnt++;
1349  if ((int)faceOrient == 7)
1350  {
1351  signarray[q] = p % 2 ? -1 : 1;
1352  }
1353  }
1354 
1355  // Edge 1 (pyramid edge 6)
1356  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 2*(order2-2);
1357  for (q = 2; q < Q; ++q)
1358  {
1359  maparray[q] = cnt++;
1360  if ((int)faceOrient == 7)
1361  {
1362  signarray[q] = q % 2 ? -1 : 1;
1363  }
1364  }
1365 
1366  // Edge 2 (pyramid edge 7)
1367  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 3*(order2-2);
1368  for (q = 2; q < Q; ++q)
1369  {
1370  maparray[Q+q-1] = cnt++;
1371  if ((int)faceOrient == 7)
1372  {
1373  signarray[Q+q-1] = q % 2 ? -1 : 1;
1374  }
1375  }
1376 
1377  // Interior
1378  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1379  + v_GetFaceIntNcoeffs(0) + v_GetFaceIntNcoeffs(1)
1380  + v_GetFaceIntNcoeffs(2);
1381  cnt2 = 2*Q + 1;
1382  for (p = 2; p < P; ++p)
1383  {
1384  for (r = 2; r < Q-p; ++r)
1385  {
1386  maparray[cnt2] = cnt++;
1387  if ((int)faceOrient == 7 && p > 1)
1388  {
1389  signarray[cnt2++] = p % 2 ? -1 : 1;
1390  }
1391  }
1392  cnt2++;
1393  }
1394  break;
1395 
1396  case 4: // Rear tri
1397  // Vertices
1398  maparray[0] = 0;
1399  maparray[1] = 4;
1400  maparray[Q] = 3;
1401 
1402  // Edge 0 (pyramid edge 3)
1403  cnt = 5 + 2*(order0-2) + (order1-2);
1404  q = 2*Q-1;
1405  for (p = 2; p < P; q += Q-p, ++p)
1406  {
1407  maparray[q] = cnt++;
1408  if ((int)faceOrient == 7)
1409  {
1410  signarray[q] = p % 2 ? -1 : 1;
1411  }
1412  }
1413 
1414  // Edge 1 (pyramid edge 7)
1415  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 3*(order2-2);
1416  for (q = 2; q < Q; ++q)
1417  {
1418  maparray[q] = cnt++;
1419  if ((int)faceOrient == 7)
1420  {
1421  signarray[q] = q % 2 ? -1 : 1;
1422  }
1423  }
1424 
1425  // Edge 2 (pyramid edge 4)
1426  cnt = 5 + 2*(order0-2) + 2*(order1-2);
1427  for (q = 2; q < Q; ++q)
1428  {
1429  maparray[Q+q-1] = cnt++;
1430  if ((int)faceOrient == 7)
1431  {
1432  signarray[Q+q-1] = q % 2 ? -1 : 1;
1433  }
1434  }
1435 
1436  // Interior
1437  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1438  + v_GetFaceIntNcoeffs(0) + v_GetFaceIntNcoeffs(1)
1439  + v_GetFaceIntNcoeffs(2) + v_GetFaceIntNcoeffs(3);
1440  cnt2 = 2*Q + 1;
1441  for (p = 2; p < P; ++p)
1442  {
1443  for (r = 2; r < Q-p; ++r)
1444  {
1445  maparray[cnt2] = cnt++;
1446  if ((int)faceOrient == 7 && p > 1)
1447  {
1448  signarray[cnt2++] = p % 2 ? -1 : 1;
1449  }
1450  }
1451  cnt2++;
1452  }
1453  break;
1454 
1455  default:
1456  ASSERTL0(false, "Face to element map unavailable.");
1457  }
1458 
1459  if (fid > 0)
1460  {
1461 
1462  if(CheckForZeroedModes)
1463  {
1464  // zero signmap and set maparray to zero if elemental
1465  // modes are not as large as face modesl
1466  int idx = 0;
1467  for (j = 0; j < P; ++j)
1468  {
1469  idx += Q-j;
1470  for (k = Q-j; k < Q-j; ++k)
1471  {
1472  signarray[idx] = 0.0;
1473  maparray[idx++] = maparray[0];
1474  }
1475  }
1476 
1477  for (j = P; j < P; ++j)
1478  {
1479  for (k = 0; k < Q-j; ++k)
1480  {
1481  signarray[idx] = 0.0;
1482  maparray[idx++] = maparray[0];
1483  }
1484  }
1485  }
1486 
1487  // Triangles only have one possible orientation (base
1488  // direction reversed); swap edge modes.
1489  if ((int)faceOrient == 7)
1490  {
1491  swap(maparray[0], maparray[P]);
1492  for (i = 1; i < P-1; ++i)
1493  {
1494  swap(maparray[i+1], maparray[P+i]);
1495  }
1496  }
1497  }
1498  else
1499  {
1500  if(CheckForZeroedModes)
1501  {
1502  // zero signmap and set maparray to zero if elemental
1503  // modes are not as large as face modesl
1504  for (j = 0; j < P; ++j)
1505  {
1506  for (k = Q; k < Q; ++k)
1507  {
1508  signarray[arrayindx[j+k*P]] = 0.0;
1509  maparray[arrayindx[j+k*P]] = maparray[0];
1510  }
1511  }
1512 
1513  for (j = P; j < P; ++j)
1514  {
1515  for (k = 0; k < Q; ++k)
1516  {
1517  signarray[arrayindx[j+k*P]] = 0.0;
1518  maparray[arrayindx[j+k*P]] = maparray[0];
1519  }
1520  }
1521  }
1522 
1523  // The code below is exactly the same as that taken from
1524  // StdHexExp and reverses the 'b' and 'a' directions as
1525  // appropriate (1st and 2nd if statements respectively) in
1526  // quadrilateral faces.
1527  if (faceOrient == 6 || faceOrient == 8 ||
1528  faceOrient == 11 || faceOrient == 12)
1529  {
1530  if (faceOrient < 9)
1531  {
1532  for (i = 3; i < Q; i += 2)
1533  {
1534  for (j = 0; j < P; j++)
1535  {
1536  signarray[arrayindx[i*P+j]] *= -1;
1537  }
1538  }
1539 
1540  for (i = 0; i < P; i++)
1541  {
1542  swap(maparray [i], maparray [i+P]);
1543  swap(signarray[i], signarray[i+P]);
1544  }
1545  }
1546  else
1547  {
1548  for (i = 0; i < Q; i++)
1549  {
1550  for (j = 3; j < P; j += 2)
1551  {
1552  signarray[arrayindx[i*P+j]] *= -1;
1553  }
1554  }
1555 
1556  for (i = 0; i < Q; i++)
1557  {
1558  swap (maparray [i], maparray [i+Q]);
1559  swap (signarray[i], signarray[i+Q]);
1560  }
1561  }
1562  }
1563 
1564  if (faceOrient == 7 || faceOrient == 8 ||
1565  faceOrient == 10 || faceOrient == 12)
1566  {
1567  if (faceOrient < 9)
1568  {
1569  for (i = 0; i < Q; i++)
1570  {
1571  for (j = 3; j < P; j += 2)
1572  {
1573  signarray[arrayindx[i*P+j]] *= -1;
1574  }
1575  }
1576 
1577  for(i = 0; i < Q; i++)
1578  {
1579  swap(maparray [i*P], maparray [i*P+1]);
1580  swap(signarray[i*P], signarray[i*P+1]);
1581  }
1582  }
1583  else
1584  {
1585  for (i = 3; i < Q; i += 2)
1586  {
1587  for (j = 0; j < P; j++)
1588  {
1589  signarray[arrayindx[i*P+j]] *= -1;
1590  }
1591  }
1592 
1593  for (i = 0; i < P; i++)
1594  {
1595  swap(maparray [i*Q], maparray [i*Q+1]);
1596  swap(signarray[i*Q], signarray[i*Q+1]);
1597  }
1598  }
1599  }
1600  }
1601  }
1602 
1603  int StdPyrExp::v_GetVertexMap(int vId, bool useCoeffPacking)
1604  {
1605  ASSERTL0(GetEdgeBasisType(vId) == LibUtilities::eModified_A ||
1606  GetEdgeBasisType(vId) == LibUtilities::eModified_A ||
1607  GetEdgeBasisType(vId) == LibUtilities::eModified_C,
1608  "Mapping not defined for this type of basis");
1609  return vId;
1610  }
1611 
1612  void StdPyrExp::v_GetEdgeInteriorMap(
1613  const int eid,
1614  const Orientation edgeOrient,
1615  Array<OneD, unsigned int> &maparray,
1616  Array<OneD, int> &signarray)
1617  {
1618  int i;
1619  bool signChange;
1620  const int P = m_base[0]->GetNumModes() - 2;
1621  const int Q = m_base[1]->GetNumModes() - 2;
1622  const int R = m_base[2]->GetNumModes() - 2;
1623  const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid) - 2;
1624 
1625  if (maparray.num_elements() != nEdgeIntCoeffs)
1626  {
1627  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1628  }
1629 
1630  if(signarray.num_elements() != nEdgeIntCoeffs)
1631  {
1632  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1633  }
1634  else
1635  {
1636  fill(signarray.get(), signarray.get()+nEdgeIntCoeffs, 1);
1637  }
1638 
1639  // If edge is oriented backwards, change sign of modes which have
1640  // degree 2n+1, n >= 1.
1641  signChange = edgeOrient == eBackwards;
1642 
1643  int offset = 5;
1644 
1645  switch (eid)
1646  {
1647  case 0:
1648  break;
1649  case 1:
1650  offset += P;
1651  break;
1652  case 2:
1653  offset += P+Q;
1654  break;
1655  case 3:
1656  offset += 2*P+Q;
1657  break;
1658  case 4:
1659  offset += 2*(P+Q);
1660  break;
1661  case 5:
1662  offset += 2*(P+Q)+R;
1663  break;
1664  case 6:
1665  offset += 2*(P+Q+R);
1666  break;
1667  case 7:
1668  offset += 2*(P+Q)+3*R;
1669  break;
1670  default:
1671  ASSERTL0(false, "Edge not defined.");
1672  break;
1673  }
1674 
1675  for (i = 0; i < nEdgeIntCoeffs; ++i)
1676  {
1677  maparray[i] = offset + i;
1678  }
1679 
1680  if (signChange)
1681  {
1682  for (i = 1; i < nEdgeIntCoeffs; i += 2)
1683  {
1684  signarray[i] = -1;
1685  }
1686  }
1687  }
1688 
1689  void StdPyrExp::v_GetFaceInteriorMap(
1690  const int fid,
1691  const Orientation faceOrient,
1692  Array<OneD, unsigned int> &maparray,
1693  Array<OneD, int> &signarray)
1694  {
1695  const int P = m_base[0]->GetNumModes() - 1;
1696  const int Q = m_base[1]->GetNumModes() - 1;
1697  const int R = m_base[2]->GetNumModes() - 1;
1698  const int nFaceIntCoeffs = v_GetFaceIntNcoeffs(fid);
1699  int p, q, r, idx = 0;
1700  int nummodesA = 0;
1701  int nummodesB = 0;
1702  int i, j;
1703 
1704  if (maparray.num_elements() != nFaceIntCoeffs)
1705  {
1706  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1707  }
1708 
1709  if (signarray.num_elements() != nFaceIntCoeffs)
1710  {
1711  signarray = Array<OneD, int>(nFaceIntCoeffs, 1);
1712  }
1713  else
1714  {
1715  fill(signarray.get(), signarray.get()+nFaceIntCoeffs, 1);
1716  }
1717 
1718  // Set up an array indexing for quad faces, since the ordering may
1719  // need to be transposed depending on orientation.
1720  Array<OneD, int> arrayindx(nFaceIntCoeffs);
1721  if (fid == 0)
1722  {
1723  nummodesA = P-1;
1724  nummodesB = Q-1;
1725 
1726  for (i = 0; i < nummodesB; i++)
1727  {
1728  for (j = 0; j < nummodesA; j++)
1729  {
1730  if (faceOrient < 9)
1731  {
1732  arrayindx[i*nummodesA+j] = i*nummodesA+j;
1733  }
1734  else
1735  {
1736  arrayindx[i*nummodesA+j] = j*nummodesB+i;
1737  }
1738  }
1739  }
1740  }
1741 
1742  int offset = 5 + 2*(P-1) + 2*(Q-1) + 4*(R-1);
1743 
1744  for (i = 0; i < fid; ++i)
1745  {
1746  offset += v_GetFaceIntNcoeffs(i);
1747  }
1748 
1749  switch (fid)
1750  {
1751  case 0:
1752  for (q = 2; q <= Q; ++q)
1753  {
1754  for (p = 2; p <= P; ++p)
1755  {
1756  maparray[arrayindx[(q-2)*nummodesA+(p-2)]]
1757  = offset + (q-2)*nummodesA+(p-2);
1758  }
1759  }
1760  break;
1761 
1762  case 1:
1763  case 3:
1764  for (p = 2; p <= P; ++p)
1765  {
1766  for (r = 1; r <= R-p; ++r, ++idx)
1767  {
1768  if ((int)faceOrient == 7)
1769  {
1770  signarray[idx] = p % 2 ? -1 : 1;
1771  }
1772  maparray[idx] = offset + idx;
1773  }
1774  }
1775  break;
1776 
1777  case 2:
1778  case 4:
1779  for (q = 2; q <= Q; ++q)
1780  {
1781  for (r = 1; r <= R-q; ++r, ++idx)
1782  {
1783  if ((int)faceOrient == 7)
1784  {
1785  signarray[idx] = q % 2 ? -1 : 1;
1786  }
1787  maparray[idx] = offset + idx;
1788  }
1789  }
1790  break;
1791 
1792  default:
1793  ASSERTL0(false, "Face interior map not available.");
1794  }
1795 
1796  // Triangular faces are processed in the above switch loop; for
1797  // remaining quad faces, set up orientation if necessary.
1798  if (fid > 0)
1799  {
1800  return;
1801  }
1802 
1803  if (faceOrient == 6 || faceOrient == 8 ||
1804  faceOrient == 11 || faceOrient == 12)
1805  {
1806  if (faceOrient < 9)
1807  {
1808  for (i = 1; i < nummodesB; i += 2)
1809  {
1810  for (j = 0; j < nummodesA; j++)
1811  {
1812  signarray[arrayindx[i*nummodesA+j]] *= -1;
1813  }
1814  }
1815  }
1816  else
1817  {
1818  for (i = 0; i < nummodesB; i++)
1819  {
1820  for (j = 1; j < nummodesA; j += 2)
1821  {
1822  signarray[arrayindx[i*nummodesA+j]] *= -1;
1823  }
1824  }
1825  }
1826  }
1827 
1828  if (faceOrient == 7 || faceOrient == 8 ||
1829  faceOrient == 10 || faceOrient == 12)
1830  {
1831  if (faceOrient < 9)
1832  {
1833  for (i = 0; i < nummodesB; i++)
1834  {
1835  for (j = 1; j < nummodesA; j += 2)
1836  {
1837  signarray[arrayindx[i*nummodesA+j]] *= -1;
1838  }
1839  }
1840  }
1841  else
1842  {
1843  for (i = 1; i < nummodesB; i += 2)
1844  {
1845  for (j = 0; j < nummodesA; j++)
1846  {
1847  signarray[arrayindx[i*nummodesA+j]] *= -1;
1848  }
1849  }
1850  }
1851  }
1852  }
1853 
1854  void StdPyrExp::v_GetInteriorMap(Array<OneD, unsigned int> &outarray)
1855  {
1856  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1857  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1858  "BasisType is not a boundary interior form");
1859  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1860  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1861  "BasisType is not a boundary interior form");
1862  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1863  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1864  "BasisType is not a boundary interior form");
1865 
1866  const int nBndCoeffs = v_NumBndryCoeffs();
1867  const int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1868 
1869  if (outarray.num_elements() != nIntCoeffs)
1870  {
1871  outarray = Array<OneD, unsigned int>(nIntCoeffs);
1872  }
1873 
1874  // Loop over all interior modes.
1875  int p, idx = 0;
1876  for (p = nBndCoeffs; p < m_ncoeffs; ++p)
1877  {
1878  outarray[idx++] = p;
1879  }
1880  }
1881 
1882  void StdPyrExp::v_GetBoundaryMap(Array<OneD, unsigned int> &maparray)
1883  {
1884  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1885  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1886  "BasisType is not a boundary interior form");
1887  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1888  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1889  "BasisType is not a boundary interior form");
1890  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1891  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1892  "BasisType is not a boundary interior form");
1893 
1894  int idx = 0, nBndry = v_NumBndryCoeffs();
1895 
1896  for (idx = 0; idx < nBndry; ++idx)
1897  {
1898  maparray[idx] = idx;
1899  }
1900  }
1901 
1902  //---------------------------------------
1903  // Wrapper functions
1904  //---------------------------------------
1905 
1906  DNekMatSharedPtr StdPyrExp::v_GenMatrix(const StdMatrixKey &mkey)
1907  {
1908  return CreateGeneralMatrix(mkey);
1909  }
1910 
1911  DNekMatSharedPtr StdPyrExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
1912  {
1913  return v_GenMatrix(mkey);
1914  }
1915 
1916  /**
1917  * @brief Number tetrahedral modes in r-direction. Much the same as
1918  * StdTetExp::GetTetMode but slightly simplified since we know that the
1919  * polynomial order is the same in each direction.
1920  */
1921  int StdPyrExp::GetTetMode(const int I, const int J, const int K)
1922  {
1923  const int R = m_base[2]->GetNumModes();
1924  int i, j, cnt = 0;
1925  for (i = 0; i < I; ++i)
1926  {
1927  cnt += (R-i)*(R-i+1)/2;
1928  }
1929 
1930  i = R-I;
1931  for (j = 0; j < J; ++j)
1932  {
1933  cnt += i;
1934  i--;
1935  }
1936 
1937  return cnt + K;
1938  }
1939 
1940  void StdPyrExp::v_MultiplyByStdQuadratureMetric(
1941  const Array<OneD, const NekDouble>& inarray,
1942  Array<OneD, NekDouble>& outarray)
1943  {
1944  int i, j;
1945 
1946  int nquad0 = m_base[0]->GetNumPoints();
1947  int nquad1 = m_base[1]->GetNumPoints();
1948  int nquad2 = m_base[2]->GetNumPoints();
1949 
1950  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1951  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1952  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
1953 
1954  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1955 
1956  // Multiply by integration constants in x-direction
1957  for(i = 0; i < nquad1*nquad2; ++i)
1958  {
1959  Vmath::Vmul(nquad0, inarray.get()+i*nquad0, 1,
1960  w0.get(), 1, outarray.get()+i*nquad0,1);
1961  }
1962 
1963  // Multiply by integration constants in y-direction
1964  for(j = 0; j < nquad2; ++j)
1965  {
1966  for(i = 0; i < nquad1; ++i)
1967  {
1968  Blas::Dscal(nquad0,w1[i], &outarray[0]+i*nquad0 +
1969  j*nquad0*nquad1,1);
1970  }
1971  }
1972 
1973  // Multiply by integration constants in z-direction; need to
1974  // incorporate factor [(1-eta_3)/2]^2 into weights, but only if
1975  // using GLL quadrature points.
1976  switch(m_base[2]->GetPointsType())
1977  {
1978  // (2,0) Jacobi inner product.
1979  case LibUtilities::eGaussRadauMAlpha2Beta0:
1980  for(i = 0; i < nquad2; ++i)
1981  {
1982  Blas::Dscal(nquad0*nquad1, 0.25*w2[i],
1983  &outarray[0]+i*nquad0*nquad1, 1);
1984  }
1985  break;
1986 
1987  default:
1988  for(i = 0; i < nquad2; ++i)
1989  {
1990  Blas::Dscal(nquad0*nquad1,0.125*(1-z2[i])*(1-z2[i])*w2[i],
1991  &outarray[0]+i*nquad0*nquad1,1);
1992  }
1993  break;
1994  }
1995  }
1996  }
1997 }
Nektar::StdRegions::StdPyrExp::v_DetShapeType
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdPyrExp.cpp:889
Nektar::StdRegions::StdPyrExp::v_GetFaceInteriorMap
virtual void v_GetFaceInteriorMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
Definition: StdPyrExp.cpp:1689
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:470
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::LibUtilities::eModified_C
Principle Modified Functions .
Definition: BasisType.h:51
Nektar::StdRegions::StdPyrExp::v_GetEdgeNcoeffs
virtual int v_GetEdgeNcoeffs(const int i) const
Definition: StdPyrExp.cpp:914
Nektar::StdRegions::StdPyrExp::v_GetCoords
virtual void v_GetCoords(Array< OneD, NekDouble > &xi_x, Array< OneD, NekDouble > &xi_y, Array< OneD, NekDouble > &xi_z)
Definition: StdPyrExp.cpp:834
Nektar::StdRegions::StdPyrExp::v_MultiplyByStdQuadratureMetric
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdPyrExp.cpp:1940
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:784
Nektar::StdRegions::StdExpansion::GetBasisNumModes
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
Definition: StdExpansion.h:178
Nektar::StdRegions::StdPyrExp::v_DetFaceBasisKey
virtual const LibUtilities::BasisKey v_DetFaceBasisKey(const int i, const int k) const
Definition: StdPyrExp.cpp:996
Nektar::StdRegions::StdExpansion::MultiplyByStdQuadratureMetric
void MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:949
Nektar::StdRegions::StdPyrExp::v_FillMode
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdPyrExp.cpp:862
Nektar::StdRegions::StdPyrExp::v_GetNedges
virtual int v_GetNedges() const
Definition: StdPyrExp.cpp:879
Nektar::StdRegions::StdExpansion::GetNumPoints
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:229
Nektar::StdRegions::EvaluateTriFaceBasisKey
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Definition: StdExpansion3D.cpp:406
Nektar::StdRegions::StdPyrExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Definition: StdPyrExp.cpp:1906
Nektar::LibUtilities::eModified_A
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::StdRegions::StdPyrExp::v_GetEdgeInteriorMap
virtual void v_GetEdgeInteriorMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
Definition: StdPyrExp.cpp:1612
Nektar::LibUtilities::StdPyrData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:266
Nektar::StdRegions::StdPyrExp::v_FwdTrans
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Forward transform from physical quadrature space stored in inarray and evaluate the expansion coeffic...
Definition: StdPyrExp.cpp:499
Nektar::StdRegions::StdPyrExp::v_IProductWRTDerivBase
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdPyrExp.cpp:668
Nektar::DNekMatSharedPtr
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:700
Nektar::StdRegions::StdPyrExp::v_PhysDeriv
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Calculate the derivative of the physical points.
Definition: StdPyrExp.cpp:251
Nektar::StdRegions::StdPyrExp::v_GetFaceNcoeffs
virtual int v_GetFaceNcoeffs(const int i) const
Definition: StdPyrExp.cpp:932
Nektar::StdRegions::StdPyrExp::m_map
map< Mode, unsigned int, cmpop > m_map
Definition: StdPyrExp.h:255
Nektar::StdRegions::StdPyrExp::v_CalcNumberOfCoefficients
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdPyrExp.cpp:1034
Nektar::StdRegions::StdPyrExp::v_IProductWRTBase_SumFac
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: StdPyrExp.cpp:550
Nektar::StdRegions::StdExpansion::GetEdgeBasisType
LibUtilities::BasisType GetEdgeBasisType(const int i) const
This function returns the type of expansion basis on the i-th edge.
Definition: StdExpansion.h:413
Nektar::StdRegions::StdExpansion3D::PhysTensorDeriv
void PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d1, Array< OneD, NekDouble > &outarray_d2, Array< OneD, NekDouble > &outarray_d3)
Calculate the 3D derivative in the local tensor/collapsed coordinate at the physical points...
Definition: StdExpansion3D.cpp:68
Nektar::StdRegions::StdPyrExp::v_GetNfaces
virtual int v_GetNfaces() const
Definition: StdPyrExp.cpp:884
Nektar::StdRegions::StdExpansion
The base class for all shapes.
Definition: StdExpansion.h:69
Nektar::StdRegions::StdExpansion3D::IProductWRTBase_SumFacKernel
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
Definition: StdExpansion3D.cpp:127
Nektar::StdRegions::StdPyrExp::v_GetBoundaryMap
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdPyrExp.cpp:1882
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:51
Nektar::StdRegions::StdPyrExp::v_GetFaceIntNcoeffs
virtual int v_GetFaceIntNcoeffs(const int i) const
Definition: StdPyrExp.cpp:952
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
Nektar::StdRegions::EvaluateQuadFaceBasisKey
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Definition: StdExpansion3D.cpp:349
Nektar::StdRegions::StdPyrExp::v_LocCoordToLocCollapsed
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Definition: StdPyrExp.cpp:814
ManagerAccess.h
Nektar::StdRegions::Mode
boost::tuple< unsigned int, unsigned int, unsigned int, unsigned int > Mode
Definition: StdPyrExp.h:49
StdPyrExp.h
Nektar::StdRegions::StdPyrExp::v_BwdTrans_SumFac
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdPyrExp.cpp:388
Nektar::StdRegions::StdPyrExp::v_GetVertexMap
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdPyrExp.cpp:1603
Nektar::StdRegions::StdPyrExp::v_NumBndryCoeffs
virtual int v_NumBndryCoeffs() const
Definition: StdPyrExp.cpp:894
Nektar::StdRegions::StdExpansion::CreateGeneralMatrix
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
Definition: StdExpansion.cpp:322
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::StdRegions::StdPyrExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdPyrExp.cpp:1854
Nektar::StdRegions::StdPyrExp::v_IProductWRTBase_SumFacKernel
virtual void v_IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
Definition: StdPyrExp.cpp:579
Nektar::StdRegions::StdPyrExp::m_idxMap
map< int, map< int, map< int, pair< int, int > > > > m_idxMap
Definition: StdPyrExp.h:256
Nektar::StdRegions::StdPyrExp::v_BwdTrans_SumFacKernel
virtual void v_BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
Definition: StdPyrExp.cpp:400
Nektar::StdRegions::StdExpansion::GetBasisType
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:165
Nektar::StdRegions::StdPyrExp::v_GetFaceToElementMap
virtual void v_GetFaceToElementMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int nummodesA=-1, int nummodesB=-1)
Definition: StdPyrExp.cpp:1069
Nektar::StdRegions::StdPyrExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdPyrExp.cpp:1911
Nektar::LibUtilities::StdPyrData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:232
Nektar::iterator
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Definition: MatrixOperations.cpp:96
Nektar::StdRegions::StdExpansion::GetFaceNcoeffs
int GetFaceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th face.
Definition: StdExpansion.h:354
Nektar::StdRegions::StdPyrExp::GetTetMode
int GetTetMode(int I, int J, int K)
Number tetrahedral modes in r-direction. Much the same as StdTetExp::GetTetMode but slightly simplifi...
Definition: StdPyrExp.cpp:1921
Nektar::StdRegions::StdPyrExp::v_GetFaceNumPoints
virtual int v_GetFaceNumPoints(const int i) const
Definition: StdPyrExp.cpp:974
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:213
Nektar::LibUtilities::eGaussRadauMAlpha2Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
Nektar::StdRegions::StdPyrExp::v_GetEdgeBasisType
virtual LibUtilities::BasisType v_GetEdgeBasisType(const int i) const
Definition: StdPyrExp.cpp:1047
Nektar::StdRegions::StdPyrExp::v_IProductWRTBase
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Inner product of inarray over region with respect to the expansion basis m_base[0]->GetBdata(),m_base[1]->GetBdata(), m_base[2]->GetBdata() and return in outarray.
Definition: StdPyrExp.cpp:534
Nektar::StdRegions::find
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:315
Nektar::LibUtilities::NullBasisKey
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
Nektar::StdRegions::StdPyrExp::v_StdPhysDeriv
virtual void v_StdPhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Definition: StdPyrExp.cpp:327
Nektar::LibUtilities::BasisKey::GetNumModes
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
Nektar::LibUtilities::eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:53
Nektar::StdRegions::StdExpansion::GetPointsType
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:216
Nektar::StdRegions::StdPyrExp::v_IProductWRTDerivBase_SumFac
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdPyrExp.cpp:682
Nektar::StdRegions::StdPyrExp::v_GetNverts
virtual int v_GetNverts() const
Definition: StdPyrExp.cpp:874
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1422
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Nektar::StdRegions::StdPyrExp::v_BwdTrans
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Backward transformation is evaluated at the quadrature points.
Definition: StdPyrExp.cpp:367
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
Generated on Wed Mar 23 2016 20:14:23 for Nektar++ by   doxygen 1.8.8