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