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