Nektar++
StdHexExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdHexExp.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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Heaxhedral methods
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <StdRegions/StdHexExp.h>
36 
37 #ifdef max
38 #undef max
39 #endif
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45  namespace StdRegions
46  {
47 
48  StdHexExp::StdHexExp()
49  {
50  }
51 
52 
53  StdHexExp::StdHexExp(const LibUtilities::BasisKey &Ba,
54  const LibUtilities::BasisKey &Bb,
55  const LibUtilities::BasisKey &Bc):
56  StdExpansion(Ba.GetNumModes()*Bb.GetNumModes()*Bc.GetNumModes(), 3,
57  Ba, Bb, Bc),
58  StdExpansion3D(Ba.GetNumModes()*Bb.GetNumModes()*Bc.GetNumModes(),
59  Ba, Bb, Bc)
60  {
61  }
62 
63 
64  StdHexExp::StdHexExp(const StdHexExp &T):
65  StdExpansion(T),
66  StdExpansion3D(T)
67  {
68  }
69 
70 
71  StdHexExp::~StdHexExp()
72  {
73  }
74 
75  bool StdHexExp::v_IsBoundaryInteriorExpansion()
76  {
77  return
78  (m_base[0]->GetBasisType() == LibUtilities::eModified_A &&
79  m_base[1]->GetBasisType() == LibUtilities::eModified_A &&
80  m_base[2]->GetBasisType() == LibUtilities::eModified_A) ||
81  (m_base[0]->GetBasisType() == LibUtilities::eGLL_Lagrange &&
82  m_base[1]->GetBasisType() == LibUtilities::eGLL_Lagrange &&
83  m_base[1]->GetBasisType() == LibUtilities::eGLL_Lagrange);
84  }
85 
86 
87  ///////////////////////////////
88  /// Differentiation Methods
89  ///////////////////////////////
90  /**
91  * For Hexahedral region can use the PhysTensorDeriv function defined
92  * under StdExpansion. Following tenserproduct:
93  */
94  void StdHexExp::v_PhysDeriv(const Array<OneD, const NekDouble>& inarray,
95  Array<OneD, NekDouble> &out_d0,
96  Array<OneD, NekDouble> &out_d1,
97  Array<OneD, NekDouble> &out_d2)
98  {
99  StdExpansion3D::PhysTensorDeriv(inarray, out_d0, out_d1, out_d2);
100  }
101 
102 
103  /**
104  * @param dir Direction in which to compute derivative.
105  * Valid values are 0, 1, 2.
106  * @param inarray Input array.
107  * @param outarray Output array.
108  */
109  void StdHexExp::v_PhysDeriv(const int dir,
110  const Array<OneD, const NekDouble>& inarray,
111  Array<OneD, NekDouble>& outarray)
112  {
113  switch(dir)
114  {
115  case 0:
116  {
117  PhysDeriv(inarray, outarray, NullNekDouble1DArray,
118  NullNekDouble1DArray);
119  }
120  break;
121  case 1:
122  {
123  PhysDeriv(inarray, NullNekDouble1DArray, outarray,
124  NullNekDouble1DArray);
125  }
126  break;
127  case 2:
128  {
129  PhysDeriv(inarray, NullNekDouble1DArray,
130  NullNekDouble1DArray, outarray);
131  }
132  break;
133  default:
134  {
135  ASSERTL1(false,"input dir is out of range");
136  }
137  break;
138  }
139  }
140 
141  void StdHexExp::v_StdPhysDeriv(
142  const Array<OneD, const NekDouble> &inarray,
143  Array<OneD, NekDouble> &out_d0,
144  Array<OneD, NekDouble> &out_d1,
145  Array<OneD, NekDouble> &out_d2)
146  {
147  StdHexExp::v_PhysDeriv(inarray, out_d0, out_d1, out_d2);
148  }
149 
150 
151  void StdHexExp::v_StdPhysDeriv(const int dir,
152  const Array<OneD, const NekDouble>& inarray,
153  Array<OneD, NekDouble>& outarray)
154  {
155  StdHexExp::v_PhysDeriv(dir, inarray, outarray);
156  }
157 
158  /**
159  * Backward transformation is three dimensional tensorial expansion
160  * \f$ u (\xi_{1i}, \xi_{2j}, \xi_{3k})
161  * = \sum_{p=0}^{Q_x} \psi_p^a (\xi_{1i})
162  * \lbrace { \sum_{q=0}^{Q_y} \psi_{q}^a (\xi_{2j})
163  * \lbrace { \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{r}^a (\xi_{3k})
164  * \rbrace}
165  * \rbrace}. \f$
166  * And sumfactorizing step of the form is as:\\
167  * \f$ f_{r} (\xi_{3k})
168  * = \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{r}^a (\xi_{3k}),\\
169  * g_{p} (\xi_{2j}, \xi_{3k})
170  * = \sum_{r=0}^{Q_y} \psi_{p}^a (\xi_{2j}) f_{r} (\xi_{3k}),\\
171  * u(\xi_{1i}, \xi_{2j}, \xi_{3k})
172  * = \sum_{p=0}^{Q_x} \psi_{p}^a (\xi_{1i}) g_{p} (\xi_{2j}, \xi_{3k}).
173  * \f$
174  *
175  * @param inarray ?
176  * @param outarray ?
177  */
178  void StdHexExp::v_BwdTrans(
179  const Array<OneD, const NekDouble>& inarray,
180  Array<OneD, NekDouble> &outarray)
181  {
182  ASSERTL1( (m_base[1]->GetBasisType() != LibUtilities::eOrtho_B) ||
183  (m_base[1]->GetBasisType() != LibUtilities::eModified_B),
184  "Basis[1] is not a general tensor type");
185 
186  ASSERTL1( (m_base[2]->GetBasisType() != LibUtilities::eOrtho_C) ||
187  (m_base[2]->GetBasisType() != LibUtilities::eModified_C),
188  "Basis[2] is not a general tensor type");
189 
190  if(m_base[0]->Collocation() && m_base[1]->Collocation()
191  && m_base[2]->Collocation())
192  {
193  Vmath::Vcopy(m_base[0]->GetNumPoints()
194  * m_base[1]->GetNumPoints()
195  * m_base[2]->GetNumPoints(),
196  inarray, 1, outarray, 1);
197  }
198  else
199  {
200  StdHexExp::BwdTrans_SumFac(inarray,outarray);
201  }
202  }
203 
204 
205  /**
206  *
207  */
208  void StdHexExp::v_BwdTrans_SumFac(const Array<OneD, const NekDouble>& inarray,
209  Array<OneD, NekDouble> &outarray)
210  {
211  Array<OneD, NekDouble> wsp(m_base[0]->GetNumPoints()*
212  m_base[2]->GetNumModes()*
213  (m_base[1]->GetNumModes() + m_base[1]->GetNumPoints())); // FIX THIS
214 
215  BwdTrans_SumFacKernel(m_base[0]->GetBdata(),
216  m_base[1]->GetBdata(),
217  m_base[2]->GetBdata(),
218  inarray,outarray,wsp,true,true,true);
219  }
220 
221 
222  /**
223  * @param base0 x-dirn basis matrix
224  * @param base1 y-dirn basis matrix
225  * @param base2 z-dirn basis matrix
226  * @param inarray Input vector of modes.
227  * @param outarray Output vector of physical space data.
228  * @param wsp Workspace of size Q_x*P_z*(P_y+Q_y)
229  * @param doCheckCollDir0 Check for collocation of basis.
230  * @param doCheckCollDir1 Check for collocation of basis.
231  * @param doCheckCollDir2 Check for collocation of basis.
232  * @todo Account for some directions being collocated. See
233  * StdQuadExp as an example.
234  */
235  void StdHexExp::v_BwdTrans_SumFacKernel(
236  const Array<OneD, const NekDouble>& base0,
237  const Array<OneD, const NekDouble>& base1,
238  const Array<OneD, const NekDouble>& base2,
239  const Array<OneD, const NekDouble>& inarray,
240  Array<OneD, NekDouble> &outarray,
241  Array<OneD, NekDouble> &wsp,
242  bool doCheckCollDir0,
243  bool doCheckCollDir1,
244  bool doCheckCollDir2)
245  {
246  int nquad0 = m_base[0]->GetNumPoints();
247  int nquad1 = m_base[1]->GetNumPoints();
248  int nquad2 = m_base[2]->GetNumPoints();
249  int nmodes0 = m_base[0]->GetNumModes();
250  int nmodes1 = m_base[1]->GetNumModes();
251  int nmodes2 = m_base[2]->GetNumModes();
252 
253  // Check if using collocation, if requested.
254  bool colldir0 = doCheckCollDir0?(m_base[0]->Collocation()):false;
255  bool colldir1 = doCheckCollDir1?(m_base[1]->Collocation()):false;
256  bool colldir2 = doCheckCollDir2?(m_base[2]->Collocation()):false;
257 
258  // If collocation in all directions, Physical values at quadrature
259  // points is just a copy of the modes.
260  if(colldir0 && colldir1 && colldir2)
261  {
262  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,outarray.get(),1);
263  }
264  else
265  {
266  // Check sufficiently large workspace.
267  ASSERTL1(wsp.num_elements()>=nquad0*nmodes2*(nmodes1+nquad1),
268  "Workspace size is not sufficient");
269 
270  // Assign second half of workspace for 2nd DGEMM operation.
271  Array<OneD, NekDouble> wsp2 = wsp + nquad0*nmodes1*nmodes2;
272 
273  // BwdTrans in each direction using DGEMM
274  Blas::Dgemm('T','T', nmodes1*nmodes2, nquad0, nmodes0,
275  1.0, &inarray[0], nmodes0,
276  base0.get(), nquad0,
277  0.0, &wsp[0], nmodes1*nmodes2);
278  Blas::Dgemm('T','T', nquad0*nmodes2, nquad1, nmodes1,
279  1.0, &wsp[0], nmodes1,
280  base1.get(), nquad1,
281  0.0, &wsp2[0], nquad0*nmodes2);
282  Blas::Dgemm('T','T', nquad0*nquad1, nquad2, nmodes2,
283  1.0, &wsp2[0], nmodes2,
284  base2.get(), nquad2,
285  0.0, &outarray[0], nquad0*nquad1);
286  }
287  }
288 
289 
290  /**
291  * Solves the system
292  * \f$ \mathbf{B^{\top}WB\hat{u}}=\mathbf{B^{\top}Wu^{\delta}} \f$
293  *
294  * @param inarray array of physical quadrature points to be
295  * transformed, \f$ \mathbf{u^{\delta}} \f$.
296  * @param outarray array of expansion coefficients,
297  * \f$ \mathbf{\hat{u}} \f$.
298  */
299  void StdHexExp::v_FwdTrans(
300  const Array<OneD, const NekDouble>& inarray,
301  Array<OneD, NekDouble> &outarray)
302  {
303  // If using collocation expansion, coefficients match physical
304  // data points so just do a direct copy.
305  if( (m_base[0]->Collocation())
306  &&(m_base[1]->Collocation())
307  &&(m_base[2]->Collocation()) )
308  {
309  Vmath::Vcopy(GetNcoeffs(), &inarray[0], 1, &outarray[0], 1);
310  }
311  else
312  {
313  // Compute B^TWu
314  IProductWRTBase(inarray,outarray);
315 
316  // get Mass matrix inverse
317  StdMatrixKey masskey(eInvMass,DetShapeType(),*this);
318  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
319 
320  // copy inarray in case inarray == outarray
321  DNekVec in (m_ncoeffs,outarray);
322  DNekVec out(m_ncoeffs,outarray,eWrapper);
323 
324  // Solve for coefficients.
325  out = (*matsys)*in;
326 
327  }
328  }
329 
330  /**
331  * \f$
332  * \begin{array}{rcl}
333  * I_{pqr} = (\phi_{pqr}, u)_{\delta} & = &
334  * \sum_{i=0}^{nq_0} \sum_{j=0}^{nq_1} \sum_{k=0}^{nq_2}
335  * \psi_{p}^{a}(\xi_{1i}) \psi_{q}^{a}(\xi_{2j}) \psi_{r}^{a}(\xi_{3k})
336  * w_i w_j w_k u(\xi_{1,i} \xi_{2,j} \xi_{3,k})
337  *
338  * J_{i,j,k}\\ & = & \sum_{i=0}^{nq_0} \psi_p^a(\xi_{1,i})
339  * \sum_{j=0}^{nq_1} \psi_{q}^a(\xi_{2,j})
340  * \sum_{k=0}^{nq_2} \psi_{r}^a
341  * u(\xi_{1i},\xi_{2j},\xi_{3k}) J_{i,j,k}
342  * \end{array} \f$ \n
343  * where
344  * \f$ \phi_{pqr} (\xi_1 , \xi_2 , \xi_3)
345  * = \psi_p^a( \xi_1) \psi_{q}^a(\xi_2) \psi_{r}^a(\xi_3) \f$ \n
346  * which can be implemented as \n
347  * \f$f_{r} (\xi_{3k})
348  * = \sum_{k=0}^{nq_3} \psi_{r}^a u(\xi_{1i},\xi_{2j}, \xi_{3k})
349  * J_{i,j,k} = {\bf B_3 U} \f$ \n
350  * \f$ g_{q} (\xi_{3k})
351  * = \sum_{j=0}^{nq_1} \psi_{q}^a(\xi_{2j}) f_{r}(\xi_{3k})
352  * = {\bf B_2 F} \f$ \n
353  * \f$ (\phi_{pqr}, u)_{\delta}
354  * = \sum_{k=0}^{nq_0} \psi_{p}^a (\xi_{3k}) g_{q} (\xi_{3k})
355  * = {\bf B_1 G} \f$
356  *
357  * @param inarray ?
358  * @param outarray ?
359  */
360  void StdHexExp::v_IProductWRTBase(
361  const Array<OneD, const NekDouble> &inarray,
362  Array<OneD, NekDouble> &outarray)
363  {
364  if(m_base[0]->Collocation() &&
365  m_base[1]->Collocation() &&
366  m_base[2]->Collocation())
367  {
368  MultiplyByQuadratureMetric(inarray,outarray);
369  }
370  else
371  {
372  StdHexExp::v_IProductWRTBase_SumFac(inarray,outarray);
373  }
374  }
375 
376  /**
377  * Implementation of the local matrix inner product operation.
378  */
379  void StdHexExp::v_IProductWRTBase_MatOp(const Array<OneD, const NekDouble>& inarray,
380  Array<OneD, NekDouble> &outarray)
381  {
382  int nq = GetTotPoints();
383  StdMatrixKey iprodmatkey(eIProductWRTBase,DetShapeType(),*this);
384  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
385 
386  Blas::Dgemv('N',m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
387  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
388  }
389 
390  /**
391  * Implementation of the sum-factorization inner product operation.
392  */
393  void StdHexExp::v_IProductWRTBase_SumFac(
394  const Array<OneD, const NekDouble>& inarray,
395  Array<OneD, NekDouble> &outarray,
396  bool multiplybyweights)
397  {
398  int nquad0 = m_base[0]->GetNumPoints();
399  int nquad1 = m_base[1]->GetNumPoints();
400  int nquad2 = m_base[2]->GetNumPoints();
401  int order0 = m_base[0]->GetNumModes();
402  int order1 = m_base[1]->GetNumModes();
403 
404  Array<OneD, NekDouble> wsp(nquad0*nquad1*(nquad2+order0) +
405  order0*order1*nquad2);
406 
407  if(multiplybyweights)
408  {
409  Array<OneD, NekDouble> tmp(inarray.num_elements());
410  MultiplyByQuadratureMetric(inarray,tmp);
411 
412  StdHexExp::IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
413  m_base[1]->GetBdata(),
414  m_base[2]->GetBdata(),
415  tmp,outarray,wsp,true,true,true);
416  }
417  else
418  {
419  StdHexExp::IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
420  m_base[1]->GetBdata(),
421  m_base[2]->GetBdata(),
422  inarray,outarray,wsp,true,true,true);
423  }
424  }
425 
426 
427  /**
428  * Implementation of the sum-factorisation inner product operation.
429  * @todo Implement cases where only some directions are collocated.
430  */
431  void StdHexExp::v_IProductWRTBase_SumFacKernel(const Array<OneD, const NekDouble>& base0,
432  const Array<OneD, const NekDouble>& base1,
433  const Array<OneD, const NekDouble>& base2,
434  const Array<OneD, const NekDouble>& inarray,
435  Array<OneD, NekDouble> &outarray,
436  Array<OneD, NekDouble> &wsp,
437  bool doCheckCollDir0,
438  bool doCheckCollDir1,
439  bool doCheckCollDir2)
440  {
441  int nquad0 = m_base[0]->GetNumPoints();
442  int nquad1 = m_base[1]->GetNumPoints();
443  int nquad2 = m_base[2]->GetNumPoints();
444  int nmodes0 = m_base[0]->GetNumModes();
445  int nmodes1 = m_base[1]->GetNumModes();
446  int nmodes2 = m_base[2]->GetNumModes();
447 
448  bool colldir0 = doCheckCollDir0?(m_base[0]->Collocation()):false;
449  bool colldir1 = doCheckCollDir1?(m_base[1]->Collocation()):false;
450  bool colldir2 = doCheckCollDir2?(m_base[2]->Collocation()):false;
451 
452  if(colldir0 && colldir1 && colldir2)
453  {
454  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,outarray.get(),1);
455  }
456  else
457  {
458  ASSERTL1(wsp.num_elements() >= nmodes0*nquad2*(nquad1+nmodes1),
459  "Insufficient workspace size");
460 
461  Array<OneD, NekDouble> tmp0 = wsp;
462  Array<OneD, NekDouble> tmp1 = wsp + nmodes0*nquad1*nquad2;
463 
464 
465  if(colldir0)
466  {
467  // reshuffle data for next operation.
468  for(int n = 0; n < nmodes0; ++n)
469  {
470  Vmath::Vcopy(nquad1*nquad2,inarray.get()+n,nquad0,
471  tmp0.get()+nquad1*nquad2*n,1);
472  }
473  }
474  else
475  {
476  Blas::Dgemm('T', 'N', nquad1*nquad2, nmodes0, nquad0,
477  1.0, inarray.get(), nquad0,
478  base0.get(), nquad0,
479  0.0, tmp0.get(), nquad1*nquad2);
480  }
481 
482  if(colldir1)
483  {
484  // reshuffle data for next operation.
485  for(int n = 0; n < nmodes1; ++n)
486  {
487  Vmath::Vcopy(nquad2*nmodes0,tmp0.get()+n,nquad1,
488  tmp1.get()+nquad2*nmodes0*n,1);
489  }
490  }
491  else
492  {
493  Blas::Dgemm('T', 'N', nquad2*nmodes0, nmodes1, nquad1,
494  1.0, tmp0.get(), nquad1,
495  base1.get(), nquad1,
496  0.0, tmp1.get(), nquad2*nmodes0);
497  }
498 
499  if(colldir2)
500  {
501  // reshuffle data for next operation.
502  for(int n = 0; n < nmodes2; ++n)
503  {
504  Vmath::Vcopy(nmodes0*nmodes1,tmp1.get()+n,nquad2,
505  outarray.get()+nmodes0*nmodes1*n,1);
506  }
507  }
508  else
509  {
510  Blas::Dgemm('T', 'N', nmodes0*nmodes1, nmodes2, nquad2,
511  1.0, tmp1.get(), nquad2,
512  base2.get(), nquad2,
513  0.0, outarray.get(), nmodes0*nmodes1);
514  }
515  }
516  }
517 
518  void StdHexExp::v_IProductWRTDerivBase(const int dir,
519  const Array<OneD, const NekDouble>& inarray,
520  Array<OneD, NekDouble> & outarray)
521  {
522  StdHexExp::IProductWRTDerivBase_SumFac(dir,inarray,outarray);
523  }
524 
525 
526  void StdHexExp::v_IProductWRTDerivBase_MatOp(const int dir,
527  const Array<OneD, const NekDouble>& inarray,
528  Array<OneD, NekDouble> &outarray)
529  {
530  ASSERTL0((dir==0)||(dir==1)||(dir==2),"input dir is out of range");
531 
532  int nq = GetTotPoints();
533  MatrixType mtype = eIProductWRTDerivBase0;
534 
535  switch (dir)
536  {
537  case 0:
538  mtype = eIProductWRTDerivBase0;
539  break;
540  case 1:
541  mtype = eIProductWRTDerivBase1;
542  break;
543  case 2:
544  mtype = eIProductWRTDerivBase2;
545  break;
546  }
547 
548  StdMatrixKey iprodmatkey(mtype,DetShapeType(),*this);
549  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
550 
551  Blas::Dgemv('N',m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
552  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
553  }
554 
555 
556  void StdHexExp::v_IProductWRTDerivBase_SumFac(const int dir,
557  const Array<OneD, const NekDouble>& inarray,
558  Array<OneD, NekDouble> &outarray)
559  {
560  ASSERTL0((dir==0)||(dir==1)||(dir==2),"input dir is out of range");
561 
562  int nquad1 = m_base[1]->GetNumPoints();
563  int nquad2 = m_base[2]->GetNumPoints();
564  int order0 = m_base[0]->GetNumModes();
565  int order1 = m_base[1]->GetNumModes();
566 
567  // If outarray > inarray then no need for temporary storage.
568  Array<OneD, NekDouble> tmp = outarray;
569  if (outarray.num_elements() < inarray.num_elements())
570  {
571  tmp = Array<OneD, NekDouble>(inarray.num_elements());
572  }
573 
574  // Need workspace for sumfackernel though
575  Array<OneD, NekDouble> wsp(order0*nquad2*(nquad1+order1));
576 
577  // multiply by integration constants
578  MultiplyByQuadratureMetric(inarray,tmp);
579 
580  // perform sum-factorisation
581  switch (dir)
582  {
583  case 0:
584  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
585  m_base[1]->GetBdata(),
586  m_base[2]->GetBdata(),
587  tmp,outarray,wsp,
588  false,true,true);
589  break;
590  case 1:
591  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
592  m_base[1]->GetDbdata(),
593  m_base[2]->GetBdata(),
594  tmp,outarray,wsp,
595  true,false,true);
596  break;
597  case 2:
598  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
599  m_base[1]->GetBdata(),
600  m_base[2]->GetDbdata(),
601  tmp,outarray,wsp,
602  true,true,false);
603  break;
604  }
605  }
606 
607  void StdHexExp::v_LocCoordToLocCollapsed(const Array<OneD, const NekDouble>& xi,
608  Array<OneD, NekDouble>& eta)
609  {
610  eta[0] = xi[0];
611  eta[1] = xi[1];
612  eta[2] = xi[2];
613  }
614 
615  /**
616  * @note for hexahedral expansions _base[0] (i.e. p) modes run fastest.
617  */
618  void StdHexExp::v_FillMode(const int mode,
619  Array<OneD, NekDouble> &outarray)
620  {
621  int i,j;
622  int nquad0 = m_base[0]->GetNumPoints();
623  int nquad1 = m_base[1]->GetNumPoints();
624  int nquad2 = m_base[2]->GetNumPoints();
625 
626  Array<OneD, const NekDouble> base0 = m_base[0]->GetBdata();
627  Array<OneD, const NekDouble> base1 = m_base[1]->GetBdata();
628  Array<OneD, const NekDouble> base2 = m_base[2]->GetBdata();
629 
630  int btmp0 = m_base[0]->GetNumModes();
631  int btmp1 = m_base[1]->GetNumModes();
632  int mode2 = mode/(btmp0*btmp1);
633  int mode1 = (mode-mode2*btmp0*btmp1)/btmp0;
634  int mode0 = (mode-mode2*btmp0*btmp1)%btmp0;
635 
636  ASSERTL2(mode2 == (int)floor((1.0*mode)/(btmp0*btmp1)),
637  "Integer Truncation not Equiv to Floor");
638  ASSERTL2(mode1 == (int)floor((1.0*mode-mode2*btmp0*btmp1)
639  /(btmp0*btmp1)),
640  "Integer Truncation not Equiv to Floor");
641  ASSERTL2(m_ncoeffs <= mode,
642  "calling argument mode is larger than total expansion "
643  "order");
644 
645  for(i = 0; i < nquad1*nquad2; ++i)
646  {
647  Vmath::Vcopy(nquad0,(NekDouble *)(base0.get() + mode0*nquad0),1,
648  &outarray[0]+i*nquad0, 1);
649  }
650 
651  for(j = 0; j < nquad2; ++j)
652  {
653  for(i = 0; i < nquad0; ++i)
654  {
655  Vmath::Vmul(nquad1,(NekDouble *)(base1.get() + mode1*nquad1),1,
656  &outarray[0]+i+j*nquad0*nquad1, nquad0,
657  &outarray[0]+i+j*nquad0*nquad1, nquad0);
658  }
659  }
660 
661  for(i = 0; i < nquad2; i++)
662  {
663  Blas::Dscal(nquad0*nquad1,base2[mode2*nquad2+i],
664  &outarray[0]+i*nquad0*nquad1,1);
665  }
666  }
667 
668 
669  int StdHexExp::v_GetNverts() const
670  {
671  return 8;
672  }
673 
674 
675  int StdHexExp::v_GetNedges() const
676  {
677  return 12;
678  }
679 
680 
681  int StdHexExp::v_GetNfaces() const
682  {
683  return 6;
684  }
685 
686 
687  LibUtilities::ShapeType StdHexExp::v_DetShapeType() const
688  {
689  return LibUtilities::eHexahedron;
690  }
691 
692 
693  int StdHexExp::v_NumBndryCoeffs() const
694  {
695  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
696  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
697  "BasisType is not a boundary interior form");
698  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
699  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
700  "BasisType is not a boundary interior form");
701  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
702  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
703  "BasisType is not a boundary interior form");
704 
705  int nmodes0 = m_base[0]->GetNumModes();
706  int nmodes1 = m_base[1]->GetNumModes();
707  int nmodes2 = m_base[2]->GetNumModes();
708 
709  return ( 2*( nmodes0*nmodes1 + nmodes0*nmodes2
710  + nmodes1*nmodes2)
711  - 4*( nmodes0 + nmodes1 + nmodes2 ) + 8 );
712  }
713 
714  int StdHexExp::v_NumDGBndryCoeffs() const
715  {
716  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
717  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
718  "BasisType is not a boundary interior form");
719  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
720  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
721  "BasisType is not a boundary interior form");
722  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
723  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
724  "BasisType is not a boundary interior form");
725 
726  int nmodes0 = m_base[0]->GetNumModes();
727  int nmodes1 = m_base[1]->GetNumModes();
728  int nmodes2 = m_base[2]->GetNumModes();
729 
730  return 2*( nmodes0*nmodes1 + nmodes0*nmodes2
731  + nmodes1*nmodes2 );
732  }
733 
734  int StdHexExp::v_GetEdgeNcoeffs(const int i) const
735  {
736  ASSERTL2((i >= 0)&&(i <= 11),"edge id is out of range");
737 
738  if((i == 0)||(i == 2)||(i == 8)||(i == 10))
739  {
740  return GetBasisNumModes(0);
741  }
742  else if((i == 1)||(i == 3)||(i == 9)||(i == 11))
743  {
744  return GetBasisNumModes(1);
745  }
746  else
747  {
748  return GetBasisNumModes(2);
749  }
750  }
751 
752  int StdHexExp::v_GetTotalEdgeIntNcoeffs() const
753  {
754  return 4*(GetBasisNumModes(0)+GetBasisNumModes(1)+GetBasisNumModes(2));
755  }
756 
757 
758  int StdHexExp::v_GetFaceNcoeffs(const int i) const
759  {
760  ASSERTL2((i >= 0) && (i <= 5), "face id is out of range");
761  if((i == 0) || (i == 5))
762  {
763  return GetBasisNumModes(0)*GetBasisNumModes(1);
764  }
765  else if((i == 1) || (i == 3))
766  {
767  return GetBasisNumModes(0)*GetBasisNumModes(2);
768  }
769  else
770  {
771  return GetBasisNumModes(1)*GetBasisNumModes(2);
772  }
773  }
774 
775 
776  int StdHexExp::v_GetFaceIntNcoeffs(const int i) const
777  {
778  ASSERTL2((i >= 0) && (i <= 5), "face id is out of range");
779  if((i == 0) || (i == 5))
780  {
781  return (GetBasisNumModes(0)-2)*(GetBasisNumModes(1)-2);
782  }
783  else if((i == 1) || (i == 3))
784  {
785  return (GetBasisNumModes(0)-2)*(GetBasisNumModes(2)-2);
786  }
787  else
788  {
789  return (GetBasisNumModes(1)-2)*(GetBasisNumModes(2)-2);
790  }
791 
792  }
793 
794  int StdHexExp::v_GetTotalFaceIntNcoeffs() const
795  {
796  return 2*((GetBasisNumModes(0)-2)*(GetBasisNumModes(1)-2)+
797  (GetBasisNumModes(0)-2)*(GetBasisNumModes(2)-2)+
798  (GetBasisNumModes(1)-2)*(GetBasisNumModes(2)-2));
799  }
800 
801  int StdHexExp::v_GetFaceNumPoints(const int i) const
802  {
803  ASSERTL2(i >= 0 && i <= 5, "face id is out of range");
804 
805  if (i == 0 || i == 5)
806  {
807  return m_base[0]->GetNumPoints()*
808  m_base[1]->GetNumPoints();
809  }
810  else if (i == 1 || i == 3)
811  {
812  return m_base[0]->GetNumPoints()*
813  m_base[2]->GetNumPoints();
814  }
815  else
816  {
817  return m_base[1]->GetNumPoints()*
818  m_base[2]->GetNumPoints();
819  }
820  }
821 
822  LibUtilities::PointsKey StdHexExp::v_GetFacePointsKey(
823  const int i, const int j) const
824  {
825  ASSERTL2(i >= 0 && i <= 5, "face id is out of range");
826  ASSERTL2(j == 0 || j == 1, "face direction is out of range");
827 
828  if (i == 0 || i == 5)
829  {
830  return m_base[j]->GetPointsKey();
831  }
832  else if (i == 1 || i == 3)
833  {
834  return m_base[2*j]->GetPointsKey();
835  }
836  else
837  {
838  return m_base[j+1]->GetPointsKey();
839  }
840  }
841 
842  int StdHexExp::v_CalcNumberOfCoefficients(const std::vector<unsigned int> &nummodes, int &modes_offset)
843  {
844  int nmodes = nummodes[modes_offset]*nummodes[modes_offset+1]*nummodes[modes_offset+2];
845  modes_offset += 3;
846 
847  return nmodes;
848  }
849 
850 
851  const LibUtilities::BasisKey StdHexExp::v_DetFaceBasisKey(
852  const int i, const int k) const
853  {
854  ASSERTL2(i >= 0 && i <= 5, "face id is out of range");
855  ASSERTL2(k >= 0 && k <= 1, "basis key id is out of range");
856 
857  int dir = k;
858  switch(i)
859  {
860  case 0:
861  case 5:
862  dir = k;
863  break;
864  case 1:
865  case 3:
866  dir = 2*k;
867  break;
868  case 2:
869  case 4:
870  dir = k+1;
871  break;
872  }
873 
874  return EvaluateQuadFaceBasisKey(k,
875  m_base[dir]->GetBasisType(),
876  m_base[dir]->GetNumPoints(),
877  m_base[dir]->GetNumModes());
878  }
879 
880  LibUtilities::BasisType StdHexExp::v_GetEdgeBasisType(const int i) const
881  {
882  ASSERTL2((i >= 0)&&(i <= 11),"edge id is out of range");
883 
884  if((i == 0)||(i == 2)||(i==8)||(i==10))
885  {
886  return GetBasisType(0);
887  }
888  else if((i == 1)||(i == 3)||(i == 9)||(i == 11))
889  {
890  return GetBasisType(1);
891  }
892  else
893  {
894  return GetBasisType(2);
895  }
896  }
897 
898  void StdHexExp::v_GetCoords( Array<OneD, NekDouble> & xi_x,
899  Array<OneD, NekDouble> & xi_y,
900  Array<OneD, NekDouble> & xi_z)
901  {
902  Array<OneD, const NekDouble> eta_x = m_base[0]->GetZ();
903  Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
904  Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
905  int Qx = GetNumPoints(0);
906  int Qy = GetNumPoints(1);
907  int Qz = GetNumPoints(2);
908 
909  // Convert collapsed coordinates into cartesian coordinates:
910  // eta --> xi
911  for( int k = 0; k < Qz; ++k ) {
912  for( int j = 0; j < Qy; ++j ) {
913  for( int i = 0; i < Qx; ++i ) {
914  int s = i + Qx*(j + Qy*k);
915  xi_x[s] = eta_x[i];
916  xi_y[s] = eta_y[j];
917  xi_z[s] = eta_z[k];
918 
919  }
920  }
921  }
922  }
923 
924  void StdHexExp::v_GetFaceNumModes(
925  const int fid,
926  const Orientation faceOrient,
927  int &numModes0,
928  int &numModes1)
929  {
930  int nummodes [3] = {m_base[0]->GetNumModes(),
931  m_base[1]->GetNumModes(),
932  m_base[2]->GetNumModes()};
933  switch(fid)
934  {
935  case 0:
936  case 5:
937  {
938  numModes0 = nummodes[0];
939  numModes1 = nummodes[1];
940  }
941  break;
942  case 1:
943  case 3:
944  {
945  numModes0 = nummodes[0];
946  numModes1 = nummodes[2];
947  }
948  break;
949  case 2:
950  case 4:
951  {
952  numModes0 = nummodes[1];
953  numModes1 = nummodes[2];
954  }
955  break;
956  default:
957  {
958  ASSERTL0(false,"fid out of range");
959  }
960  break;
961  }
962 
963  if ( faceOrient >= eDir1FwdDir2_Dir2FwdDir1 )
964  {
965  std::swap(numModes0, numModes1);
966  }
967  }
968 
969  /**
970  * Only for basis type Modified_A or GLL_LAGRANGE in all directions.
971  */
972  void StdHexExp::v_GetFaceToElementMap(
973  const int fid,
974  const Orientation faceOrient,
975  Array<OneD, unsigned int> &maparray,
976  Array<OneD, int> &signarray,
977  int P,
978  int Q)
979  {
980  int i,j;
981  int nummodesA=0, nummodesB=0;
982 
983  ASSERTL1(GetEdgeBasisType(0) == GetEdgeBasisType(1) &&
984  GetEdgeBasisType(0) == GetEdgeBasisType(2),
985  "Method only implemented if BasisType is indentical in "
986  "all directions");
987  ASSERTL1(GetEdgeBasisType(0) == LibUtilities::eModified_A ||
988  GetEdgeBasisType(0) == LibUtilities::eGLL_Lagrange,
989  "Method only implemented for Modified_A or GLL_Lagrange BasisType");
990 
991  const int nummodes0 = m_base[0]->GetNumModes();
992  const int nummodes1 = m_base[1]->GetNumModes();
993  const int nummodes2 = m_base[2]->GetNumModes();
994 
995  switch(fid)
996  {
997  case 0:
998  case 5:
999  nummodesA = nummodes0;
1000  nummodesB = nummodes1;
1001  break;
1002  case 1:
1003  case 3:
1004  nummodesA = nummodes0;
1005  nummodesB = nummodes2;
1006  break;
1007  case 2:
1008  case 4:
1009  nummodesA = nummodes1;
1010  nummodesB = nummodes2;
1011  break;
1012  default:
1013  ASSERTL0(false,"fid must be between 0 and 5");
1014  }
1015 
1016  bool CheckForZeroedModes = false;
1017 
1018  if (P == -1)
1019  {
1020  P = nummodesA;
1021  Q = nummodesB;
1022  }
1023 
1024  if((P != nummodesA)||(Q != nummodesB))
1025  {
1026  CheckForZeroedModes = true;
1027  }
1028 
1029  bool modified = (GetEdgeBasisType(0) == LibUtilities::eModified_A);
1030  int nFaceCoeffs = P*Q;
1031 
1032  if(maparray.num_elements() != nFaceCoeffs)
1033  {
1034  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1035  }
1036 
1037  if(signarray.num_elements() != nFaceCoeffs)
1038  {
1039  signarray = Array<OneD, int>(nFaceCoeffs,1);
1040  }
1041  else
1042  {
1043  fill( signarray.get() , signarray.get()+nFaceCoeffs, 1 );
1044  }
1045 
1046  Array<OneD, int> arrayindx(nFaceCoeffs);
1047 
1048  for(i = 0; i < Q; i++)
1049  {
1050  for(j = 0; j < P; j++)
1051  {
1052  if( faceOrient < eDir1FwdDir2_Dir2FwdDir1 )
1053  {
1054  arrayindx[i*P+j] = i*P+j;
1055  }
1056  else
1057  {
1058  arrayindx[i*P+j] = j*Q+i;
1059  }
1060  }
1061  }
1062 
1063  int offset = 0;
1064  int jump1 = 1;
1065  int jump2 = 1;
1066 
1067  switch(fid)
1068  {
1069  case 5:
1070  {
1071  if (modified)
1072  {
1073  offset = nummodes0*nummodes1;
1074  }
1075  else
1076  {
1077  offset = (nummodes2-1)*nummodes0*nummodes1;
1078  jump1 = nummodes0;
1079  }
1080  }
1081  /* Falls through. */
1082  case 0:
1083  {
1084  jump1 = nummodes0;
1085  break;
1086  }
1087  case 3:
1088  {
1089  if (modified)
1090  {
1091  offset = nummodes0;
1092  }
1093  else
1094  {
1095  offset = nummodes0*(nummodes1-1);
1096  jump1 = nummodes0*nummodes1;
1097  }
1098  }
1099  /* Falls through. */
1100  case 1:
1101  {
1102  jump1 = nummodes0*nummodes1;
1103  break;
1104  }
1105  case 2:
1106  {
1107  if (modified)
1108  {
1109  offset = 1;
1110  }
1111  else
1112  {
1113  offset = nummodes0-1;
1114  jump1 = nummodes0*nummodes1;
1115  jump2 = nummodes0;
1116 
1117  }
1118  }
1119  /* Falls through. */
1120  case 4:
1121  {
1122  jump1 = nummodes0*nummodes1;
1123  jump2 = nummodes0;
1124  break;
1125  }
1126  default:
1127  ASSERTL0(false,"fid must be between 0 and 5");
1128  }
1129 
1130  for(i = 0; i < Q; i++)
1131  {
1132  for(j = 0; j < P; j++)
1133  {
1134  maparray[ arrayindx[i*P+j] ]
1135  = i*jump1 + j*jump2 + offset;
1136  }
1137  }
1138 
1139 
1140  if(CheckForZeroedModes)
1141  {
1142  if(modified)
1143  {
1144  // zero signmap and set maparray to zero if elemental
1145  // modes are not as large as face modesl
1146  for(i = 0; i < nummodesB; i++)
1147  {
1148  for(j = nummodesA; j < P; j++)
1149  {
1150  signarray[arrayindx[i*P+j]] = 0.0;
1151  maparray[arrayindx[i*P+j]] = maparray[0];
1152  }
1153  }
1154 
1155  for(i = nummodesB; i < Q; i++)
1156  {
1157  for(j = 0; j < P; j++)
1158  {
1159  signarray[arrayindx[i*P+j]] = 0.0;
1160  maparray[arrayindx[i*P+j]] = maparray[0];
1161  }
1162  }
1163  }
1164  else
1165  {
1166  ASSERTL0(false,"Different trace space face dimention and element face dimention not possible for GLL-Lagrange bases");
1167  }
1168  }
1169 
1170  if( (faceOrient==eDir1FwdDir1_Dir2BwdDir2) ||
1171  (faceOrient==eDir1BwdDir1_Dir2BwdDir2) ||
1172  (faceOrient==eDir1BwdDir2_Dir2FwdDir1) ||
1173  (faceOrient==eDir1BwdDir2_Dir2BwdDir1) )
1174  {
1175  if(faceOrient<eDir1FwdDir2_Dir2FwdDir1)
1176  {
1177  if (modified)
1178  {
1179  for(i = 3; i < Q; i+=2)
1180  {
1181  for(j = 0; j < P; j++)
1182  {
1183  signarray[ arrayindx[i*P+j] ] *= -1;
1184  }
1185  }
1186 
1187  for(i = 0; i < P; i++)
1188  {
1189  swap( maparray[i] , maparray[i+P] );
1190  swap( signarray[i] , signarray[i+P] );
1191  }
1192 
1193  }
1194  else
1195  {
1196  for(i = 0; i < P; i++)
1197  {
1198  for(j = 0; j < Q/2; j++)
1199  {
1200  swap( maparray[i + j*P],
1201  maparray[i+P*Q
1202  -P -j*P] );
1203  swap( signarray[i + j*P],
1204  signarray[i+P*Q
1205  -P -j*P]);
1206  }
1207  }
1208  }
1209  }
1210  else
1211  {
1212  if (modified)
1213  {
1214  for(i = 0; i < Q; i++)
1215  {
1216  for(j = 3; j < P; j+=2)
1217  {
1218  signarray[ arrayindx[i*P+j] ] *= -1;
1219  }
1220  }
1221 
1222  for(i = 0; i < Q; i++)
1223  {
1224  swap( maparray[i] , maparray[i+Q] );
1225  swap( signarray[i] , signarray[i+Q] );
1226  }
1227 
1228  }
1229  else
1230  {
1231  for(i = 0; i < P; i++)
1232  {
1233  for(j = 0; j < Q/2; j++)
1234  {
1235  swap( maparray[i*Q + j],
1236  maparray[i*Q + Q -1 -j]);
1237  swap( signarray[i*Q + j],
1238  signarray[i*Q + Q -1 -j]);
1239  }
1240  }
1241  }
1242  }
1243  }
1244 
1245  if( (faceOrient==eDir1BwdDir1_Dir2FwdDir2) ||
1246  (faceOrient==eDir1BwdDir1_Dir2BwdDir2) ||
1247  (faceOrient==eDir1FwdDir2_Dir2BwdDir1) ||
1248  (faceOrient==eDir1BwdDir2_Dir2BwdDir1) )
1249  {
1250  if(faceOrient<eDir1FwdDir2_Dir2FwdDir1)
1251  {
1252  if (modified)
1253  {
1254  for(i = 0; i < Q; i++)
1255  {
1256  for(j = 3; j < P; j+=2)
1257  {
1258  signarray[ arrayindx[i*P+j] ] *= -1;
1259  }
1260  }
1261 
1262  for(i = 0; i < Q; i++)
1263  {
1264  swap( maparray[i*P],
1265  maparray[i*P+1]);
1266  swap( signarray[i*P],
1267  signarray[i*P+1]);
1268  }
1269  }
1270  else
1271  {
1272  for(i = 0; i < Q; i++)
1273  {
1274  for(j = 0; j < P/2; j++)
1275  {
1276  swap( maparray[i*P + j],
1277  maparray[i*P + P -1 -j]);
1278  swap( signarray[i*P + j],
1279  signarray[i*P + P -1 -j]);
1280  }
1281  }
1282  }
1283 
1284 
1285 
1286  }
1287  else
1288  {
1289  if (modified)
1290  {
1291  for(i = 3; i < Q; i+=2)
1292  {
1293  for(j = 0; j < P; j++)
1294  {
1295  signarray[ arrayindx[i*P+j] ] *= -1;
1296  }
1297  }
1298 
1299  for(i = 0; i < P; i++)
1300  {
1301  swap( maparray[i*Q],
1302  maparray[i*Q+1]);
1303  swap( signarray[i*Q],
1304  signarray[i*Q+1]);
1305  }
1306  }
1307  else
1308  {
1309  for(i = 0; i < Q; i++)
1310  {
1311  for(j = 0; j < P/2; j++)
1312  {
1313  swap( maparray[i + j*Q] ,
1314  maparray[i+P*Q - Q -j*Q] );
1315  swap( signarray[i + j*Q] ,
1316  signarray[i+P*Q - Q -j*Q] );
1317  }
1318  }
1319  }
1320  }
1321  }
1322  }
1323 
1324 
1325 
1326  /**
1327  * Expansions in each of the three dimensions must be of type
1328  * LibUtilities#eModified_A or LibUtilities#eGLL_Lagrange.
1329  *
1330  * @param localVertexId ID of vertex (0..7)
1331  * @returns Position of vertex in local numbering scheme.
1332  */
1333  int StdHexExp::v_GetVertexMap(const int localVertexId, bool useCoeffPacking)
1334  {
1335  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1336  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1337  "BasisType is not a boundary interior form");
1338  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1339  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1340  "BasisType is not a boundary interior form");
1341  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
1342  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1343  "BasisType is not a boundary interior form");
1344 
1345  ASSERTL1((localVertexId>=0)&&(localVertexId<8),
1346  "local vertex id must be between 0 and 7");
1347 
1348  int p = 0;
1349  int q = 0;
1350  int r = 0;
1351 
1352  // Retrieve the number of modes in each dimension.
1353  int nummodes [3] = {m_base[0]->GetNumModes(),
1354  m_base[1]->GetNumModes(),
1355  m_base[2]->GetNumModes()};
1356 
1357  if(useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1358  {
1359  if(localVertexId > 3)
1360  {
1361  if( GetBasisType(2) == LibUtilities::eGLL_Lagrange)
1362  {
1363  r = nummodes[2]-1;
1364  }
1365  else
1366  {
1367  r = 1;
1368  }
1369  }
1370 
1371  switch(localVertexId % 4)
1372  {
1373  case 0:
1374  break;
1375  case 1:
1376  {
1377  if( GetBasisType(0) == LibUtilities::eGLL_Lagrange)
1378  {
1379  p = nummodes[0]-1;
1380  }
1381  else
1382  {
1383  p = 1;
1384  }
1385  }
1386  break;
1387  case 2:
1388  {
1389  if( GetBasisType(1) == LibUtilities::eGLL_Lagrange)
1390  {
1391  q = nummodes[1]-1;
1392  }
1393  else
1394  {
1395  q = 1;
1396  }
1397  }
1398  break;
1399  case 3:
1400  {
1401  if( GetBasisType(1) == LibUtilities::eGLL_Lagrange)
1402  {
1403  p = nummodes[0]-1;
1404  q = nummodes[1]-1;
1405  }
1406  else
1407  {
1408  p = 1;
1409  q = 1;
1410  }
1411  }
1412  break;
1413  }
1414  }
1415  else
1416  {
1417  // Right face (vertices 1,2,5,6)
1418  if( (localVertexId % 4) % 3 > 0 )
1419  {
1420  if( GetBasisType(0) == LibUtilities::eGLL_Lagrange)
1421  {
1422  p = nummodes[0]-1;
1423  }
1424  else
1425  {
1426  p = 1;
1427  }
1428  }
1429 
1430  // Back face (vertices 2,3,6,7)
1431  if( localVertexId % 4 > 1 )
1432  {
1433  if( GetBasisType(1) == LibUtilities::eGLL_Lagrange)
1434  {
1435  q = nummodes[1]-1;
1436  }
1437  else
1438  {
1439  q = 1;
1440  }
1441  }
1442 
1443  // Top face (vertices 4,5,6,7)
1444  if( localVertexId > 3)
1445  {
1446  if( GetBasisType(2) == LibUtilities::eGLL_Lagrange)
1447  {
1448  r = nummodes[2]-1;
1449  }
1450  else
1451  {
1452  r = 1;
1453  }
1454  }
1455  }
1456  // Compute the local number.
1457  return r*nummodes[0]*nummodes[1] + q*nummodes[0] + p;
1458  }
1459 
1460 
1461  /**
1462  * @param eid The edge to compute the numbering for.
1463  * @param edgeOrient Orientation of the edge.
1464  * @param maparray Storage for computed mapping array.
1465  * @param signarray ?
1466  */
1467  void StdHexExp::v_GetEdgeInteriorMap(const int eid,
1468  const Orientation edgeOrient,
1469  Array<OneD, unsigned int> &maparray,
1470  Array<OneD, int> &signarray)
1471  {
1472  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1473  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1474  "BasisType is not a boundary interior form");
1475  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1476  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1477  "BasisType is not a boundary interior form");
1478  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
1479  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1480  "BasisType is not a boundary interior form");
1481 
1482  ASSERTL1((eid>=0)&&(eid<12),
1483  "local edge id must be between 0 and 11");
1484 
1485  int nEdgeIntCoeffs = GetEdgeNcoeffs(eid)-2;
1486 
1487  if(maparray.num_elements()!=nEdgeIntCoeffs)
1488  {
1489  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1490  }
1491 
1492  if(signarray.num_elements() != nEdgeIntCoeffs)
1493  {
1494  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1495  }
1496  else
1497  {
1498  fill( signarray.get() , signarray.get()+nEdgeIntCoeffs, 1 );
1499  }
1500 
1501  int nummodes [3] = {m_base[0]->GetNumModes(),
1502  m_base[1]->GetNumModes(),
1503  m_base[2]->GetNumModes()};
1504 
1505  const LibUtilities::BasisType bType [3] = {GetBasisType(0),
1506  GetBasisType(1),
1507  GetBasisType(2)};
1508 
1509  bool reverseOrdering = false;
1510  bool signChange = false;
1511 
1512  int IdxRange [3][2] = {{0,0},{0,0},{0,0}};
1513 
1514  switch(eid)
1515  {
1516  case 0:
1517  case 1:
1518  case 2:
1519  case 3:
1520  {
1521  IdxRange[2][0] = 0;
1522  IdxRange[2][1] = 1;
1523  }
1524  break;
1525  case 8:
1526  case 9:
1527  case 10:
1528  case 11:
1529  {
1530  if( bType[2] == LibUtilities::eGLL_Lagrange)
1531  {
1532  IdxRange[2][0] = nummodes[2] - 1;
1533  IdxRange[2][1] = nummodes[2];
1534  }
1535  else
1536  {
1537  IdxRange[2][0] = 1;
1538  IdxRange[2][1] = 2;
1539  }
1540  }
1541  break;
1542  case 4:
1543  case 5:
1544  case 6:
1545  case 7:
1546  {
1547  if( bType[2] == LibUtilities::eGLL_Lagrange)
1548  {
1549  IdxRange[2][0] = 1;
1550  IdxRange[2][1] = nummodes[2] - 1;
1551 
1552  if(edgeOrient==eBackwards)
1553  {
1554  reverseOrdering = true;
1555  }
1556  }
1557  else
1558  {
1559  IdxRange[2][0] = 2;
1560  IdxRange[2][1] = nummodes[2];
1561 
1562  if(edgeOrient==eBackwards)
1563  {
1564  signChange = true;
1565  }
1566  }
1567  }
1568  break;
1569  }
1570 
1571  switch(eid)
1572  {
1573  case 0:
1574  case 4:
1575  case 5:
1576  case 8:
1577  {
1578  IdxRange[1][0] = 0;
1579  IdxRange[1][1] = 1;
1580  }
1581  break;
1582  case 2:
1583  case 6:
1584  case 7:
1585  case 10:
1586  {
1587  if( bType[1] == LibUtilities::eGLL_Lagrange)
1588  {
1589  IdxRange[1][0] = nummodes[1] - 1;
1590  IdxRange[1][1] = nummodes[1];
1591  }
1592  else
1593  {
1594  IdxRange[1][0] = 1;
1595  IdxRange[1][1] = 2;
1596  }
1597  }
1598  break;
1599  case 1:
1600  case 9:
1601  {
1602  if( bType[1] == LibUtilities::eGLL_Lagrange)
1603  {
1604  IdxRange[1][0] = 1;
1605  IdxRange[1][1] = nummodes[1] - 1;
1606 
1607  if(edgeOrient==eBackwards)
1608  {
1609  reverseOrdering = true;
1610  }
1611  }
1612  else
1613  {
1614  IdxRange[1][0] = 2;
1615  IdxRange[1][1] = nummodes[1];
1616 
1617  if(edgeOrient==eBackwards)
1618  {
1619  signChange = true;
1620  }
1621  }
1622  }
1623  break;
1624  case 3:
1625  case 11:
1626  {
1627  if( bType[1] == LibUtilities::eGLL_Lagrange)
1628  {
1629  IdxRange[1][0] = 1;
1630  IdxRange[1][1] = nummodes[1] - 1;
1631 
1632  if(edgeOrient==eForwards)
1633  {
1634  reverseOrdering = true;
1635  }
1636  }
1637  else
1638  {
1639  IdxRange[1][0] = 2;
1640  IdxRange[1][1] = nummodes[1];
1641 
1642  if(edgeOrient==eForwards)
1643  {
1644  signChange = true;
1645  }
1646  }
1647  }
1648  break;
1649  }
1650 
1651  switch(eid)
1652  {
1653  case 3:
1654  case 4:
1655  case 7:
1656  case 11:
1657  {
1658  IdxRange[0][0] = 0;
1659  IdxRange[0][1] = 1;
1660  }
1661  break;
1662  case 1:
1663  case 5:
1664  case 6:
1665  case 9:
1666  {
1667  if( bType[0] == LibUtilities::eGLL_Lagrange)
1668  {
1669  IdxRange[0][0] = nummodes[0] - 1;
1670  IdxRange[0][1] = nummodes[0];
1671  }
1672  else
1673  {
1674  IdxRange[0][0] = 1;
1675  IdxRange[0][1] = 2;
1676  }
1677  }
1678  break;
1679  case 0:
1680  case 8:
1681  {
1682  if( bType[0] == LibUtilities::eGLL_Lagrange)
1683  {
1684  IdxRange[0][0] = 1;
1685  IdxRange[0][1] = nummodes[0] - 1;
1686 
1687  if(edgeOrient==eBackwards)
1688  {
1689  reverseOrdering = true;
1690  }
1691  }
1692  else
1693  {
1694  IdxRange[0][0] = 2;
1695  IdxRange[0][1] = nummodes[0];
1696 
1697  if(edgeOrient==eBackwards)
1698  {
1699  signChange = true;
1700  }
1701  }
1702  }
1703  break;
1704  case 2:
1705  case 10:
1706  {
1707  if( bType[0] == LibUtilities::eGLL_Lagrange)
1708  {
1709  IdxRange[0][0] = 1;
1710  IdxRange[0][1] = nummodes[0] - 1;
1711 
1712  if(edgeOrient==eForwards)
1713  {
1714  reverseOrdering = true;
1715  }
1716  }
1717  else
1718  {
1719  IdxRange[0][0] = 2;
1720  IdxRange[0][1] = nummodes[0];
1721 
1722  if(edgeOrient==eForwards)
1723  {
1724  signChange = true;
1725  }
1726  }
1727  }
1728  break;
1729  }
1730 
1731  int p,q,r;
1732  int cnt = 0;
1733 
1734  for(r = IdxRange[2][0]; r < IdxRange[2][1]; r++)
1735  {
1736  for(q = IdxRange[1][0]; q < IdxRange[1][1]; q++)
1737  {
1738  for(p = IdxRange[0][0]; p < IdxRange[0][1]; p++)
1739  {
1740  maparray[cnt++]
1741  = r*nummodes[0]*nummodes[1] + q*nummodes[0] + p;
1742  }
1743  }
1744  }
1745 
1746  if( reverseOrdering )
1747  {
1748  reverse( maparray.get() , maparray.get()+nEdgeIntCoeffs );
1749  }
1750 
1751  if( signChange )
1752  {
1753  for(p = 1; p < nEdgeIntCoeffs; p+=2)
1754  {
1755  signarray[p] = -1;
1756  }
1757  }
1758  }
1759 
1760 
1761  /**
1762  * Generate mapping describing which elemental modes lie on the
1763  * interior of a given face. Accounts for face orientation.
1764  */
1765  void StdHexExp::v_GetFaceInteriorMap(const int fid,
1766  const Orientation faceOrient,
1767  Array<OneD, unsigned int> &maparray,
1768  Array<OneD, int>& signarray)
1769  {
1770  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1771  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1772  "BasisType is not a boundary interior form");
1773  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1774  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1775  "BasisType is not a boundary interior form");
1776  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
1777  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1778  "BasisType is not a boundary interior form");
1779 
1780  ASSERTL1((fid>=0)&&(fid<6),
1781  "local face id must be between 0 and 5");
1782 
1783  int nFaceIntCoeffs = GetFaceIntNcoeffs(fid);
1784 
1785  if(maparray.num_elements()!=nFaceIntCoeffs)
1786  {
1787  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1788  }
1789 
1790  if(signarray.num_elements() != nFaceIntCoeffs)
1791  {
1792  signarray = Array<OneD, int>(nFaceIntCoeffs,1);
1793  }
1794  else
1795  {
1796  fill( signarray.get() , signarray.get()+nFaceIntCoeffs, 1 );
1797  }
1798 
1799  int nummodes [3] = {m_base[0]->GetNumModes(),
1800  m_base[1]->GetNumModes(),
1801  m_base[2]->GetNumModes()};
1802 
1803  const LibUtilities::BasisType bType [3] = {GetBasisType(0),
1804  GetBasisType(1),
1805  GetBasisType(2)};
1806 
1807  int nummodesA = 0;
1808  int nummodesB = 0;
1809 
1810  // Determine the number of modes in face directions A & B based
1811  // on the face index given.
1812  switch(fid)
1813  {
1814  case 0:
1815  case 5:
1816  {
1817  nummodesA = nummodes[0];
1818  nummodesB = nummodes[1];
1819  }
1820  break;
1821  case 1:
1822  case 3:
1823  {
1824  nummodesA = nummodes[0];
1825  nummodesB = nummodes[2];
1826  }
1827  break;
1828  case 2:
1829  case 4:
1830  {
1831  nummodesA = nummodes[1];
1832  nummodesB = nummodes[2];
1833  }
1834  }
1835 
1836  int i,j;
1837  Array<OneD, int> arrayindx(nFaceIntCoeffs);
1838 
1839  // Create a mapping array to account for transposition of the
1840  // coordinates due to face orientation.
1841  for(i = 0; i < (nummodesB-2); i++)
1842  {
1843  for(j = 0; j < (nummodesA-2); j++)
1844  {
1845  if( faceOrient < eDir1FwdDir2_Dir2FwdDir1 )
1846  {
1847  arrayindx[i*(nummodesA-2)+j] = i*(nummodesA-2)+j;
1848  }
1849  else
1850  {
1851  arrayindx[i*(nummodesA-2)+j] = j*(nummodesB-2)+i;
1852  }
1853  }
1854  }
1855 
1856  int IdxRange [3][2];
1857  int Incr[3];
1858 
1859  Array<OneD, int> sign0(nummodes[0], 1);
1860  Array<OneD, int> sign1(nummodes[1], 1);
1861  Array<OneD, int> sign2(nummodes[2], 1);
1862 
1863  // Set the upper and lower bounds, and increment for the faces
1864  // involving the first coordinate direction.
1865  switch(fid)
1866  {
1867  case 0: // bottom face
1868  {
1869  IdxRange[2][0] = 0;
1870  IdxRange[2][1] = 1;
1871  Incr[2] = 1;
1872  }
1873  break;
1874  case 5: // top face
1875  {
1876  if( bType[2] == LibUtilities::eGLL_Lagrange)
1877  {
1878  IdxRange[2][0] = nummodes[2] - 1;
1879  IdxRange[2][1] = nummodes[2];
1880  Incr[2] = 1;
1881  }
1882  else
1883  {
1884  IdxRange[2][0] = 1;
1885  IdxRange[2][1] = 2;
1886  Incr[2] = 1;
1887  }
1888 
1889  }
1890  break;
1891  default: // all other faces
1892  {
1893  if( bType[2] == LibUtilities::eGLL_Lagrange)
1894  {
1895  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 2 )
1896  {
1897  IdxRange[2][0] = nummodes[2] - 2;
1898  IdxRange[2][1] = 0;
1899  Incr[2] = -1;
1900 
1901  }
1902  else
1903  {
1904  IdxRange[2][0] = 1;
1905  IdxRange[2][1] = nummodes[2] - 1;
1906  Incr[2] = 1;
1907  }
1908  }
1909  else
1910  {
1911  IdxRange[2][0] = 2;
1912  IdxRange[2][1] = nummodes[2];
1913  Incr[2] = 1;
1914 
1915  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 2 )
1916  {
1917  for(i = 3; i < nummodes[2]; i+=2)
1918  {
1919  sign2[i] = -1;
1920  }
1921  }
1922  }
1923  }
1924  }
1925 
1926  // Set the upper and lower bounds, and increment for the faces
1927  // involving the second coordinate direction.
1928  switch(fid)
1929  {
1930  case 1:
1931  {
1932  IdxRange[1][0] = 0;
1933  IdxRange[1][1] = 1;
1934  Incr[1] = 1;
1935  }
1936  break;
1937  case 3:
1938  {
1939  if( bType[1] == LibUtilities::eGLL_Lagrange)
1940  {
1941  IdxRange[1][0] = nummodes[1] - 1;
1942  IdxRange[1][1] = nummodes[1];
1943  Incr[1] = 1;
1944  }
1945  else
1946  {
1947  IdxRange[1][0] = 1;
1948  IdxRange[1][1] = 2;
1949  Incr[1] = 1;
1950  }
1951  }
1952  break;
1953  case 0:
1954  case 5:
1955  {
1956  if( bType[1] == LibUtilities::eGLL_Lagrange)
1957  {
1958  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 2 )
1959  {
1960  IdxRange[1][0] = nummodes[1] - 2;
1961  IdxRange[1][1] = 0;
1962  Incr[1] = -1;
1963 
1964  }
1965  else
1966  {
1967  IdxRange[1][0] = 1;
1968  IdxRange[1][1] = nummodes[1] - 1;
1969  Incr[1] = 1;
1970  }
1971  }
1972  else
1973  {
1974  IdxRange[1][0] = 2;
1975  IdxRange[1][1] = nummodes[1];
1976  Incr[1] = 1;
1977 
1978  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 2 )
1979  {
1980  for(i = 3; i < nummodes[1]; i+=2)
1981  {
1982  sign1[i] = -1;
1983  }
1984  }
1985  }
1986  }
1987  break;
1988  default: // case2: case4:
1989  {
1990  if( bType[1] == LibUtilities::eGLL_Lagrange)
1991  {
1992  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1 )
1993  {
1994  IdxRange[1][0] = nummodes[1] - 2;
1995  IdxRange[1][1] = 0;
1996  Incr[1] = -1;
1997 
1998  }
1999  else
2000  {
2001  IdxRange[1][0] = 1;
2002  IdxRange[1][1] = nummodes[1] - 1;
2003  Incr[1] = 1;
2004  }
2005  }
2006  else
2007  {
2008  IdxRange[1][0] = 2;
2009  IdxRange[1][1] = nummodes[1];
2010  Incr[1] = 1;
2011 
2012  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1 )
2013  {
2014  for(i = 3; i < nummodes[1]; i+=2)
2015  {
2016  sign1[i] = -1;
2017  }
2018  }
2019  }
2020  }
2021  }
2022 
2023  switch(fid)
2024  {
2025  case 4:
2026  {
2027  IdxRange[0][0] = 0;
2028  IdxRange[0][1] = 1;
2029  Incr[0] = 1;
2030  }
2031  break;
2032  case 2:
2033  {
2034  if( bType[0] == LibUtilities::eGLL_Lagrange)
2035  {
2036  IdxRange[0][0] = nummodes[0] - 1;
2037  IdxRange[0][1] = nummodes[0];
2038  Incr[0] = 1;
2039  }
2040  else
2041  {
2042  IdxRange[0][0] = 1;
2043  IdxRange[0][1] = 2;
2044  Incr[0] = 1;
2045  }
2046  }
2047  break;
2048  default:
2049  {
2050  if( bType[0] == LibUtilities::eGLL_Lagrange)
2051  {
2052  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1 )
2053  {
2054  IdxRange[0][0] = nummodes[0] - 2;
2055  IdxRange[0][1] = 0;
2056  Incr[0] = -1;
2057 
2058  }
2059  else
2060  {
2061  IdxRange[0][0] = 1;
2062  IdxRange[0][1] = nummodes[0] - 1;
2063  Incr[0] = 1;
2064  }
2065  }
2066  else
2067  {
2068  IdxRange[0][0] = 2;
2069  IdxRange[0][1] = nummodes[0];
2070  Incr[0] = 1;
2071 
2072  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1 )
2073  {
2074  for(i = 3; i < nummodes[0]; i+=2)
2075  {
2076  sign0[i] = -1;
2077  }
2078  }
2079  }
2080  }
2081  }
2082 
2083  int p,q,r;
2084  int cnt = 0;
2085 
2086  for(r = IdxRange[2][0]; r != IdxRange[2][1]; r+=Incr[2])
2087  {
2088  for(q = IdxRange[1][0]; q != IdxRange[1][1]; q+=Incr[1])
2089  {
2090  for(p = IdxRange[0][0]; p != IdxRange[0][1]; p+=Incr[0])
2091  {
2092  maparray [ arrayindx[cnt ] ]
2093  = r*nummodes[0]*nummodes[1] + q*nummodes[0] + p;
2094  signarray[ arrayindx[cnt++] ]
2095  = sign0[p] * sign1[q] * sign2[r];
2096  }
2097  }
2098  }
2099  }
2100 
2101 
2102  /**
2103  * @param outarray Storage area for computed map.
2104  */
2105  void StdHexExp::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
2106  {
2107  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
2108  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
2109  "BasisType is not a boundary interior form");
2110  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
2111  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
2112  "BasisType is not a boundary interior form");
2113  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
2114  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
2115  "BasisType is not a boundary interior form");
2116 
2117  int i;
2118  int nummodes [3] = {m_base[0]->GetNumModes(),
2119  m_base[1]->GetNumModes(),
2120  m_base[2]->GetNumModes()};
2121 
2122  int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
2123 
2124  if(outarray.num_elements() != nIntCoeffs)
2125  {
2126  outarray = Array<OneD, unsigned int>(nIntCoeffs);
2127  }
2128 
2129  const LibUtilities::BasisType Btype [3] = {GetBasisType(0),
2130  GetBasisType(1),
2131  GetBasisType(2)};
2132 
2133  int p,q,r;
2134  int cnt = 0;
2135 
2136  int IntIdx [3][2];
2137 
2138  for(i = 0; i < 3; i++)
2139  {
2140  if( Btype[i] == LibUtilities::eModified_A)
2141  {
2142  IntIdx[i][0] = 2;
2143  IntIdx[i][1] = nummodes[i];
2144  }
2145  else
2146  {
2147  IntIdx[i][0] = 1;
2148  IntIdx[i][1] = nummodes[i]-1;
2149  }
2150  }
2151 
2152  for(r = IntIdx[2][0]; r < IntIdx[2][1]; r++)
2153  {
2154  for( q = IntIdx[1][0]; q < IntIdx[1][1]; q++)
2155  {
2156  for( p = IntIdx[0][0]; p < IntIdx[0][1]; p++)
2157  {
2158  outarray[cnt++] = r*nummodes[0]*nummodes[1] +
2159  q*nummodes[0] + p;
2160  }
2161  }
2162  }
2163  }
2164 
2165 
2166  /**
2167  * @param outarray Storage for computed map.
2168  */
2169  void StdHexExp::v_GetBoundaryMap(Array<OneD, unsigned int>& outarray)
2170  {
2171  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
2172  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
2173  "BasisType is not a boundary interior form");
2174  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
2175  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
2176  "BasisType is not a boundary interior form");
2177  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
2178  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
2179  "BasisType is not a boundary interior form");
2180 
2181  int i;
2182  int nummodes [3] = {m_base[0]->GetNumModes(),
2183  m_base[1]->GetNumModes(),
2184  m_base[2]->GetNumModes()};
2185 
2186  int nBndCoeffs = NumBndryCoeffs();
2187 
2188  if(outarray.num_elements()!=nBndCoeffs)
2189  {
2190  outarray = Array<OneD, unsigned int>(nBndCoeffs);
2191  }
2192 
2193  const LibUtilities::BasisType Btype [3] = {GetBasisType(0),
2194  GetBasisType(1),
2195  GetBasisType(2)};
2196 
2197  int p,q,r;
2198  int cnt = 0;
2199 
2200  int BndIdx [3][2];
2201  int IntIdx [3][2];
2202 
2203  for(i = 0; i < 3; i++)
2204  {
2205  BndIdx[i][0] = 0;
2206 
2207  if( Btype[i] == LibUtilities::eModified_A)
2208  {
2209  BndIdx[i][1] = 1;
2210  IntIdx[i][0] = 2;
2211  IntIdx[i][1] = nummodes[i];
2212  }
2213  else
2214  {
2215  BndIdx[i][1] = nummodes[i]-1;
2216  IntIdx[i][0] = 1;
2217  IntIdx[i][1] = nummodes[i]-1;
2218  }
2219  }
2220 
2221 
2222  for(i = 0; i < 2; i++)
2223  {
2224  r = BndIdx[2][i];
2225  for( q = 0; q < nummodes[1]; q++)
2226  {
2227  for( p = 0; p < nummodes[0]; p++)
2228  {
2229  outarray[cnt++] = r*nummodes[0]*nummodes[1]+q*nummodes[0] + p;
2230  }
2231  }
2232  }
2233 
2234  for(r = IntIdx[2][0]; r < IntIdx[2][1]; r++)
2235  {
2236  for( i = 0; i < 2; i++)
2237  {
2238  q = BndIdx[1][i];
2239  for( p = 0; p < nummodes[0]; p++)
2240  {
2241  outarray[cnt++] = r*nummodes[0]*nummodes[1] +
2242  q*nummodes[0] + p;
2243  }
2244  }
2245 
2246  for( q = IntIdx[1][0]; q < IntIdx[1][1]; q++)
2247  {
2248  for( i = 0; i < 2; i++)
2249  {
2250  p = BndIdx[0][i];
2251  outarray[cnt++] = r*nummodes[0]*nummodes[1] +
2252  q*nummodes[0] + p;
2253  }
2254  }
2255  }
2256 
2257  sort(outarray.get(), outarray.get() + nBndCoeffs);
2258  }
2259 
2260  DNekMatSharedPtr StdHexExp::v_GenMatrix(const StdMatrixKey &mkey)
2261  {
2262  return StdExpansion::CreateGeneralMatrix(mkey);
2263  }
2264 
2265 
2266  DNekMatSharedPtr StdHexExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
2267  {
2268  return StdExpansion::CreateGeneralMatrix(mkey);
2269  }
2270 
2271 
2272  void StdHexExp::v_MassMatrixOp(
2273  const Array<OneD, const NekDouble> &inarray,
2274  Array<OneD,NekDouble> &outarray,
2275  const StdMatrixKey &mkey)
2276  {
2277  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
2278  }
2279 
2280 
2281  void StdHexExp::v_LaplacianMatrixOp(
2282  const Array<OneD, const NekDouble> &inarray,
2283  Array<OneD,NekDouble> &outarray,
2284  const StdMatrixKey &mkey)
2285  {
2286  StdHexExp::v_LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
2287  }
2288 
2289 
2290  void StdHexExp::v_LaplacianMatrixOp(const int k1, const int k2,
2291  const Array<OneD, const NekDouble> &inarray,
2292  Array<OneD,NekDouble> &outarray,
2293  const StdMatrixKey &mkey)
2294  {
2295  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,
2296  mkey);
2297  }
2298 
2299 
2300  void StdHexExp::v_WeakDerivMatrixOp(const int i,
2301  const Array<OneD, const NekDouble> &inarray,
2302  Array<OneD,NekDouble> &outarray,
2303  const StdMatrixKey &mkey)
2304  {
2305  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,
2306  mkey);
2307  }
2308 
2309  void StdHexExp::v_HelmholtzMatrixOp(
2310  const Array<OneD, const NekDouble> &inarray,
2311  Array<OneD,NekDouble> &outarray,
2312  const StdMatrixKey &mkey)
2313  {
2314  StdHexExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
2315  }
2316 
2317 
2318  void StdHexExp::v_GeneralMatrixOp_MatOp(
2319  const Array<OneD, const NekDouble> &inarray,
2320  Array<OneD,NekDouble> &outarray,
2321  const StdMatrixKey &mkey)
2322  {
2323  DNekMatSharedPtr mat = m_stdMatrixManager[mkey];
2324 
2325  if(inarray.get() == outarray.get())
2326  {
2327  Array<OneD,NekDouble> tmp(m_ncoeffs);
2328  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
2329 
2330  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
2331  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
2332  }
2333  else
2334  {
2335  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
2336  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
2337  }
2338  }
2339 
2340 
2341  void StdHexExp::v_MultiplyByStdQuadratureMetric(const Array<OneD, const NekDouble>& inarray,
2342  Array<OneD, NekDouble> &outarray)
2343  {
2344  int i;
2345  int nquad0 = m_base[0]->GetNumPoints();
2346  int nquad1 = m_base[1]->GetNumPoints();
2347  int nquad2 = m_base[2]->GetNumPoints();
2348  int nq01 = nquad0*nquad1;
2349  int nq12 = nquad1*nquad2;
2350 
2351  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
2352  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
2353  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
2354 
2355  for(i = 0; i < nq12; ++i)
2356  {
2357  Vmath::Vmul(nquad0, inarray.get()+i*nquad0, 1,
2358  w0.get(), 1, outarray.get()+i*nquad0,1);
2359  }
2360 
2361  for(i = 0; i < nq12; ++i)
2362  {
2363  Vmath::Smul(nquad0, w1[i%nquad1], outarray.get()+i*nquad0, 1,
2364  outarray.get()+i*nquad0, 1);
2365  }
2366 
2367  for(i = 0; i < nquad2; ++i)
2368  {
2369  Vmath::Smul(nq01, w2[i], outarray.get()+i*nq01, 1,
2370  outarray.get()+i*nq01, 1);
2371  }
2372  }
2373 
2374  void StdHexExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
2375  const StdMatrixKey &mkey)
2376  {
2377  // Generate an orthonogal expansion
2378  int qa = m_base[0]->GetNumPoints();
2379  int qb = m_base[1]->GetNumPoints();
2380  int qc = m_base[2]->GetNumPoints();
2381  int nmodes_a = m_base[0]->GetNumModes();
2382  int nmodes_b = m_base[1]->GetNumModes();
2383  int nmodes_c = m_base[2]->GetNumModes();
2384  // Declare orthogonal basis.
2385  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
2386  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
2387  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
2388 
2389  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
2390  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A,nmodes_b,pb);
2391  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_A,nmodes_c,pc);
2392  StdHexExp OrthoExp(Ba,Bb,Bc);
2393 
2394  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2395  int i,j,k,cnt=0;
2396 
2397  // project onto modal space.
2398  OrthoExp.FwdTrans(array,orthocoeffs);
2399 
2400  if(mkey.ConstFactorExists(eFactorSVVPowerKerDiffCoeff))
2401  {
2402  // Rodrigo's power kernel
2403  NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
2404  NekDouble SvvDiffCoeff =
2405  mkey.GetConstFactor(eFactorSVVPowerKerDiffCoeff)*
2406  mkey.GetConstFactor(eFactorSVVDiffCoeff);
2407 
2408  for(int i = 0; i < nmodes_a; ++i)
2409  {
2410  for(int j = 0; j < nmodes_b; ++j)
2411  {
2412  NekDouble fac1 = std::max(
2413  pow((1.0*i)/(nmodes_a-1),cutoff*nmodes_a),
2414  pow((1.0*j)/(nmodes_b-1),cutoff*nmodes_b));
2415 
2416  for(int k = 0; k < nmodes_c; ++k)
2417  {
2418  NekDouble fac = std::max(fac1,
2419  pow((1.0*k)/(nmodes_c-1),cutoff*nmodes_c));
2420 
2421  orthocoeffs[cnt]
2422  *= SvvDiffCoeff * fac;
2423  cnt++;
2424  }
2425  }
2426  }
2427  }
2428  else if(mkey.ConstFactorExists(eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
2429  {
2430  NekDouble SvvDiffCoeff =
2431  mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff)*
2432  mkey.GetConstFactor(eFactorSVVDiffCoeff);
2433 
2434  int max_abc = max(nmodes_a-kSVVDGFiltermodesmin,
2435  nmodes_b-kSVVDGFiltermodesmin);
2436  max_abc = max(max_abc, nmodes_c-kSVVDGFiltermodesmin);
2437  // clamp max_abc
2438  max_abc = max(max_abc,0);
2439  max_abc = min(max_abc,kSVVDGFiltermodesmax-kSVVDGFiltermodesmin);
2440 
2441  for(int i = 0; i < nmodes_a; ++i)
2442  {
2443  for(int j = 0; j < nmodes_b; ++j)
2444  {
2445  int maxij = max(i,j);
2446 
2447  for(int k = 0; k < nmodes_c; ++k)
2448  {
2449  int maxijk = max(maxij,k);
2450  maxijk = min(maxijk,kSVVDGFiltermodesmax-1);
2451 
2452  orthocoeffs[cnt] *= SvvDiffCoeff *
2453  kSVVDGFilter[max_abc][maxijk];
2454  cnt++;
2455  }
2456  }
2457  }
2458  }
2459  else
2460  {
2461 
2462  int cutoff = (int) (mkey.GetConstFactor(eFactorSVVCutoffRatio)*min(nmodes_a,nmodes_b));
2463  NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
2464  // Filter just trilinear space
2465  int nmodes = max(nmodes_a,nmodes_b);
2466  nmodes = max(nmodes,nmodes_c);
2467 
2468  Array<OneD, NekDouble> fac(nmodes,1.0);
2469  for(j = cutoff; j < nmodes; ++j)
2470  {
2471  fac[j] = fabs((j-nmodes)/((NekDouble) (j-cutoff+1.0)));
2472  fac[j] *= fac[j]; //added this line to conform with equation
2473  }
2474 
2475  for(i = 0; i < nmodes_a; ++i)
2476  {
2477  for(j = 0; j < nmodes_b; ++j)
2478  {
2479  for(k = 0; k < nmodes_c; ++k)
2480  {
2481  if((i >= cutoff)||(j >= cutoff)||(k >= cutoff))
2482  {
2483  orthocoeffs[i*nmodes_a*nmodes_b +
2484  j*nmodes_c + k] *=
2485  (SvvDiffCoeff*exp(-(fac[i]+fac[j]+fac[k])));
2486  }
2487  else
2488  {
2489  orthocoeffs[i*nmodes_a*nmodes_b + j*nmodes_c + k] *= 0.0;
2490  }
2491  }
2492  }
2493  }
2494  }
2495 
2496  // backward transform to physical space
2497  OrthoExp.BwdTrans(orthocoeffs,array);
2498  }
2499 
2500  void StdHexExp::v_ExponentialFilter(
2501  Array<OneD, NekDouble> &array,
2502  const NekDouble alpha,
2503  const NekDouble exponent,
2504  const NekDouble cutoff)
2505  {
2506  // Generate an orthogonal expansion
2507  int qa = m_base[0]->GetNumPoints();
2508  int qb = m_base[1]->GetNumPoints();
2509  int qc = m_base[2]->GetNumPoints();
2510  int nmodesA = m_base[0]->GetNumModes();
2511  int nmodesB = m_base[1]->GetNumModes();
2512  int nmodesC = m_base[2]->GetNumModes();
2513  int P = nmodesA - 1;
2514  int Q = nmodesB - 1;
2515  int R = nmodesC - 1;
2516 
2517  // Declare orthogonal basis.
2518  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
2519  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
2520  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
2521 
2522  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, nmodesA, pa);
2523  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, nmodesB, pb);
2524  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_A, nmodesC, pc);
2525  StdHexExp OrthoExp(Ba,Bb,Bc);
2526 
2527  // Cutoff
2528  int Pcut = cutoff*P;
2529  int Qcut = cutoff*Q;
2530  int Rcut = cutoff*R;
2531 
2532  // Project onto orthogonal space.
2533  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2534  OrthoExp.FwdTrans(array,orthocoeffs);
2535 
2536  //
2537  NekDouble fac, fac1, fac2, fac3;
2538  int index = 0;
2539  for(int i = 0; i < nmodesA; ++i)
2540  {
2541  for(int j = 0; j < nmodesB; ++j)
2542  {
2543  for(int k = 0; k < nmodesC; ++k, ++index)
2544  {
2545  //to filter out only the "high-modes"
2546  if(i > Pcut || j > Qcut || k > Rcut)
2547  {
2548  fac1 = (NekDouble) (i - Pcut)/( (NekDouble)(P - Pcut) );
2549  fac2 = (NekDouble) (j - Qcut)/( (NekDouble)(Q - Qcut) );
2550  fac3 = (NekDouble) (k - Rcut)/( (NekDouble)(R - Rcut) );
2551  fac = max( max(fac1, fac2), fac3);
2552  fac = pow(fac, exponent);
2553  orthocoeffs[index] *= exp(-alpha*fac);
2554  }
2555  }
2556  }
2557  }
2558 
2559  // backward transform to physical space
2560  OrthoExp.BwdTrans(orthocoeffs,array);
2561  }
2562 
2563  }
2564 }
Nektar::StdRegions::StdHexExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdHexExp.cpp:2105
Nektar::StdRegions::StdHexExp::v_LocCoordToLocCollapsed
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Definition: StdHexExp.cpp:607
Nektar::StdRegions::StdHexExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2260
Nektar::StdRegions::StdHexExp::v_GetCoords
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_x, Array< OneD, NekDouble > &coords_y, Array< OneD, NekDouble > &coords_z)
Definition: StdHexExp.cpp:898
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::StdRegions::StdHexExp::v_HelmholtzMatrixOp
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2309
Nektar::StdRegions::StdHexExp::v_GetFaceNumPoints
virtual int v_GetFaceNumPoints(const int i) const
Definition: StdHexExp.cpp:801
Nektar::StdRegions::StdHexExp::v_IProductWRTDerivBase_SumFac
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:556
Nektar::LibUtilities::eModified_C
Principle Modified Functions .
Definition: BasisType.h:50
Nektar::StdRegions::kSVVDGFiltermodesmax
const int kSVVDGFiltermodesmax
Definition: StdRegions.hpp:385
Nektar::StdRegions::StdExpansion3D::BwdTrans_SumFacKernel
void BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
Definition: StdExpansion3D.cpp:114
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:787
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:469
Nektar::StdRegions::StdExpansion3D::v_LaplacianMatrixOp_MatFree
virtual void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: StdExpansion3D.cpp:207
Nektar::StdRegions::StdExpansion::IProductWRTDerivBase_SumFac
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:1474
Nektar::StdRegions::StdExpansion::GetNumPoints
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:228
Nektar::StdRegions::StdExpansion::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:945
Nektar::StdRegions::StdExpansion::GetBasisNumModes
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
Definition: StdExpansion.h:177
Nektar::StdRegions::StdHexExp::v_GetFaceIntNcoeffs
virtual int v_GetFaceIntNcoeffs(const int i) const
Definition: StdHexExp.cpp:776
Nektar::StdRegions::StdExpansion::IProductWRTBase
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
this function calculates the inner product of a given function f with the different modes of the expa...
Definition: StdExpansion.h:634
Nektar::StdRegions::StdHexExp::v_IProductWRTBase
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:360
Nektar::StdRegions::StdExpansion::LaplacianMatrixOp_MatFree
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:1508
Nektar::StdRegions::StdHexExp::v_FillMode
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:618
Nektar::StdRegions::StdHexExp::v_GetBoundaryMap
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdHexExp.cpp:2169
Nektar::LibUtilities::eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
std
STL namespace.
Nektar::StdRegions::StdHexExp::v_CalcNumberOfCoefficients
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdHexExp.cpp:842
Nektar::StdRegions::StdHexExp::v_DetFaceBasisKey
virtual const LibUtilities::BasisKey v_DetFaceBasisKey(const int i, const int k) const
Definition: StdHexExp.cpp:851
Nektar::StdRegions::StdHexExp::v_FwdTrans
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:299
Nektar::StdRegions::StdHexExp::v_GetFaceNumModes
virtual void v_GetFaceNumModes(const int fid, const Orientation faceOrient, int &numModes0, int &numModes1)
Definition: StdHexExp.cpp:924
Nektar::StdRegions::kSVVDGFiltermodesmin
const int kSVVDGFiltermodesmin
Definition: StdRegions.hpp:384
Nektar::StdRegions::StdHexExp::v_MassMatrixOp
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2272
Nektar::StdRegions::StdHexExp::v_WeakDerivMatrixOp
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2300
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
Blas::Dgemm
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where A[m x n], B[n x k], C[m x k].
Definition: Blas.hpp:213
Nektar::StdRegions::StdHexExp::v_IsBoundaryInteriorExpansion
virtual bool v_IsBoundaryInteriorExpansion()
Definition: StdHexExp.cpp:75
Nektar::StdRegions::StdHexExp::v_IProductWRTBase_SumFacKernel
virtual void v_IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
Definition: StdHexExp.cpp:431
class_topology.P
P
Definition: class_topology.py:290
Nektar::StdRegions::StdHexExp::v_GeneralMatrixOp_MatOp
virtual void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2318
Nektar::StdRegions::StdExpansion::PhysDeriv
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Definition: StdExpansion.h:1065
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:714
Nektar::LibUtilities::eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:46
Nektar::StdRegions::StdHexExp::v_GetFacePointsKey
virtual LibUtilities::PointsKey v_GetFacePointsKey(const int i, const int j) const
Definition: StdHexExp.cpp:822
Nektar::StdRegions::StdHexExp::v_IProductWRTDerivBase_MatOp
virtual void v_IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:526
Nektar::StdRegions::StdExpansion3D::PhysTensorDeriv
void PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d1, Array< OneD, NekDouble > &outarray_d2, Array< OneD, NekDouble > &outarray_d3)
Calculate the 3D derivative in the local tensor/collapsed coordinate at the physical points...
Definition: StdExpansion3D.cpp:69
Nektar::StdRegions::StdHexExp::v_MultiplyByStdQuadratureMetric
virtual void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:2341
Nektar::StdRegions::StdHexExp::v_GetEdgeInteriorMap
virtual void v_GetEdgeInteriorMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
Definition: StdHexExp.cpp:1467
Nektar::StdRegions::StdExpansion
The base class for all shapes.
Definition: StdExpansion.h:68
Nektar::StdRegions::StdExpansion3D::IProductWRTBase_SumFacKernel
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
Definition: StdExpansion3D.cpp:128
Nektar::StdRegions::StdExpansion::GetEdgeNcoeffs
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
Definition: StdExpansion.h:286
Nektar::StdRegions::StdHexExp::v_NumBndryCoeffs
virtual int v_NumBndryCoeffs() const
Definition: StdHexExp.cpp:693
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
Nektar::StdRegions::StdHexExp::v_GetEdgeNcoeffs
virtual int v_GetEdgeNcoeffs(const int i) const
Definition: StdHexExp.cpp:734
StdHexExp.h
Nektar::StdRegions::EvaluateQuadFaceBasisKey
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Definition: StdExpansion3D.cpp:351
Nektar::StdRegions::StdExpansion::GetFaceIntNcoeffs
int GetFaceIntNcoeffs(const int i) const
Definition: StdExpansion.h:358
Nektar::StdRegions::StdHexExp::v_LaplacianMatrixOp
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2281
Nektar::LibUtilities::eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::StdRegions::StdHexExp
Class representing a hexehedral element in reference space.
Definition: StdHexExp.h:47
Nektar::StdRegions::StdHexExp::v_GetEdgeBasisType
virtual LibUtilities::BasisType v_GetEdgeBasisType(const int i) const
Definition: StdHexExp.cpp:880
Nektar::StdRegions::StdHexExp::v_ExponentialFilter
virtual void v_ExponentialFilter(Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff)
Definition: StdHexExp.cpp:2500
Nektar::StdRegions::StdHexExp::v_GetFaceNcoeffs
virtual int v_GetFaceNcoeffs(const int i) const
Definition: StdHexExp.cpp:758
Nektar::StdRegions::StdHexExp::v_DetShapeType
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdHexExp.cpp:687
Nektar::StdRegions::StdExpansion::WeakDerivMatrixOp_MatFree
void WeakDerivMatrixOp_MatFree(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.cpp:786
Nektar::LibUtilities::eOrtho_C
Principle Orthogonal Functions .
Definition: BasisType.h:47
Nektar::StdRegions::StdHexExp::v_IProductWRTBase_MatOp
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:379
Nektar::StdRegions::StdExpansion::CreateGeneralMatrix
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
Definition: StdExpansion.cpp:324
Nektar::LibUtilities::eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:45
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:59
Nektar::StdRegions::StdHexExp::v_BwdTrans_SumFacKernel
virtual void v_BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
Definition: StdHexExp.cpp:235
Nektar::StdRegions::StdHexExp::v_GetTotalFaceIntNcoeffs
virtual int v_GetTotalFaceIntNcoeffs() const
Definition: StdHexExp.cpp:794
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::StdRegions::StdHexExp::v_BwdTrans_SumFac
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:208
Blas::Dgemv
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = A x where A[m x n].
Definition: Blas.hpp:168
Nektar::StdRegions::StdHexExp::v_GetVertexMap
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdHexExp.cpp:1333
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:121
Nektar::StdRegions::StdHexExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2374
Nektar::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:130
Nektar::StdRegions::kSVVDGFilter
const NekDouble kSVVDGFilter[9][11]
Definition: StdRegions.hpp:387
Nektar::StdRegions::StdHexExp::v_GetFaceToElementMap
virtual void v_GetFaceToElementMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int nummodesA=-1, int nummodesB=-1)
Definition: StdHexExp.cpp:972
Blas::Dscal
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:125
Nektar::StdRegions::StdExpansion::GetPointsType
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:215
Nektar::StdRegions::StdHexExp::v_BwdTrans
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:178
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:130
Nektar::StdRegions::StdHexExp::v_IProductWRTDerivBase
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:518
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:274
Nektar::StdRegions::StdExpansion::GetBasisType
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:164
Nektar::StdRegions::StdExpansion::MassMatrixOp_MatFree
void MassMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.cpp:669
Nektar::StdRegions::StdHexExp::v_GetNedges
virtual int v_GetNedges() const
Definition: StdHexExp.cpp:675
Nektar::StdRegions::StdHexExp::v_GetTotalEdgeIntNcoeffs
virtual int v_GetTotalEdgeIntNcoeffs() const
Definition: StdHexExp.cpp:752
Nektar::StdRegions::StdExpansion::BwdTrans
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space...
Definition: StdExpansion.h:530
Nektar::StdRegions::StdHexExp::v_NumDGBndryCoeffs
virtual int v_NumDGBndryCoeffs() const
Definition: StdHexExp.cpp:714
Nektar::StdRegions::StdHexExp::v_GetNfaces
virtual int v_GetNfaces() const
Definition: StdHexExp.cpp:681
Nektar::StdRegions::StdExpansion::m_stdMatrixManager
LibUtilities::NekManager< StdMatrixKey, DNekMat, StdMatrixKey::opLess > m_stdMatrixManager
Definition: StdExpansion.h:1432
Nektar::StdRegions::StdHexExp::v_PhysDeriv
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Differentiation Methods.
Definition: StdHexExp.cpp:94
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:140
Nektar::LibUtilities::eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:54
Nektar::StdRegions::StdHexExp::v_GetFaceInteriorMap
virtual void v_GetFaceInteriorMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
Definition: StdHexExp.cpp:1765
Nektar::StdRegions::StdExpansion3D::v_HelmholtzMatrixOp_MatFree
virtual void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: StdExpansion3D.cpp:249
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250
Nektar::StdRegions::StdHexExp::v_GetNverts
virtual int v_GetNverts() const
Definition: StdHexExp.cpp:669
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1429
Nektar::StdRegions::StdHexExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2266
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Nektar::StdRegions::StdHexExp::v_StdPhysDeriv
virtual void v_StdPhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Definition: StdHexExp.cpp:141
Nektar::StdRegions::StdExpansion::GetEdgeBasisType
LibUtilities::BasisType GetEdgeBasisType(const int i) const
This function returns the type of expansion basis on the i-th edge.
Definition: StdExpansion.h:412
Nektar::StdRegions::StdHexExp::v_IProductWRTBase_SumFac
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multbyweights=true)
Definition: StdHexExp.cpp:393
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:49
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186
Nektar::StdRegions::StdExpansion::BwdTrans_SumFac
void BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:1467
Nektar::StdRegions::StdExpansion::FwdTrans
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space...
Definition: StdExpansion.h:1957
Generated on Wed Jan 20 2021 03:33:41 for Nektar++ by   doxygen 1.8.13