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