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