Nektar++
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 nummodesA,
1075  int nummodesB)
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, p, q, r, nFaceCoeffs;
1087 
1088  int order0 = m_base[0]->GetNumModes();
1089  int order1 = m_base[1]->GetNumModes();
1090  int order2 = m_base[2]->GetNumModes();
1091 
1092  if (nummodesA == -1)
1093  {
1094  switch (fid)
1095  {
1096  case 0:
1097  nummodesA = m_base[0]->GetNumModes();
1098  nummodesB = m_base[1]->GetNumModes();
1099  break;
1100  case 1:
1101  case 3:
1102  nummodesA = m_base[0]->GetNumModes();
1103  nummodesB = m_base[2]->GetNumModes();
1104  break;
1105  case 2:
1106  case 4:
1107  nummodesA = m_base[1]->GetNumModes();
1108  nummodesB = m_base[2]->GetNumModes();
1109  break;
1110  }
1111  nFaceCoeffs = GetFaceNcoeffs(fid);
1112  }
1113  else if (fid > 0)
1114  {
1115  nFaceCoeffs = nummodesB + (nummodesA-1)*(1+2*(nummodesB-1)-(nummodesA-1))/2;
1116  }
1117  else
1118  {
1119  nFaceCoeffs = nummodesA*nummodesB;
1120  }
1121 
1122  // Allocate the map array and sign array; set sign array to ones (+)
1123  if (maparray.num_elements() != nFaceCoeffs)
1124  {
1125  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1126  }
1127 
1128  if (signarray.num_elements() != nFaceCoeffs)
1129  {
1130  signarray = Array<OneD, int>(nFaceCoeffs,1);
1131  }
1132  else
1133  {
1134  fill(signarray.get(), signarray.get() + nFaceCoeffs, 1);
1135  }
1136 
1137  // Set up an array indexing for quads, since the ordering may need
1138  // to be transposed.
1139  Array<OneD, int> arrayindx(nFaceCoeffs,-1);
1140 
1141  if (fid == 0)
1142  {
1143  for (i = 0; i < nummodesB; i++)
1144  {
1145  for (j = 0; j < nummodesA; j++)
1146  {
1147  if (faceOrient < 9)
1148  {
1149  arrayindx[i*nummodesA+j] = i*nummodesA+j;
1150  }
1151  else
1152  {
1153  arrayindx[i*nummodesA+j] = j*nummodesB+i;
1154  }
1155  }
1156  }
1157  }
1158 
1159  // Set up ordering inside each 2D face. Also for triangular faces,
1160  // populate signarray.
1161  int cnt = 0, cnt2;
1162  switch (fid)
1163  {
1164  case 0: // Bottom quad
1165 
1166  // Fill in vertices
1167  maparray[arrayindx[0]] = 0;
1168  maparray[arrayindx[1]] = 1;
1169  maparray[arrayindx[nummodesA+1]] = 2;
1170  maparray[arrayindx[nummodesA]] = 3;
1171 
1172  // Edge 0
1173  cnt = 5;
1174  for (p = 2; p < nummodesA; ++p)
1175  {
1176  maparray[arrayindx[p]] = p-2 + cnt;
1177  }
1178 
1179  // Edge 1
1180  cnt += nummodesA-2;
1181  for (q = 2; q < nummodesB; ++q)
1182  {
1183  maparray[arrayindx[q*nummodesA+1]] = q-2 + cnt;
1184  }
1185 
1186  // Edge 2
1187  cnt += nummodesB-2;
1188  for (p = 2; p < nummodesA; ++p)
1189  {
1190  maparray[arrayindx[nummodesA+p]] = p-2 + cnt;
1191  }
1192 
1193  // Edge 3
1194  cnt += nummodesA-2;
1195  for (q = 2; q < nummodesB; ++q)
1196  {
1197  maparray[arrayindx[q*nummodesA]] = q-2 + cnt;
1198  }
1199 
1200  // Interior
1201  cnt += nummodesB-2 + 4*(nummodesA-2);
1202  for (q = 2; q < nummodesB; ++q)
1203  {
1204  for (p = 2; p < nummodesA; ++p)
1205  {
1206  maparray[arrayindx[q*nummodesA+p]] = cnt + (q-2)*nummodesA+(p-2);
1207  }
1208  }
1209  break;
1210 
1211  case 1: // Left triangle
1212  // Vertices
1213  maparray[0] = 0;
1214  maparray[1] = 4;
1215  maparray[nummodesB] = 1;
1216 
1217  // Edge 0 (pyramid edge 0)
1218  cnt = 5;
1219  q = 2*nummodesB-1;
1220  for (p = 2; p < nummodesA; q += nummodesB-p, ++p)
1221  {
1222  maparray[q] = cnt++;
1223  if ((int)faceOrient == 7)
1224  {
1225  signarray[q] = p % 2 ? -1 : 1;
1226  }
1227  }
1228 
1229  // Edge 1 (pyramid edge 5)
1230  cnt = 5 + 2*(order0-2) + 2*(order1-2) + (order2-2);
1231  for (q = 2; q < nummodesB; ++q)
1232  {
1233  maparray[q] = cnt++;
1234  if ((int)faceOrient == 7)
1235  {
1236  signarray[q] = q % 2 ? -1 : 1;
1237  }
1238  }
1239 
1240  // Edge 2 (pyramid edge 4)
1241  cnt = 5 + 2*(order0-2) + 2*(order1-2);
1242  for (q = 2; q < nummodesB; ++q)
1243  {
1244  maparray[nummodesB+q-1] = cnt++;
1245  if ((int)faceOrient == 7)
1246  {
1247  signarray[nummodesB+q-1] = q % 2 ? -1 : 1;
1248  }
1249  }
1250 
1251  // Interior
1252  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1253  + v_GetFaceIntNcoeffs(0);
1254  cnt2 = 2*nummodesB + 1;
1255  for (p = 2; p < nummodesA; ++p)
1256  {
1257  for (r = 2; r < nummodesB-p; ++r)
1258  {
1259  maparray[cnt2] = cnt++;
1260  if ((int)faceOrient == 7 && p > 1)
1261  {
1262  signarray[cnt2++] = p % 2 ? -1 : 1;
1263  }
1264  }
1265  cnt2++;
1266  }
1267  break;
1268 
1269  case 2:
1270  // Vertices
1271  maparray[0] = 1;
1272  maparray[1] = 4;
1273  maparray[nummodesB] = 2;
1274 
1275  // Edge 0 (pyramid edge 1)
1276  cnt = 5 + (order0-2);
1277  q = 2*nummodesB-1;
1278  for (p = 2; p < nummodesA; q += nummodesB-p, ++p)
1279  {
1280  maparray[q] = cnt++;
1281  if ((int)faceOrient == 7)
1282  {
1283  signarray[q] = p % 2 ? -1 : 1;
1284  }
1285  }
1286 
1287  // Edge 1 (pyramid edge 6)
1288  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 2*(order2-2);
1289  for (q = 2; q < nummodesB; ++q)
1290  {
1291  maparray[q] = cnt++;
1292  if ((int)faceOrient == 7)
1293  {
1294  signarray[q] = q % 2 ? -1 : 1;
1295  }
1296  }
1297 
1298  // Edge 2 (pyramid edge 5)
1299  cnt = 5 + 2*(order0-2) + 2*(order1-2) + (order2-2);
1300  for (q = 2; q < nummodesB; ++q)
1301  {
1302  maparray[nummodesB+q-1] = cnt++;
1303  if ((int)faceOrient == 7)
1304  {
1305  signarray[nummodesB+q-1] = q % 2 ? -1 : 1;
1306  }
1307  }
1308 
1309  // Interior
1310  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1311  + v_GetFaceIntNcoeffs(0) + v_GetFaceIntNcoeffs(1);
1312  cnt2 = 2*nummodesB + 1;
1313  for (p = 2; p < nummodesA; ++p)
1314  {
1315  for (r = 2; r < nummodesB-p; ++r)
1316  {
1317  maparray[cnt2] = cnt++;
1318  if ((int)faceOrient == 7 && p > 1)
1319  {
1320  signarray[cnt2++] = p % 2 ? -1 : 1;
1321  }
1322  }
1323  cnt2++;
1324  }
1325  break;
1326 
1327  case 3: // Right triangle
1328  // Vertices
1329  maparray[0] = 3;
1330  maparray[1] = 4;
1331  maparray[nummodesB] = 2;
1332 
1333  // Edge 0 (pyramid edge 2)
1334  cnt = 5 + (order0-2) + (order1-2);
1335  q = 2*nummodesB-1;
1336  for (p = 2; p < nummodesA; q += nummodesB-p, ++p)
1337  {
1338  maparray[q] = cnt++;
1339  if ((int)faceOrient == 7)
1340  {
1341  signarray[q] = p % 2 ? -1 : 1;
1342  }
1343  }
1344 
1345  // Edge 1 (pyramid edge 6)
1346  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 2*(order2-2);
1347  for (q = 2; q < nummodesB; ++q)
1348  {
1349  maparray[q] = cnt++;
1350  if ((int)faceOrient == 7)
1351  {
1352  signarray[q] = q % 2 ? -1 : 1;
1353  }
1354  }
1355 
1356  // Edge 2 (pyramid edge 7)
1357  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 3*(order2-2);
1358  for (q = 2; q < nummodesB; ++q)
1359  {
1360  maparray[nummodesB+q-1] = cnt++;
1361  if ((int)faceOrient == 7)
1362  {
1363  signarray[nummodesB+q-1] = q % 2 ? -1 : 1;
1364  }
1365  }
1366 
1367  // Interior
1368  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1369  + v_GetFaceIntNcoeffs(0) + v_GetFaceIntNcoeffs(1)
1370  + v_GetFaceIntNcoeffs(2);
1371  cnt2 = 2*nummodesB + 1;
1372  for (p = 2; p < nummodesA; ++p)
1373  {
1374  for (r = 2; r < nummodesB-p; ++r)
1375  {
1376  maparray[cnt2] = cnt++;
1377  if ((int)faceOrient == 7 && p > 1)
1378  {
1379  signarray[cnt2++] = p % 2 ? -1 : 1;
1380  }
1381  }
1382  cnt2++;
1383  }
1384  break;
1385 
1386  case 4: // Rear quad
1387  // Vertices
1388  maparray[0] = 0;
1389  maparray[1] = 4;
1390  maparray[nummodesB] = 3;
1391 
1392  // Edge 0 (pyramid edge 3)
1393  cnt = 5 + 2*(order0-2) + (order1-2);
1394  q = 2*nummodesB-1;
1395  for (p = 2; p < nummodesA; q += nummodesB-p, ++p)
1396  {
1397  maparray[q] = cnt++;
1398  if ((int)faceOrient == 7)
1399  {
1400  signarray[q] = p % 2 ? -1 : 1;
1401  }
1402  }
1403 
1404  // Edge 1 (pyramid edge 7)
1405  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 3*(order2-2);
1406  for (q = 2; q < nummodesB; ++q)
1407  {
1408  maparray[q] = cnt++;
1409  if ((int)faceOrient == 7)
1410  {
1411  signarray[q] = q % 2 ? -1 : 1;
1412  }
1413  }
1414 
1415  // Edge 2 (pyramid edge 4)
1416  cnt = 5 + 2*(order0-2) + 2*(order1-2);
1417  for (q = 2; q < nummodesB; ++q)
1418  {
1419  maparray[nummodesB+q-1] = cnt++;
1420  if ((int)faceOrient == 7)
1421  {
1422  signarray[nummodesB+q-1] = q % 2 ? -1 : 1;
1423  }
1424  }
1425 
1426  // Interior
1427  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1428  + v_GetFaceIntNcoeffs(0) + v_GetFaceIntNcoeffs(1)
1429  + v_GetFaceIntNcoeffs(2) + v_GetFaceIntNcoeffs(3);
1430  cnt2 = 2*nummodesB + 1;
1431  for (p = 2; p < nummodesA; ++p)
1432  {
1433  for (r = 2; r < nummodesB-p; ++r)
1434  {
1435  maparray[cnt2] = cnt++;
1436  if ((int)faceOrient == 7 && p > 1)
1437  {
1438  signarray[cnt2++] = p % 2 ? -1 : 1;
1439  }
1440  }
1441  cnt2++;
1442  }
1443  break;
1444 
1445  default:
1446  ASSERTL0(false, "Face to element map unavailable.");
1447  }
1448 
1449  if (fid > 0)
1450  {
1451  // Triangles only have one possible orientation (base
1452  // direction reversed); swap edge modes.
1453  if ((int)faceOrient == 7)
1454  {
1455  swap(maparray[0], maparray[nummodesA]);
1456  for (i = 1; i < nummodesA-1; ++i)
1457  {
1458  swap(maparray[i+1], maparray[nummodesA+i]);
1459  }
1460  }
1461  }
1462  else
1463  {
1464  // The code below is exactly the same as that taken from
1465  // StdHexExp and reverses the 'b' and 'a' directions as
1466  // appropriate (1st and 2nd if statements respectively) in
1467  // quadrilateral faces.
1468  if (faceOrient == 6 || faceOrient == 8 ||
1469  faceOrient == 11 || faceOrient == 12)
1470  {
1471  if (faceOrient < 9)
1472  {
1473  for (i = 3; i < nummodesB; i += 2)
1474  {
1475  for (j = 0; j < nummodesA; j++)
1476  {
1477  signarray[arrayindx[i*nummodesA+j]] *= -1;
1478  }
1479  }
1480 
1481  for (i = 0; i < nummodesA; i++)
1482  {
1483  swap(maparray [i], maparray [i+nummodesA]);
1484  swap(signarray[i], signarray[i+nummodesA]);
1485  }
1486  }
1487  else
1488  {
1489  for (i = 0; i < nummodesB; i++)
1490  {
1491  for (j = 3; j < nummodesA; j += 2)
1492  {
1493  signarray[arrayindx[i*nummodesA+j]] *= -1;
1494  }
1495  }
1496 
1497  for (i = 0; i < nummodesB; i++)
1498  {
1499  swap (maparray [i], maparray [i+nummodesB]);
1500  swap (signarray[i], signarray[i+nummodesB]);
1501  }
1502  }
1503  }
1504 
1505  if (faceOrient == 7 || faceOrient == 8 ||
1506  faceOrient == 10 || faceOrient == 12)
1507  {
1508  if (faceOrient < 9)
1509  {
1510  for (i = 0; i < nummodesB; i++)
1511  {
1512  for (j = 3; j < nummodesA; j += 2)
1513  {
1514  signarray[arrayindx[i*nummodesA+j]] *= -1;
1515  }
1516  }
1517 
1518  for(i = 0; i < nummodesB; i++)
1519  {
1520  swap(maparray [i*nummodesA], maparray [i*nummodesA+1]);
1521  swap(signarray[i*nummodesA], signarray[i*nummodesA+1]);
1522  }
1523  }
1524  else
1525  {
1526  for (i = 3; i < nummodesB; i += 2)
1527  {
1528  for (j = 0; j < nummodesA; j++)
1529  {
1530  signarray[arrayindx[i*nummodesA+j]] *= -1;
1531  }
1532  }
1533 
1534  for (i = 0; i < nummodesA; i++)
1535  {
1536  swap(maparray [i*nummodesB], maparray [i*nummodesB+1]);
1537  swap(signarray[i*nummodesB], signarray[i*nummodesB+1]);
1538  }
1539  }
1540  }
1541  }
1542  }
1543 
1544  int StdPyrExp::v_GetVertexMap(int vId, bool useCoeffPacking)
1545  {
1546  ASSERTL0(GetEdgeBasisType(vId) == LibUtilities::eModified_A ||
1547  GetEdgeBasisType(vId) == LibUtilities::eModified_A ||
1548  GetEdgeBasisType(vId) == LibUtilities::eModified_C,
1549  "Mapping not defined for this type of basis");
1550  return vId;
1551  }
1552 
1553  void StdPyrExp::v_GetEdgeInteriorMap(
1554  const int eid,
1555  const Orientation edgeOrient,
1556  Array<OneD, unsigned int> &maparray,
1557  Array<OneD, int> &signarray)
1558  {
1559  int i;
1560  bool signChange;
1561  const int P = m_base[0]->GetNumModes() - 2;
1562  const int Q = m_base[1]->GetNumModes() - 2;
1563  const int R = m_base[2]->GetNumModes() - 2;
1564  const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid) - 2;
1565 
1566  if (maparray.num_elements() != nEdgeIntCoeffs)
1567  {
1568  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1569  }
1570 
1571  if(signarray.num_elements() != nEdgeIntCoeffs)
1572  {
1573  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1574  }
1575  else
1576  {
1577  fill(signarray.get(), signarray.get()+nEdgeIntCoeffs, 1);
1578  }
1579 
1580  // If edge is oriented backwards, change sign of modes which have
1581  // degree 2n+1, n >= 1.
1582  signChange = edgeOrient == eBackwards;
1583 
1584  int offset = 5;
1585 
1586  switch (eid)
1587  {
1588  case 0:
1589  break;
1590  case 1:
1591  offset += P;
1592  break;
1593  case 2:
1594  offset += P+Q;
1595  break;
1596  case 3:
1597  offset += 2*P+Q;
1598  break;
1599  case 4:
1600  offset += 2*(P+Q);
1601  break;
1602  case 5:
1603  offset += 2*(P+Q)+R;
1604  break;
1605  case 6:
1606  offset += 2*(P+Q+R);
1607  break;
1608  case 7:
1609  offset += 2*(P+Q)+3*R;
1610  break;
1611  default:
1612  ASSERTL0(false, "Edge not defined.");
1613  break;
1614  }
1615 
1616  for (i = 0; i < nEdgeIntCoeffs; ++i)
1617  {
1618  maparray[i] = offset + i;
1619  }
1620 
1621  if (signChange)
1622  {
1623  for (i = 1; i < nEdgeIntCoeffs; i += 2)
1624  {
1625  signarray[i] = -1;
1626  }
1627  }
1628  }
1629 
1630  void StdPyrExp::v_GetFaceInteriorMap(
1631  const int fid,
1632  const Orientation faceOrient,
1633  Array<OneD, unsigned int> &maparray,
1634  Array<OneD, int> &signarray)
1635  {
1636  const int P = m_base[0]->GetNumModes() - 1;
1637  const int Q = m_base[1]->GetNumModes() - 1;
1638  const int R = m_base[2]->GetNumModes() - 1;
1639  const int nFaceIntCoeffs = v_GetFaceIntNcoeffs(fid);
1640  int p, q, r, idx = 0;
1641  int nummodesA = 0;
1642  int nummodesB = 0;
1643  int i, j;
1644 
1645  if (maparray.num_elements() != nFaceIntCoeffs)
1646  {
1647  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1648  }
1649 
1650  if (signarray.num_elements() != nFaceIntCoeffs)
1651  {
1652  signarray = Array<OneD, int>(nFaceIntCoeffs, 1);
1653  }
1654  else
1655  {
1656  fill(signarray.get(), signarray.get()+nFaceIntCoeffs, 1);
1657  }
1658 
1659  // Set up an array indexing for quad faces, since the ordering may
1660  // need to be transposed depending on orientation.
1661  Array<OneD, int> arrayindx(nFaceIntCoeffs);
1662  if (fid == 0)
1663  {
1664  nummodesA = P-1;
1665  nummodesB = Q-1;
1666 
1667  for (i = 0; i < nummodesB; i++)
1668  {
1669  for (j = 0; j < nummodesA; j++)
1670  {
1671  if (faceOrient < 9)
1672  {
1673  arrayindx[i*nummodesA+j] = i*nummodesA+j;
1674  }
1675  else
1676  {
1677  arrayindx[i*nummodesA+j] = j*nummodesB+i;
1678  }
1679  }
1680  }
1681  }
1682 
1683  int offset = 5 + 2*(P-1) + 2*(Q-1) + 4*(R-1);
1684 
1685  for (i = 0; i < fid; ++i)
1686  {
1687  offset += v_GetFaceIntNcoeffs(i);
1688  }
1689 
1690  switch (fid)
1691  {
1692  case 0:
1693  for (q = 2; q <= Q; ++q)
1694  {
1695  for (p = 2; p <= P; ++p)
1696  {
1697  maparray[arrayindx[(q-2)*nummodesA+(p-2)]]
1698  = offset + (q-2)*nummodesA+(p-2);
1699  }
1700  }
1701  break;
1702 
1703  case 1:
1704  case 3:
1705  for (p = 2; p <= P; ++p)
1706  {
1707  for (r = 1; r <= R-p; ++r, ++idx)
1708  {
1709  if ((int)faceOrient == 7)
1710  {
1711  signarray[idx] = p % 2 ? -1 : 1;
1712  }
1713  maparray[idx] = offset + idx;
1714  }
1715  }
1716  break;
1717 
1718  case 2:
1719  case 4:
1720  for (q = 2; q <= Q; ++q)
1721  {
1722  for (r = 1; r <= R-q; ++r, ++idx)
1723  {
1724  if ((int)faceOrient == 7)
1725  {
1726  signarray[idx] = q % 2 ? -1 : 1;
1727  }
1728  maparray[idx] = offset + idx;
1729  }
1730  }
1731  break;
1732 
1733  default:
1734  ASSERTL0(false, "Face interior map not available.");
1735  }
1736 
1737  // Triangular faces are processed in the above switch loop; for
1738  // remaining quad faces, set up orientation if necessary.
1739  if (fid > 0)
1740  {
1741  return;
1742  }
1743 
1744  if (faceOrient == 6 || faceOrient == 8 ||
1745  faceOrient == 11 || faceOrient == 12)
1746  {
1747  if (faceOrient < 9)
1748  {
1749  for (i = 1; i < nummodesB; i += 2)
1750  {
1751  for (j = 0; j < nummodesA; j++)
1752  {
1753  signarray[arrayindx[i*nummodesA+j]] *= -1;
1754  }
1755  }
1756  }
1757  else
1758  {
1759  for (i = 0; i < nummodesB; i++)
1760  {
1761  for (j = 1; j < nummodesA; j += 2)
1762  {
1763  signarray[arrayindx[i*nummodesA+j]] *= -1;
1764  }
1765  }
1766  }
1767  }
1768 
1769  if (faceOrient == 7 || faceOrient == 8 ||
1770  faceOrient == 10 || faceOrient == 12)
1771  {
1772  if (faceOrient < 9)
1773  {
1774  for (i = 0; i < nummodesB; i++)
1775  {
1776  for (j = 1; j < nummodesA; j += 2)
1777  {
1778  signarray[arrayindx[i*nummodesA+j]] *= -1;
1779  }
1780  }
1781  }
1782  else
1783  {
1784  for (i = 1; i < nummodesB; i += 2)
1785  {
1786  for (j = 0; j < nummodesA; j++)
1787  {
1788  signarray[arrayindx[i*nummodesA+j]] *= -1;
1789  }
1790  }
1791  }
1792  }
1793  }
1794 
1795  void StdPyrExp::v_GetInteriorMap(Array<OneD, unsigned int> &outarray)
1796  {
1797  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1798  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1799  "BasisType is not a boundary interior form");
1800  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1801  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1802  "BasisType is not a boundary interior form");
1803  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1804  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1805  "BasisType is not a boundary interior form");
1806 
1807  const int nBndCoeffs = v_NumBndryCoeffs();
1808  const int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1809 
1810  if (outarray.num_elements() != nIntCoeffs)
1811  {
1812  outarray = Array<OneD, unsigned int>(nIntCoeffs);
1813  }
1814 
1815  // Loop over all interior modes.
1816  int p, idx = 0;
1817  for (p = nBndCoeffs; p < m_ncoeffs; ++p)
1818  {
1819  outarray[idx++] = p;
1820  }
1821  }
1822 
1823  void StdPyrExp::v_GetBoundaryMap(Array<OneD, unsigned int> &maparray)
1824  {
1825  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1826  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1827  "BasisType is not a boundary interior form");
1828  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1829  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1830  "BasisType is not a boundary interior form");
1831  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1832  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1833  "BasisType is not a boundary interior form");
1834 
1835  int idx = 0, nBndry = v_NumBndryCoeffs();
1836 
1837  for (idx = 0; idx < nBndry; ++idx)
1838  {
1839  maparray[idx] = idx;
1840  }
1841  }
1842 
1843  //---------------------------------------
1844  // Wrapper functions
1845  //---------------------------------------
1846 
1847  DNekMatSharedPtr StdPyrExp::v_GenMatrix(const StdMatrixKey &mkey)
1848  {
1849  return CreateGeneralMatrix(mkey);
1850  }
1851 
1852  DNekMatSharedPtr StdPyrExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
1853  {
1854  return v_GenMatrix(mkey);
1855  }
1856 
1857  /**
1858  * @brief Number tetrahedral modes in r-direction. Much the same as
1859  * StdTetExp::GetTetMode but slightly simplified since we know that the
1860  * polynomial order is the same in each direction.
1861  */
1862  int StdPyrExp::GetTetMode(const int I, const int J, const int K)
1863  {
1864  const int R = m_base[2]->GetNumModes();
1865  int i, j, cnt = 0;
1866  for (i = 0; i < I; ++i)
1867  {
1868  cnt += (R-i)*(R-i+1)/2;
1869  }
1870 
1871  i = R-I;
1872  for (j = 0; j < J; ++j)
1873  {
1874  cnt += i;
1875  i--;
1876  }
1877 
1878  return cnt + K;
1879  }
1880 
1881  void StdPyrExp::v_MultiplyByStdQuadratureMetric(
1882  const Array<OneD, const NekDouble>& inarray,
1883  Array<OneD, NekDouble>& outarray)
1884  {
1885  int i, j;
1886 
1887  int nquad0 = m_base[0]->GetNumPoints();
1888  int nquad1 = m_base[1]->GetNumPoints();
1889  int nquad2 = m_base[2]->GetNumPoints();
1890 
1891  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1892  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1893  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
1894 
1895  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1896 
1897  // Multiply by integration constants in x-direction
1898  for(i = 0; i < nquad1*nquad2; ++i)
1899  {
1900  Vmath::Vmul(nquad0, inarray.get()+i*nquad0, 1,
1901  w0.get(), 1, outarray.get()+i*nquad0,1);
1902  }
1903 
1904  // Multiply by integration constants in y-direction
1905  for(j = 0; j < nquad2; ++j)
1906  {
1907  for(i = 0; i < nquad1; ++i)
1908  {
1909  Blas::Dscal(nquad0,w1[i], &outarray[0]+i*nquad0 +
1910  j*nquad0*nquad1,1);
1911  }
1912  }
1913 
1914  // Multiply by integration constants in z-direction; need to
1915  // incorporate factor [(1-eta_3)/2]^2 into weights, but only if
1916  // using GLL quadrature points.
1917  switch(m_base[2]->GetPointsType())
1918  {
1919  // (2,0) Jacobi inner product.
1920  case LibUtilities::eGaussRadauMAlpha2Beta0:
1921  for(i = 0; i < nquad2; ++i)
1922  {
1923  Blas::Dscal(nquad0*nquad1, 0.25*w2[i],
1924  &outarray[0]+i*nquad0*nquad1, 1);
1925  }
1926  break;
1927 
1928  default:
1929  for(i = 0; i < nquad2; ++i)
1930  {
1931  Blas::Dscal(nquad0*nquad1,0.125*(1-z2[i])*(1-z2[i])*w2[i],
1932  &outarray[0]+i*nquad0*nquad1,1);
1933  }
1934  break;
1935  }
1936  }
1937  }
1938 }
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:1630
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:454
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:135
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:1881
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:909
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:401
Nektar::StdRegions::StdPyrExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Definition: StdPyrExp.cpp:1847
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:1553
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:684
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:397
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:1823
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:344
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:1544
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:43
Nektar::StdRegions::StdPyrExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdPyrExp.cpp:1795
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:1852
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:1862
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:187
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:312
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:165
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1356
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1016
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 Tue Dec 8 2015 23:41:30 for Nektar++ by   doxygen 1.8.9.1