Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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  void StdPyrExp::v_GetFaceNumModes(
872  const int fid,
873  const Orientation faceOrient,
874  int &numModes0,
875  int &numModes1)
876  {
877  int nummodes [3] = {m_base[0]->GetNumModes(),
878  m_base[1]->GetNumModes(),
879  m_base[2]->GetNumModes()};
880  switch(fid)
881  {
882  // quad
883  case 0:
884  {
885  numModes0 = nummodes[0];
886  numModes1 = nummodes[1];
887  }
888  break;
889  case 1:
890  case 3:
891  {
892  numModes0 = nummodes[0];
893  numModes1 = nummodes[2];
894  }
895  break;
896  case 2:
897  case 4:
898  {
899  numModes0 = nummodes[1];
900  numModes1 = nummodes[2];
901  }
902  break;
903  }
904 
905  if ( faceOrient >= 9 )
906  {
907  std::swap(numModes0, numModes1);
908  }
909  }
910 
911  //---------------------------------------
912  // Helper functions
913  //---------------------------------------
914 
915  int StdPyrExp::v_GetNverts() const
916  {
917  return 5;
918  }
919 
920  int StdPyrExp::v_GetNedges() const
921  {
922  return 8;
923  }
924 
925  int StdPyrExp::v_GetNfaces() const
926  {
927  return 5;
928  }
929 
930  LibUtilities::ShapeType StdPyrExp::v_DetShapeType() const
931  {
932  return LibUtilities::ePyramid;
933  }
934 
935  int StdPyrExp::v_NumBndryCoeffs() const
936  {
937  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
938  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
939  "BasisType is not a boundary interior form");
940  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
941  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
942  "BasisType is not a boundary interior form");
943  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
944  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
945  "BasisType is not a boundary interior form");
946 
947  int P = m_base[0]->GetNumModes();
948  int Q = m_base[1]->GetNumModes();
949  int R = m_base[2]->GetNumModes();
950 
951  return LibUtilities::StdPyrData::
952  getNumberOfBndCoefficients(P, Q, R);
953  }
954 
955  int StdPyrExp::v_GetEdgeNcoeffs(const int i) const
956  {
957  ASSERTL2(i >= 0 && i <= 7, "edge id is out of range");
958 
959  if (i == 0 || i == 2)
960  {
961  return GetBasisNumModes(0);
962  }
963  else if (i == 1 || i == 3)
964  {
965  return GetBasisNumModes(1);
966  }
967  else
968  {
969  return GetBasisNumModes(2);
970  }
971  }
972 
973  int StdPyrExp::v_GetFaceNcoeffs(const int i) const
974  {
975  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
976 
977  if (i == 0)
978  {
979  return GetBasisNumModes(0)*GetBasisNumModes(1);
980  }
981  else if (i == 1 || i == 3)
982  {
983  int P = GetBasisNumModes(0)-1, Q = GetBasisNumModes(2)-1;
984  return Q+1 + (P*(1 + 2*Q - P))/2;
985  }
986  else
987  {
988  int P = GetBasisNumModes(1)-1, Q = GetBasisNumModes(2)-1;
989  return Q+1 + (P*(1 + 2*Q - P))/2;
990  }
991  }
992 
993  int StdPyrExp::v_GetFaceIntNcoeffs(const int i) const
994  {
995  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
996 
997  int P = m_base[0]->GetNumModes()-1;
998  int Q = m_base[1]->GetNumModes()-1;
999  int R = m_base[2]->GetNumModes()-1;
1000 
1001  if (i == 0)
1002  {
1003  return (P-1)*(Q-1);
1004  }
1005  else if (i == 1 || i == 3)
1006  {
1007  return (P-1) * (2*(R-1) - (P-1) - 1) / 2;
1008  }
1009  else
1010  {
1011  return (Q-1) * (2*(R-1) - (Q-1) - 1) / 2;
1012  }
1013  }
1014 
1015  int StdPyrExp::v_GetFaceNumPoints(const int i) const
1016  {
1017  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
1018 
1019  if (i == 0)
1020  {
1021  return m_base[0]->GetNumPoints()*
1022  m_base[1]->GetNumPoints();
1023  }
1024  else if (i == 1 || i == 3)
1025  {
1026  return m_base[0]->GetNumPoints()*
1027  m_base[2]->GetNumPoints();
1028  }
1029  else
1030  {
1031  return m_base[1]->GetNumPoints()*
1032  m_base[2]->GetNumPoints();
1033  }
1034  }
1035 
1036 
1037  const LibUtilities::BasisKey StdPyrExp::v_DetFaceBasisKey(
1038  const int i, const int k) const
1039  {
1040  ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
1041  ASSERTL2(k >= 0 && k <= 1, "basis key id is out of range");
1042 
1043  switch(i)
1044  {
1045  case 0:
1046  {
1047  return EvaluateQuadFaceBasisKey(k,
1048  m_base[k]->GetBasisType(),
1049  m_base[k]->GetNumPoints(),
1050  m_base[k]->GetNumModes());
1051 
1052  }
1053  case 1:
1054  case 3:
1055  {
1056  return EvaluateTriFaceBasisKey(k,
1057  m_base[2*k]->GetBasisType(),
1058  m_base[2*k]->GetNumPoints(),
1059  m_base[2*k]->GetNumModes());
1060  }
1061  case 2:
1062  case 4:
1063  {
1064  return EvaluateTriFaceBasisKey(k,
1065  m_base[k+1]->GetBasisType(),
1066  m_base[k+1]->GetNumPoints(),
1067  m_base[k+1]->GetNumModes());
1068  }
1069  }
1070 
1071  // Should never get here.
1072  return LibUtilities::NullBasisKey;
1073  }
1074 
1075  int StdPyrExp::v_CalcNumberOfCoefficients(
1076  const std::vector<unsigned int> &nummodes,
1077  int &modes_offset)
1078  {
1079  int nmodes = LibUtilities::StdPyrData::getNumberOfCoefficients(
1080  nummodes[modes_offset],
1081  nummodes[modes_offset+1],
1082  nummodes[modes_offset+2]);
1083 
1084  modes_offset += 3;
1085  return nmodes;
1086  }
1087 
1088  LibUtilities::BasisType StdPyrExp::v_GetEdgeBasisType(const int i) const
1089  {
1090  ASSERTL2(i >= 0 && i <= 7, "edge id is out of range");
1091  if (i == 0 || i == 2)
1092  {
1093  return GetBasisType(0);
1094  }
1095  else if (i == 1 || i == 3)
1096  {
1097  return GetBasisType(1);
1098  }
1099  else
1100  {
1101  return GetBasisType(2);
1102  }
1103  }
1104 
1105 
1106  //---------------------------------------
1107  // Mappings
1108  //---------------------------------------
1109 
1110  void StdPyrExp::v_GetFaceToElementMap(
1111  const int fid,
1112  const Orientation faceOrient,
1113  Array<OneD, unsigned int> &maparray,
1114  Array<OneD, int> &signarray,
1115  int P,
1116  int Q)
1117  {
1118  ASSERTL1(GetEdgeBasisType(0) == GetEdgeBasisType(1),
1119  "Method only implemented if BasisType is identical"
1120  "in x and y directions");
1121  ASSERTL1(GetEdgeBasisType(0) == LibUtilities::eModified_A &&
1122  GetEdgeBasisType(4) == LibUtilities::eModified_C,
1123  "Method only implemented for Modified_A BasisType"
1124  "(x and y direction) and Modified_C BasisType (z "
1125  "direction)");
1126 
1127  int i, j, k, p, q, r, nFaceCoeffs;
1128  int nummodesA=0, nummodesB=0;
1129 
1130  int order0 = m_base[0]->GetNumModes();
1131  int order1 = m_base[1]->GetNumModes();
1132  int order2 = m_base[2]->GetNumModes();
1133 
1134  switch (fid)
1135  {
1136  case 0:
1137  nummodesA = order0;
1138  nummodesB = order1;
1139  break;
1140  case 1:
1141  case 3:
1142  nummodesA = order0;
1143  nummodesB = order2;
1144  break;
1145  case 2:
1146  case 4:
1147  nummodesA = order1;
1148  nummodesB = order2;
1149  break;
1150  }
1151 
1152  bool CheckForZeroedModes = false;
1153 
1154  if (P == -1)
1155  {
1156  P = nummodesA;
1157  Q = nummodesB;
1158  nFaceCoeffs = GetFaceNcoeffs(fid);
1159  }
1160  else if (fid > 0)
1161  {
1162  nFaceCoeffs = P*(2*Q-P+1)/2;
1163  CheckForZeroedModes = true;
1164  }
1165  else
1166  {
1167  nFaceCoeffs = P*Q;
1168  CheckForZeroedModes = true;
1169  }
1170 
1171  // Allocate the map array and sign array; set sign array to ones (+)
1172  if (maparray.num_elements() != nFaceCoeffs)
1173  {
1174  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1175  }
1176 
1177  if (signarray.num_elements() != nFaceCoeffs)
1178  {
1179  signarray = Array<OneD, int>(nFaceCoeffs,1);
1180  }
1181  else
1182  {
1183  fill(signarray.get(), signarray.get() + nFaceCoeffs, 1);
1184  }
1185 
1186  // Set up an array indexing for quads, since the ordering may need
1187  // to be transposed.
1188  Array<OneD, int> arrayindx(nFaceCoeffs,-1);
1189 
1190  if (fid == 0)
1191  {
1192  for (i = 0; i < Q; i++)
1193  {
1194  for (j = 0; j < P; j++)
1195  {
1196  if (faceOrient < 9)
1197  {
1198  arrayindx[i*P+j] = i*P+j;
1199  }
1200  else
1201  {
1202  arrayindx[i*P+j] = j*Q+i;
1203  }
1204  }
1205  }
1206  }
1207 
1208  // Set up ordering inside each 2D face. Also for triangular faces,
1209  // populate signarray.
1210  int cnt = 0, cnt2;
1211  switch (fid)
1212  {
1213  case 0: // Bottom quad
1214 
1215  // Fill in vertices
1216  maparray[arrayindx[0]] = 0;
1217  maparray[arrayindx[1]] = 1;
1218  maparray[arrayindx[P+1]] = 2;
1219  maparray[arrayindx[P]] = 3;
1220 
1221  // Edge 0
1222  cnt = 5;
1223  for (p = 2; p < P; ++p)
1224  {
1225  maparray[arrayindx[p]] = p-2 + cnt;
1226  }
1227 
1228  // Edge 1
1229  cnt += P-2;
1230  for (q = 2; q < Q; ++q)
1231  {
1232  maparray[arrayindx[q*P+1]] = q-2 + cnt;
1233  }
1234 
1235  // Edge 2
1236  cnt += Q-2;
1237  for (p = 2; p < P; ++p)
1238  {
1239  maparray[arrayindx[P+p]] = p-2 + cnt;
1240  }
1241 
1242  // Edge 3
1243  cnt += P-2;
1244  for (q = 2; q < Q; ++q)
1245  {
1246  maparray[arrayindx[q*P]] = q-2 + cnt;
1247  }
1248 
1249  // Interior
1250  cnt += Q-2 + 4*(P-2);
1251  for (q = 2; q < Q; ++q)
1252  {
1253  for (p = 2; p < P; ++p)
1254  {
1255  maparray[arrayindx[q*P+p]] = cnt + (q-2)*P+(p-2);
1256  }
1257  }
1258 
1259 
1260  break;
1261 
1262  case 1: // Left triangle
1263  // Vertices
1264  maparray[0] = 0;
1265  maparray[1] = 4;
1266  maparray[Q] = 1;
1267 
1268  // Edge 0 (pyramid edge 0)
1269  cnt = 5;
1270  q = 2*Q-1;
1271  for (p = 2; p < P; q += Q-p, ++p)
1272  {
1273  maparray[q] = cnt++;
1274  if ((int)faceOrient == 7)
1275  {
1276  signarray[q] = p % 2 ? -1 : 1;
1277  }
1278  }
1279 
1280  // Edge 1 (pyramid edge 5)
1281  cnt = 5 + 2*(order0-2) + 2*(order1-2) + (order2-2);
1282  for (q = 2; q < Q; ++q)
1283  {
1284  maparray[q] = cnt++;
1285  if ((int)faceOrient == 7)
1286  {
1287  signarray[q] = q % 2 ? -1 : 1;
1288  }
1289  }
1290 
1291  // Edge 2 (pyramid edge 4)
1292  cnt = 5 + 2*(order0-2) + 2*(order1-2);
1293  for (q = 2; q < Q; ++q)
1294  {
1295  maparray[Q+q-1] = cnt++;
1296  if ((int)faceOrient == 7)
1297  {
1298  signarray[Q+q-1] = q % 2 ? -1 : 1;
1299  }
1300  }
1301 
1302  // Interior
1303  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1304  + v_GetFaceIntNcoeffs(0);
1305  cnt2 = 2*Q + 1;
1306  for (p = 2; p < P; ++p)
1307  {
1308  for (r = 2; r < Q-p; ++r)
1309  {
1310  maparray[cnt2] = cnt++;
1311  if ((int)faceOrient == 7 && p > 1)
1312  {
1313  signarray[cnt2++] = p % 2 ? -1 : 1;
1314  }
1315  }
1316  cnt2++;
1317  }
1318  break;
1319 
1320  case 2:
1321  // Vertices
1322  maparray[0] = 1;
1323  maparray[1] = 4;
1324  maparray[Q] = 2;
1325 
1326  // Edge 0 (pyramid edge 1)
1327  cnt = 5 + (order0-2);
1328  q = 2*Q-1;
1329  for (p = 2; p < P; q += Q-p, ++p)
1330  {
1331  maparray[q] = cnt++;
1332  if ((int)faceOrient == 7)
1333  {
1334  signarray[q] = p % 2 ? -1 : 1;
1335  }
1336  }
1337 
1338  // Edge 1 (pyramid edge 6)
1339  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 2*(order2-2);
1340  for (q = 2; q < Q; ++q)
1341  {
1342  maparray[q] = cnt++;
1343  if ((int)faceOrient == 7)
1344  {
1345  signarray[q] = q % 2 ? -1 : 1;
1346  }
1347  }
1348 
1349  // Edge 2 (pyramid edge 5)
1350  cnt = 5 + 2*(order0-2) + 2*(order1-2) + (order2-2);
1351  for (q = 2; q < Q; ++q)
1352  {
1353  maparray[Q+q-1] = cnt++;
1354  if ((int)faceOrient == 7)
1355  {
1356  signarray[Q+q-1] = q % 2 ? -1 : 1;
1357  }
1358  }
1359 
1360  // Interior
1361  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1362  + v_GetFaceIntNcoeffs(0) + v_GetFaceIntNcoeffs(1);
1363  cnt2 = 2*Q + 1;
1364  for (p = 2; p < P; ++p)
1365  {
1366  for (r = 2; r < Q-p; ++r)
1367  {
1368  maparray[cnt2] = cnt++;
1369  if ((int)faceOrient == 7 && p > 1)
1370  {
1371  signarray[cnt2++] = p % 2 ? -1 : 1;
1372  }
1373  }
1374  cnt2++;
1375  }
1376  break;
1377 
1378  case 3: // Right triangle
1379  // Vertices
1380  maparray[0] = 3;
1381  maparray[1] = 4;
1382  maparray[Q] = 2;
1383 
1384  // Edge 0 (pyramid edge 2)
1385  cnt = 5 + (order0-2) + (order1-2);
1386  q = 2*Q-1;
1387  for (p = 2; p < P; q += Q-p, ++p)
1388  {
1389  maparray[q] = cnt++;
1390  if ((int)faceOrient == 7)
1391  {
1392  signarray[q] = p % 2 ? -1 : 1;
1393  }
1394  }
1395 
1396  // Edge 1 (pyramid edge 6)
1397  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 2*(order2-2);
1398  for (q = 2; q < Q; ++q)
1399  {
1400  maparray[q] = cnt++;
1401  if ((int)faceOrient == 7)
1402  {
1403  signarray[q] = q % 2 ? -1 : 1;
1404  }
1405  }
1406 
1407  // Edge 2 (pyramid edge 7)
1408  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 3*(order2-2);
1409  for (q = 2; q < Q; ++q)
1410  {
1411  maparray[Q+q-1] = cnt++;
1412  if ((int)faceOrient == 7)
1413  {
1414  signarray[Q+q-1] = q % 2 ? -1 : 1;
1415  }
1416  }
1417 
1418  // Interior
1419  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1420  + v_GetFaceIntNcoeffs(0) + v_GetFaceIntNcoeffs(1)
1421  + v_GetFaceIntNcoeffs(2);
1422  cnt2 = 2*Q + 1;
1423  for (p = 2; p < P; ++p)
1424  {
1425  for (r = 2; r < Q-p; ++r)
1426  {
1427  maparray[cnt2] = cnt++;
1428  if ((int)faceOrient == 7 && p > 1)
1429  {
1430  signarray[cnt2++] = p % 2 ? -1 : 1;
1431  }
1432  }
1433  cnt2++;
1434  }
1435  break;
1436 
1437  case 4: // Rear tri
1438  // Vertices
1439  maparray[0] = 0;
1440  maparray[1] = 4;
1441  maparray[Q] = 3;
1442 
1443  // Edge 0 (pyramid edge 3)
1444  cnt = 5 + 2*(order0-2) + (order1-2);
1445  q = 2*Q-1;
1446  for (p = 2; p < P; q += Q-p, ++p)
1447  {
1448  maparray[q] = cnt++;
1449  if ((int)faceOrient == 7)
1450  {
1451  signarray[q] = p % 2 ? -1 : 1;
1452  }
1453  }
1454 
1455  // Edge 1 (pyramid edge 7)
1456  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 3*(order2-2);
1457  for (q = 2; q < Q; ++q)
1458  {
1459  maparray[q] = cnt++;
1460  if ((int)faceOrient == 7)
1461  {
1462  signarray[q] = q % 2 ? -1 : 1;
1463  }
1464  }
1465 
1466  // Edge 2 (pyramid edge 4)
1467  cnt = 5 + 2*(order0-2) + 2*(order1-2);
1468  for (q = 2; q < Q; ++q)
1469  {
1470  maparray[Q+q-1] = cnt++;
1471  if ((int)faceOrient == 7)
1472  {
1473  signarray[Q+q-1] = q % 2 ? -1 : 1;
1474  }
1475  }
1476 
1477  // Interior
1478  cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1479  + v_GetFaceIntNcoeffs(0) + v_GetFaceIntNcoeffs(1)
1480  + v_GetFaceIntNcoeffs(2) + v_GetFaceIntNcoeffs(3);
1481  cnt2 = 2*Q + 1;
1482  for (p = 2; p < P; ++p)
1483  {
1484  for (r = 2; r < Q-p; ++r)
1485  {
1486  maparray[cnt2] = cnt++;
1487  if ((int)faceOrient == 7 && p > 1)
1488  {
1489  signarray[cnt2++] = p % 2 ? -1 : 1;
1490  }
1491  }
1492  cnt2++;
1493  }
1494  break;
1495 
1496  default:
1497  ASSERTL0(false, "Face to element map unavailable.");
1498  }
1499 
1500  if (fid > 0)
1501  {
1502 
1503  if(CheckForZeroedModes)
1504  {
1505  // zero signmap and set maparray to zero if elemental
1506  // modes are not as large as face modesl
1507  int idx = 0;
1508  for (j = 0; j < P; ++j)
1509  {
1510  idx += Q-j;
1511  for (k = Q-j; k < Q-j; ++k)
1512  {
1513  signarray[idx] = 0.0;
1514  maparray[idx++] = maparray[0];
1515  }
1516  }
1517 
1518  for (j = P; j < P; ++j)
1519  {
1520  for (k = 0; k < Q-j; ++k)
1521  {
1522  signarray[idx] = 0.0;
1523  maparray[idx++] = maparray[0];
1524  }
1525  }
1526  }
1527 
1528  // Triangles only have one possible orientation (base
1529  // direction reversed); swap edge modes.
1530  if ((int)faceOrient == 7)
1531  {
1532  swap(maparray[0], maparray[P]);
1533  for (i = 1; i < P-1; ++i)
1534  {
1535  swap(maparray[i+1], maparray[P+i]);
1536  }
1537  }
1538  }
1539  else
1540  {
1541  if(CheckForZeroedModes)
1542  {
1543  // zero signmap and set maparray to zero if elemental
1544  // modes are not as large as face modesl
1545  for (j = 0; j < P; ++j)
1546  {
1547  for (k = Q; k < Q; ++k)
1548  {
1549  signarray[arrayindx[j+k*P]] = 0.0;
1550  maparray[arrayindx[j+k*P]] = maparray[0];
1551  }
1552  }
1553 
1554  for (j = P; j < P; ++j)
1555  {
1556  for (k = 0; k < Q; ++k)
1557  {
1558  signarray[arrayindx[j+k*P]] = 0.0;
1559  maparray[arrayindx[j+k*P]] = maparray[0];
1560  }
1561  }
1562  }
1563 
1564  // The code below is exactly the same as that taken from
1565  // StdHexExp and reverses the 'b' and 'a' directions as
1566  // appropriate (1st and 2nd if statements respectively) in
1567  // quadrilateral faces.
1568  if (faceOrient == 6 || faceOrient == 8 ||
1569  faceOrient == 11 || faceOrient == 12)
1570  {
1571  if (faceOrient < 9)
1572  {
1573  for (i = 3; i < Q; i += 2)
1574  {
1575  for (j = 0; j < P; j++)
1576  {
1577  signarray[arrayindx[i*P+j]] *= -1;
1578  }
1579  }
1580 
1581  for (i = 0; i < P; i++)
1582  {
1583  swap(maparray [i], maparray [i+P]);
1584  swap(signarray[i], signarray[i+P]);
1585  }
1586  }
1587  else
1588  {
1589  for (i = 0; i < Q; i++)
1590  {
1591  for (j = 3; j < P; j += 2)
1592  {
1593  signarray[arrayindx[i*P+j]] *= -1;
1594  }
1595  }
1596 
1597  for (i = 0; i < Q; i++)
1598  {
1599  swap (maparray [i], maparray [i+Q]);
1600  swap (signarray[i], signarray[i+Q]);
1601  }
1602  }
1603  }
1604 
1605  if (faceOrient == 7 || faceOrient == 8 ||
1606  faceOrient == 10 || faceOrient == 12)
1607  {
1608  if (faceOrient < 9)
1609  {
1610  for (i = 0; i < Q; i++)
1611  {
1612  for (j = 3; j < P; j += 2)
1613  {
1614  signarray[arrayindx[i*P+j]] *= -1;
1615  }
1616  }
1617 
1618  for(i = 0; i < Q; i++)
1619  {
1620  swap(maparray [i*P], maparray [i*P+1]);
1621  swap(signarray[i*P], signarray[i*P+1]);
1622  }
1623  }
1624  else
1625  {
1626  for (i = 3; i < Q; i += 2)
1627  {
1628  for (j = 0; j < P; j++)
1629  {
1630  signarray[arrayindx[i*P+j]] *= -1;
1631  }
1632  }
1633 
1634  for (i = 0; i < P; i++)
1635  {
1636  swap(maparray [i*Q], maparray [i*Q+1]);
1637  swap(signarray[i*Q], signarray[i*Q+1]);
1638  }
1639  }
1640  }
1641  }
1642  }
1643 
1644  int StdPyrExp::v_GetVertexMap(int vId, bool useCoeffPacking)
1645  {
1646  ASSERTL0(GetEdgeBasisType(vId) == LibUtilities::eModified_A ||
1647  GetEdgeBasisType(vId) == LibUtilities::eModified_A ||
1648  GetEdgeBasisType(vId) == LibUtilities::eModified_C,
1649  "Mapping not defined for this type of basis");
1650  return vId;
1651  }
1652 
1653  void StdPyrExp::v_GetEdgeInteriorMap(
1654  const int eid,
1655  const Orientation edgeOrient,
1656  Array<OneD, unsigned int> &maparray,
1657  Array<OneD, int> &signarray)
1658  {
1659  int i;
1660  bool signChange;
1661  const int P = m_base[0]->GetNumModes() - 2;
1662  const int Q = m_base[1]->GetNumModes() - 2;
1663  const int R = m_base[2]->GetNumModes() - 2;
1664  const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid) - 2;
1665 
1666  if (maparray.num_elements() != nEdgeIntCoeffs)
1667  {
1668  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1669  }
1670 
1671  if(signarray.num_elements() != nEdgeIntCoeffs)
1672  {
1673  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1674  }
1675  else
1676  {
1677  fill(signarray.get(), signarray.get()+nEdgeIntCoeffs, 1);
1678  }
1679 
1680  // If edge is oriented backwards, change sign of modes which have
1681  // degree 2n+1, n >= 1.
1682  signChange = edgeOrient == eBackwards;
1683 
1684  int offset = 5;
1685 
1686  switch (eid)
1687  {
1688  case 0:
1689  break;
1690  case 1:
1691  offset += P;
1692  break;
1693  case 2:
1694  offset += P+Q;
1695  break;
1696  case 3:
1697  offset += 2*P+Q;
1698  break;
1699  case 4:
1700  offset += 2*(P+Q);
1701  break;
1702  case 5:
1703  offset += 2*(P+Q)+R;
1704  break;
1705  case 6:
1706  offset += 2*(P+Q+R);
1707  break;
1708  case 7:
1709  offset += 2*(P+Q)+3*R;
1710  break;
1711  default:
1712  ASSERTL0(false, "Edge not defined.");
1713  break;
1714  }
1715 
1716  for (i = 0; i < nEdgeIntCoeffs; ++i)
1717  {
1718  maparray[i] = offset + i;
1719  }
1720 
1721  if (signChange)
1722  {
1723  for (i = 1; i < nEdgeIntCoeffs; i += 2)
1724  {
1725  signarray[i] = -1;
1726  }
1727  }
1728  }
1729 
1730  void StdPyrExp::v_GetFaceInteriorMap(
1731  const int fid,
1732  const Orientation faceOrient,
1733  Array<OneD, unsigned int> &maparray,
1734  Array<OneD, int> &signarray)
1735  {
1736  const int P = m_base[0]->GetNumModes() - 1;
1737  const int Q = m_base[1]->GetNumModes() - 1;
1738  const int R = m_base[2]->GetNumModes() - 1;
1739  const int nFaceIntCoeffs = v_GetFaceIntNcoeffs(fid);
1740  int p, q, r, idx = 0;
1741  int nummodesA = 0;
1742  int nummodesB = 0;
1743  int i, j;
1744 
1745  if (maparray.num_elements() != nFaceIntCoeffs)
1746  {
1747  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1748  }
1749 
1750  if (signarray.num_elements() != nFaceIntCoeffs)
1751  {
1752  signarray = Array<OneD, int>(nFaceIntCoeffs, 1);
1753  }
1754  else
1755  {
1756  fill(signarray.get(), signarray.get()+nFaceIntCoeffs, 1);
1757  }
1758 
1759  // Set up an array indexing for quad faces, since the ordering may
1760  // need to be transposed depending on orientation.
1761  Array<OneD, int> arrayindx(nFaceIntCoeffs);
1762  if (fid == 0)
1763  {
1764  nummodesA = P-1;
1765  nummodesB = Q-1;
1766 
1767  for (i = 0; i < nummodesB; i++)
1768  {
1769  for (j = 0; j < nummodesA; j++)
1770  {
1771  if (faceOrient < 9)
1772  {
1773  arrayindx[i*nummodesA+j] = i*nummodesA+j;
1774  }
1775  else
1776  {
1777  arrayindx[i*nummodesA+j] = j*nummodesB+i;
1778  }
1779  }
1780  }
1781  }
1782 
1783  int offset = 5 + 2*(P-1) + 2*(Q-1) + 4*(R-1);
1784 
1785  for (i = 0; i < fid; ++i)
1786  {
1787  offset += v_GetFaceIntNcoeffs(i);
1788  }
1789 
1790  switch (fid)
1791  {
1792  case 0:
1793  for (q = 2; q <= Q; ++q)
1794  {
1795  for (p = 2; p <= P; ++p)
1796  {
1797  maparray[arrayindx[(q-2)*nummodesA+(p-2)]]
1798  = offset + (q-2)*nummodesA+(p-2);
1799  }
1800  }
1801  break;
1802 
1803  case 1:
1804  case 3:
1805  for (p = 2; p <= P; ++p)
1806  {
1807  for (r = 1; r <= R-p; ++r, ++idx)
1808  {
1809  if ((int)faceOrient == 7)
1810  {
1811  signarray[idx] = p % 2 ? -1 : 1;
1812  }
1813  maparray[idx] = offset + idx;
1814  }
1815  }
1816  break;
1817 
1818  case 2:
1819  case 4:
1820  for (q = 2; q <= Q; ++q)
1821  {
1822  for (r = 1; r <= R-q; ++r, ++idx)
1823  {
1824  if ((int)faceOrient == 7)
1825  {
1826  signarray[idx] = q % 2 ? -1 : 1;
1827  }
1828  maparray[idx] = offset + idx;
1829  }
1830  }
1831  break;
1832 
1833  default:
1834  ASSERTL0(false, "Face interior map not available.");
1835  }
1836 
1837  // Triangular faces are processed in the above switch loop; for
1838  // remaining quad faces, set up orientation if necessary.
1839  if (fid > 0)
1840  {
1841  return;
1842  }
1843 
1844  if (faceOrient == 6 || faceOrient == 8 ||
1845  faceOrient == 11 || faceOrient == 12)
1846  {
1847  if (faceOrient < 9)
1848  {
1849  for (i = 1; i < nummodesB; i += 2)
1850  {
1851  for (j = 0; j < nummodesA; j++)
1852  {
1853  signarray[arrayindx[i*nummodesA+j]] *= -1;
1854  }
1855  }
1856  }
1857  else
1858  {
1859  for (i = 0; i < nummodesB; i++)
1860  {
1861  for (j = 1; j < nummodesA; j += 2)
1862  {
1863  signarray[arrayindx[i*nummodesA+j]] *= -1;
1864  }
1865  }
1866  }
1867  }
1868 
1869  if (faceOrient == 7 || faceOrient == 8 ||
1870  faceOrient == 10 || faceOrient == 12)
1871  {
1872  if (faceOrient < 9)
1873  {
1874  for (i = 0; i < nummodesB; i++)
1875  {
1876  for (j = 1; j < nummodesA; j += 2)
1877  {
1878  signarray[arrayindx[i*nummodesA+j]] *= -1;
1879  }
1880  }
1881  }
1882  else
1883  {
1884  for (i = 1; i < nummodesB; i += 2)
1885  {
1886  for (j = 0; j < nummodesA; j++)
1887  {
1888  signarray[arrayindx[i*nummodesA+j]] *= -1;
1889  }
1890  }
1891  }
1892  }
1893  }
1894 
1895  void StdPyrExp::v_GetInteriorMap(Array<OneD, unsigned int> &outarray)
1896  {
1897  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1898  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1899  "BasisType is not a boundary interior form");
1900  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1901  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1902  "BasisType is not a boundary interior form");
1903  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1904  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1905  "BasisType is not a boundary interior form");
1906 
1907  const int nBndCoeffs = v_NumBndryCoeffs();
1908  const int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1909 
1910  if (outarray.num_elements() != nIntCoeffs)
1911  {
1912  outarray = Array<OneD, unsigned int>(nIntCoeffs);
1913  }
1914 
1915  // Loop over all interior modes.
1916  int p, idx = 0;
1917  for (p = nBndCoeffs; p < m_ncoeffs; ++p)
1918  {
1919  outarray[idx++] = p;
1920  }
1921  }
1922 
1923  void StdPyrExp::v_GetBoundaryMap(Array<OneD, unsigned int> &maparray)
1924  {
1925  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1926  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1927  "BasisType is not a boundary interior form");
1928  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1929  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1930  "BasisType is not a boundary interior form");
1931  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_C ||
1932  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1933  "BasisType is not a boundary interior form");
1934 
1935  int idx = 0, nBndry = v_NumBndryCoeffs();
1936 
1937  for (idx = 0; idx < nBndry; ++idx)
1938  {
1939  maparray[idx] = idx;
1940  }
1941  }
1942 
1943  //---------------------------------------
1944  // Wrapper functions
1945  //---------------------------------------
1946 
1947  DNekMatSharedPtr StdPyrExp::v_GenMatrix(const StdMatrixKey &mkey)
1948  {
1949  return CreateGeneralMatrix(mkey);
1950  }
1951 
1952  DNekMatSharedPtr StdPyrExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
1953  {
1954  return v_GenMatrix(mkey);
1955  }
1956 
1957  /**
1958  * @brief Number tetrahedral modes in r-direction. Much the same as
1959  * StdTetExp::GetTetMode but slightly simplified since we know that the
1960  * polynomial order is the same in each direction.
1961  */
1962  int StdPyrExp::GetTetMode(const int I, const int J, const int K)
1963  {
1964  const int R = m_base[2]->GetNumModes();
1965  int i, j, cnt = 0;
1966  for (i = 0; i < I; ++i)
1967  {
1968  cnt += (R-i)*(R-i+1)/2;
1969  }
1970 
1971  i = R-I;
1972  for (j = 0; j < J; ++j)
1973  {
1974  cnt += i;
1975  i--;
1976  }
1977 
1978  return cnt + K;
1979  }
1980 
1981  void StdPyrExp::v_MultiplyByStdQuadratureMetric(
1982  const Array<OneD, const NekDouble>& inarray,
1983  Array<OneD, NekDouble>& outarray)
1984  {
1985  int i, j;
1986 
1987  int nquad0 = m_base[0]->GetNumPoints();
1988  int nquad1 = m_base[1]->GetNumPoints();
1989  int nquad2 = m_base[2]->GetNumPoints();
1990 
1991  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
1992  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
1993  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
1994 
1995  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1996 
1997  // Multiply by integration constants in x-direction
1998  for(i = 0; i < nquad1*nquad2; ++i)
1999  {
2000  Vmath::Vmul(nquad0, inarray.get()+i*nquad0, 1,
2001  w0.get(), 1, outarray.get()+i*nquad0,1);
2002  }
2003 
2004  // Multiply by integration constants in y-direction
2005  for(j = 0; j < nquad2; ++j)
2006  {
2007  for(i = 0; i < nquad1; ++i)
2008  {
2009  Blas::Dscal(nquad0,w1[i], &outarray[0]+i*nquad0 +
2010  j*nquad0*nquad1,1);
2011  }
2012  }
2013 
2014  // Multiply by integration constants in z-direction; need to
2015  // incorporate factor [(1-eta_3)/2]^2 into weights, but only if
2016  // using GLL quadrature points.
2017  switch(m_base[2]->GetPointsType())
2018  {
2019  // (2,0) Jacobi inner product.
2020  case LibUtilities::eGaussRadauMAlpha2Beta0:
2021  for(i = 0; i < nquad2; ++i)
2022  {
2023  Blas::Dscal(nquad0*nquad1, 0.25*w2[i],
2024  &outarray[0]+i*nquad0*nquad1, 1);
2025  }
2026  break;
2027 
2028  default:
2029  for(i = 0; i < nquad2; ++i)
2030  {
2031  Blas::Dscal(nquad0*nquad1,0.125*(1-z2[i])*(1-z2[i])*w2[i],
2032  &outarray[0]+i*nquad0*nquad1,1);
2033  }
2034  break;
2035  }
2036  }
2037  }
2038 }
Nektar::StdRegions::StdPyrExp::v_DetShapeType
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdPyrExp.cpp:930
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:1730
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:198
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:955
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:1981
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:785
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:1037
Nektar::StdRegions::StdExpansion::MultiplyByStdQuadratureMetric
void MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:954
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:920
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:1947
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:1653
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
class_topology.P
P
Definition: class_topology.py:290
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:706
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:973
CellMLToNektar.cellml_metadata.p
tuple p
Definition: cellml_metadata.py:350
Nektar::StdRegions::StdPyrExp::v_CalcNumberOfCoefficients
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdPyrExp.cpp:1075
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:925
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:1923
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:993
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:213
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:1644
Nektar::StdRegions::StdPyrExp::v_NumBndryCoeffs
virtual int v_NumBndryCoeffs() const
Definition: StdPyrExp.cpp:935
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:1895
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_GetFaceNumModes
virtual void v_GetFaceNumModes(const int fid, const Orientation faceOrient, int &numModes0, int &numModes1)
Definition: StdPyrExp.cpp:871
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:1110
Nektar::StdRegions::StdPyrExp::m_map
std::map< Mode, unsigned int, cmpop > m_map
Definition: StdPyrExp.h:260
Nektar::StdRegions::StdPyrExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdPyrExp.cpp:1952
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:1962
Nektar::StdRegions::StdPyrExp::v_GetFaceNumPoints
virtual int v_GetFaceNumPoints(const int i) const
Definition: StdPyrExp.cpp:1015
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:250
Nektar::LibUtilities::eGaussRadauMAlpha2Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:60
Nektar::StdRegions::StdPyrExp::v_GetEdgeBasisType
virtual LibUtilities::BasisType v_GetEdgeBasisType(const int i) const
Definition: StdPyrExp.cpp:1088
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:316
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:915
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1418
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
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:261
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:299
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:183
Generated on Sat Nov 25 2017 16:18:17 for Nektar++ by   doxygen 1.8.8