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(mode == mode2 * btmp0 * btmp1 + mode1 * btmp0 + mode0,
637  "Mode lookup failed.");
638  ASSERTL2(mode < m_ncoeffs,
639  "Calling argument mode is larger than total expansion "
640  "order");
641 
642  for(i = 0; i < nquad1*nquad2; ++i)
643  {
644  Vmath::Vcopy(nquad0,(NekDouble *)(base0.get() + mode0*nquad0),1,
645  &outarray[0]+i*nquad0, 1);
646  }
647 
648  for(j = 0; j < nquad2; ++j)
649  {
650  for(i = 0; i < nquad0; ++i)
651  {
652  Vmath::Vmul(nquad1,(NekDouble *)(base1.get() + mode1*nquad1),1,
653  &outarray[0]+i+j*nquad0*nquad1, nquad0,
654  &outarray[0]+i+j*nquad0*nquad1, nquad0);
655  }
656  }
657 
658  for(i = 0; i < nquad2; i++)
659  {
660  Blas::Dscal(nquad0*nquad1,base2[mode2*nquad2+i],
661  &outarray[0]+i*nquad0*nquad1,1);
662  }
663  }
664 
665 
666  int StdHexExp::v_GetNverts() const
667  {
668  return 8;
669  }
670 
671 
672  int StdHexExp::v_GetNedges() const
673  {
674  return 12;
675  }
676 
677 
678  int StdHexExp::v_GetNfaces() const
679  {
680  return 6;
681  }
682 
683 
684  LibUtilities::ShapeType StdHexExp::v_DetShapeType() const
685  {
686  return LibUtilities::eHexahedron;
687  }
688 
689 
690  int StdHexExp::v_NumBndryCoeffs() const
691  {
692  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
693  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
694  "BasisType is not a boundary interior form");
695  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
696  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
697  "BasisType is not a boundary interior form");
698  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
699  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
700  "BasisType is not a boundary interior form");
701 
702  int nmodes0 = m_base[0]->GetNumModes();
703  int nmodes1 = m_base[1]->GetNumModes();
704  int nmodes2 = m_base[2]->GetNumModes();
705 
706  return ( 2*( nmodes0*nmodes1 + nmodes0*nmodes2
707  + nmodes1*nmodes2)
708  - 4*( nmodes0 + nmodes1 + nmodes2 ) + 8 );
709  }
710 
711  int StdHexExp::v_NumDGBndryCoeffs() const
712  {
713  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
714  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
715  "BasisType is not a boundary interior form");
716  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
717  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
718  "BasisType is not a boundary interior form");
719  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
720  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
721  "BasisType is not a boundary interior form");
722 
723  int nmodes0 = m_base[0]->GetNumModes();
724  int nmodes1 = m_base[1]->GetNumModes();
725  int nmodes2 = m_base[2]->GetNumModes();
726 
727  return 2*( nmodes0*nmodes1 + nmodes0*nmodes2
728  + nmodes1*nmodes2 );
729  }
730 
731  int StdHexExp::v_GetEdgeNcoeffs(const int i) const
732  {
733  ASSERTL2((i >= 0)&&(i <= 11),"edge id is out of range");
734 
735  if((i == 0)||(i == 2)||(i == 8)||(i == 10))
736  {
737  return GetBasisNumModes(0);
738  }
739  else if((i == 1)||(i == 3)||(i == 9)||(i == 11))
740  {
741  return GetBasisNumModes(1);
742  }
743  else
744  {
745  return GetBasisNumModes(2);
746  }
747  }
748 
749  int StdHexExp::v_GetTotalEdgeIntNcoeffs() const
750  {
751  return 4*(GetBasisNumModes(0)+GetBasisNumModes(1)+GetBasisNumModes(2));
752  }
753 
754 
755  int StdHexExp::v_GetFaceNcoeffs(const int i) const
756  {
757  ASSERTL2((i >= 0) && (i <= 5), "face id is out of range");
758  if((i == 0) || (i == 5))
759  {
760  return GetBasisNumModes(0)*GetBasisNumModes(1);
761  }
762  else if((i == 1) || (i == 3))
763  {
764  return GetBasisNumModes(0)*GetBasisNumModes(2);
765  }
766  else
767  {
768  return GetBasisNumModes(1)*GetBasisNumModes(2);
769  }
770  }
771 
772 
773  int StdHexExp::v_GetFaceIntNcoeffs(const int i) const
774  {
775  ASSERTL2((i >= 0) && (i <= 5), "face id is out of range");
776  if((i == 0) || (i == 5))
777  {
778  return (GetBasisNumModes(0)-2)*(GetBasisNumModes(1)-2);
779  }
780  else if((i == 1) || (i == 3))
781  {
782  return (GetBasisNumModes(0)-2)*(GetBasisNumModes(2)-2);
783  }
784  else
785  {
786  return (GetBasisNumModes(1)-2)*(GetBasisNumModes(2)-2);
787  }
788 
789  }
790 
791  int StdHexExp::v_GetTotalFaceIntNcoeffs() const
792  {
793  return 2*((GetBasisNumModes(0)-2)*(GetBasisNumModes(1)-2)+
794  (GetBasisNumModes(0)-2)*(GetBasisNumModes(2)-2)+
795  (GetBasisNumModes(1)-2)*(GetBasisNumModes(2)-2));
796  }
797 
798  int StdHexExp::v_GetFaceNumPoints(const int i) const
799  {
800  ASSERTL2(i >= 0 && i <= 5, "face id is out of range");
801 
802  if (i == 0 || i == 5)
803  {
804  return m_base[0]->GetNumPoints()*
805  m_base[1]->GetNumPoints();
806  }
807  else if (i == 1 || i == 3)
808  {
809  return m_base[0]->GetNumPoints()*
810  m_base[2]->GetNumPoints();
811  }
812  else
813  {
814  return m_base[1]->GetNumPoints()*
815  m_base[2]->GetNumPoints();
816  }
817  }
818 
819  LibUtilities::PointsKey StdHexExp::v_GetFacePointsKey(
820  const int i, const int j) const
821  {
822  ASSERTL2(i >= 0 && i <= 5, "face id is out of range");
823  ASSERTL2(j == 0 || j == 1, "face direction is out of range");
824 
825  if (i == 0 || i == 5)
826  {
827  return m_base[j]->GetPointsKey();
828  }
829  else if (i == 1 || i == 3)
830  {
831  return m_base[2*j]->GetPointsKey();
832  }
833  else
834  {
835  return m_base[j+1]->GetPointsKey();
836  }
837  }
838 
839  int StdHexExp::v_CalcNumberOfCoefficients(const std::vector<unsigned int> &nummodes, int &modes_offset)
840  {
841  int nmodes = nummodes[modes_offset]*nummodes[modes_offset+1]*nummodes[modes_offset+2];
842  modes_offset += 3;
843 
844  return nmodes;
845  }
846 
847 
848  const LibUtilities::BasisKey StdHexExp::v_DetFaceBasisKey(
849  const int i, const int k) const
850  {
851  ASSERTL2(i >= 0 && i <= 5, "face id is out of range");
852  ASSERTL2(k >= 0 && k <= 1, "basis key id is out of range");
853 
854  int dir = k;
855  switch(i)
856  {
857  case 0:
858  case 5:
859  dir = k;
860  break;
861  case 1:
862  case 3:
863  dir = 2*k;
864  break;
865  case 2:
866  case 4:
867  dir = k+1;
868  break;
869  }
870 
871  return EvaluateQuadFaceBasisKey(k,
872  m_base[dir]->GetBasisType(),
873  m_base[dir]->GetNumPoints(),
874  m_base[dir]->GetNumModes());
875  }
876 
877  LibUtilities::BasisType StdHexExp::v_GetEdgeBasisType(const int i) const
878  {
879  ASSERTL2((i >= 0)&&(i <= 11),"edge id is out of range");
880 
881  if((i == 0)||(i == 2)||(i==8)||(i==10))
882  {
883  return GetBasisType(0);
884  }
885  else if((i == 1)||(i == 3)||(i == 9)||(i == 11))
886  {
887  return GetBasisType(1);
888  }
889  else
890  {
891  return GetBasisType(2);
892  }
893  }
894 
895  void StdHexExp::v_GetCoords( Array<OneD, NekDouble> & xi_x,
896  Array<OneD, NekDouble> & xi_y,
897  Array<OneD, NekDouble> & xi_z)
898  {
899  Array<OneD, const NekDouble> eta_x = m_base[0]->GetZ();
900  Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
901  Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
902  int Qx = GetNumPoints(0);
903  int Qy = GetNumPoints(1);
904  int Qz = GetNumPoints(2);
905 
906  // Convert collapsed coordinates into cartesian coordinates:
907  // eta --> xi
908  for( int k = 0; k < Qz; ++k ) {
909  for( int j = 0; j < Qy; ++j ) {
910  for( int i = 0; i < Qx; ++i ) {
911  int s = i + Qx*(j + Qy*k);
912  xi_x[s] = eta_x[i];
913  xi_y[s] = eta_y[j];
914  xi_z[s] = eta_z[k];
915 
916  }
917  }
918  }
919  }
920 
921  void StdHexExp::v_GetFaceNumModes(
922  const int fid,
923  const Orientation faceOrient,
924  int &numModes0,
925  int &numModes1)
926  {
927  int nummodes [3] = {m_base[0]->GetNumModes(),
928  m_base[1]->GetNumModes(),
929  m_base[2]->GetNumModes()};
930  switch(fid)
931  {
932  case 0:
933  case 5:
934  {
935  numModes0 = nummodes[0];
936  numModes1 = nummodes[1];
937  }
938  break;
939  case 1:
940  case 3:
941  {
942  numModes0 = nummodes[0];
943  numModes1 = nummodes[2];
944  }
945  break;
946  case 2:
947  case 4:
948  {
949  numModes0 = nummodes[1];
950  numModes1 = nummodes[2];
951  }
952  break;
953  default:
954  {
955  ASSERTL0(false,"fid out of range");
956  }
957  break;
958  }
959 
960  if ( faceOrient >= eDir1FwdDir2_Dir2FwdDir1 )
961  {
962  std::swap(numModes0, numModes1);
963  }
964  }
965 
966  /**
967  * Only for basis type Modified_A or GLL_LAGRANGE in all directions.
968  */
969  void StdHexExp::v_GetFaceToElementMap(
970  const int fid,
971  const Orientation faceOrient,
972  Array<OneD, unsigned int> &maparray,
973  Array<OneD, int> &signarray,
974  int P,
975  int Q)
976  {
977  int i,j;
978  int nummodesA=0, nummodesB=0;
979 
980  ASSERTL1(GetEdgeBasisType(0) == GetEdgeBasisType(1) &&
981  GetEdgeBasisType(0) == GetEdgeBasisType(2),
982  "Method only implemented if BasisType is indentical in "
983  "all directions");
984  ASSERTL1(GetEdgeBasisType(0) == LibUtilities::eModified_A ||
985  GetEdgeBasisType(0) == LibUtilities::eGLL_Lagrange,
986  "Method only implemented for Modified_A or GLL_Lagrange BasisType");
987 
988  const int nummodes0 = m_base[0]->GetNumModes();
989  const int nummodes1 = m_base[1]->GetNumModes();
990  const int nummodes2 = m_base[2]->GetNumModes();
991 
992  switch(fid)
993  {
994  case 0:
995  case 5:
996  nummodesA = nummodes0;
997  nummodesB = nummodes1;
998  break;
999  case 1:
1000  case 3:
1001  nummodesA = nummodes0;
1002  nummodesB = nummodes2;
1003  break;
1004  case 2:
1005  case 4:
1006  nummodesA = nummodes1;
1007  nummodesB = nummodes2;
1008  break;
1009  default:
1010  ASSERTL0(false,"fid must be between 0 and 5");
1011  }
1012 
1013  bool CheckForZeroedModes = false;
1014 
1015  if (P == -1)
1016  {
1017  P = nummodesA;
1018  Q = nummodesB;
1019  }
1020 
1021  if((P != nummodesA)||(Q != nummodesB))
1022  {
1023  CheckForZeroedModes = true;
1024  }
1025 
1026  bool modified = (GetEdgeBasisType(0) == LibUtilities::eModified_A);
1027  int nFaceCoeffs = P*Q;
1028 
1029  if(maparray.num_elements() != nFaceCoeffs)
1030  {
1031  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1032  }
1033 
1034  if(signarray.num_elements() != nFaceCoeffs)
1035  {
1036  signarray = Array<OneD, int>(nFaceCoeffs,1);
1037  }
1038  else
1039  {
1040  fill( signarray.get() , signarray.get()+nFaceCoeffs, 1 );
1041  }
1042 
1043  Array<OneD, int> arrayindx(nFaceCoeffs);
1044 
1045  for(i = 0; i < Q; i++)
1046  {
1047  for(j = 0; j < P; j++)
1048  {
1049  if( faceOrient < eDir1FwdDir2_Dir2FwdDir1 )
1050  {
1051  arrayindx[i*P+j] = i*P+j;
1052  }
1053  else
1054  {
1055  arrayindx[i*P+j] = j*Q+i;
1056  }
1057  }
1058  }
1059 
1060  int offset = 0;
1061  int jump1 = 1;
1062  int jump2 = 1;
1063 
1064  switch(fid)
1065  {
1066  case 5:
1067  {
1068  if (modified)
1069  {
1070  offset = nummodes0*nummodes1;
1071  }
1072  else
1073  {
1074  offset = (nummodes2-1)*nummodes0*nummodes1;
1075  jump1 = nummodes0;
1076  }
1077  }
1078  /* Falls through. */
1079  case 0:
1080  {
1081  jump1 = nummodes0;
1082  break;
1083  }
1084  case 3:
1085  {
1086  if (modified)
1087  {
1088  offset = nummodes0;
1089  }
1090  else
1091  {
1092  offset = nummodes0*(nummodes1-1);
1093  jump1 = nummodes0*nummodes1;
1094  }
1095  }
1096  /* Falls through. */
1097  case 1:
1098  {
1099  jump1 = nummodes0*nummodes1;
1100  break;
1101  }
1102  case 2:
1103  {
1104  if (modified)
1105  {
1106  offset = 1;
1107  }
1108  else
1109  {
1110  offset = nummodes0-1;
1111  jump1 = nummodes0*nummodes1;
1112  jump2 = nummodes0;
1113 
1114  }
1115  }
1116  /* Falls through. */
1117  case 4:
1118  {
1119  jump1 = nummodes0*nummodes1;
1120  jump2 = nummodes0;
1121  break;
1122  }
1123  default:
1124  ASSERTL0(false,"fid must be between 0 and 5");
1125  }
1126 
1127  for(i = 0; i < Q; i++)
1128  {
1129  for(j = 0; j < P; j++)
1130  {
1131  maparray[ arrayindx[i*P+j] ]
1132  = i*jump1 + j*jump2 + offset;
1133  }
1134  }
1135 
1136 
1137  if(CheckForZeroedModes)
1138  {
1139  if(modified)
1140  {
1141  // zero signmap and set maparray to zero if elemental
1142  // modes are not as large as face modesl
1143  for(i = 0; i < nummodesB; i++)
1144  {
1145  for(j = nummodesA; j < P; j++)
1146  {
1147  signarray[arrayindx[i*P+j]] = 0.0;
1148  maparray[arrayindx[i*P+j]] = maparray[0];
1149  }
1150  }
1151 
1152  for(i = nummodesB; i < Q; i++)
1153  {
1154  for(j = 0; j < P; j++)
1155  {
1156  signarray[arrayindx[i*P+j]] = 0.0;
1157  maparray[arrayindx[i*P+j]] = maparray[0];
1158  }
1159  }
1160  }
1161  else
1162  {
1163  ASSERTL0(false,"Different trace space face dimention and element face dimention not possible for GLL-Lagrange bases");
1164  }
1165  }
1166 
1167  if( (faceOrient==eDir1FwdDir1_Dir2BwdDir2) ||
1168  (faceOrient==eDir1BwdDir1_Dir2BwdDir2) ||
1169  (faceOrient==eDir1BwdDir2_Dir2FwdDir1) ||
1170  (faceOrient==eDir1BwdDir2_Dir2BwdDir1) )
1171  {
1172  if(faceOrient<eDir1FwdDir2_Dir2FwdDir1)
1173  {
1174  if (modified)
1175  {
1176  for(i = 3; i < Q; i+=2)
1177  {
1178  for(j = 0; j < P; j++)
1179  {
1180  signarray[ arrayindx[i*P+j] ] *= -1;
1181  }
1182  }
1183 
1184  for(i = 0; i < P; i++)
1185  {
1186  swap( maparray[i] , maparray[i+P] );
1187  swap( signarray[i] , signarray[i+P] );
1188  }
1189 
1190  }
1191  else
1192  {
1193  for(i = 0; i < P; i++)
1194  {
1195  for(j = 0; j < Q/2; j++)
1196  {
1197  swap( maparray[i + j*P],
1198  maparray[i+P*Q
1199  -P -j*P] );
1200  swap( signarray[i + j*P],
1201  signarray[i+P*Q
1202  -P -j*P]);
1203  }
1204  }
1205  }
1206  }
1207  else
1208  {
1209  if (modified)
1210  {
1211  for(i = 0; i < Q; i++)
1212  {
1213  for(j = 3; j < P; j+=2)
1214  {
1215  signarray[ arrayindx[i*P+j] ] *= -1;
1216  }
1217  }
1218 
1219  for(i = 0; i < Q; i++)
1220  {
1221  swap( maparray[i] , maparray[i+Q] );
1222  swap( signarray[i] , signarray[i+Q] );
1223  }
1224 
1225  }
1226  else
1227  {
1228  for(i = 0; i < P; i++)
1229  {
1230  for(j = 0; j < Q/2; j++)
1231  {
1232  swap( maparray[i*Q + j],
1233  maparray[i*Q + Q -1 -j]);
1234  swap( signarray[i*Q + j],
1235  signarray[i*Q + Q -1 -j]);
1236  }
1237  }
1238  }
1239  }
1240  }
1241 
1242  if( (faceOrient==eDir1BwdDir1_Dir2FwdDir2) ||
1243  (faceOrient==eDir1BwdDir1_Dir2BwdDir2) ||
1244  (faceOrient==eDir1FwdDir2_Dir2BwdDir1) ||
1245  (faceOrient==eDir1BwdDir2_Dir2BwdDir1) )
1246  {
1247  if(faceOrient<eDir1FwdDir2_Dir2FwdDir1)
1248  {
1249  if (modified)
1250  {
1251  for(i = 0; i < Q; i++)
1252  {
1253  for(j = 3; j < P; j+=2)
1254  {
1255  signarray[ arrayindx[i*P+j] ] *= -1;
1256  }
1257  }
1258 
1259  for(i = 0; i < Q; i++)
1260  {
1261  swap( maparray[i*P],
1262  maparray[i*P+1]);
1263  swap( signarray[i*P],
1264  signarray[i*P+1]);
1265  }
1266  }
1267  else
1268  {
1269  for(i = 0; i < Q; i++)
1270  {
1271  for(j = 0; j < P/2; j++)
1272  {
1273  swap( maparray[i*P + j],
1274  maparray[i*P + P -1 -j]);
1275  swap( signarray[i*P + j],
1276  signarray[i*P + P -1 -j]);
1277  }
1278  }
1279  }
1280 
1281 
1282 
1283  }
1284  else
1285  {
1286  if (modified)
1287  {
1288  for(i = 3; i < Q; i+=2)
1289  {
1290  for(j = 0; j < P; j++)
1291  {
1292  signarray[ arrayindx[i*P+j] ] *= -1;
1293  }
1294  }
1295 
1296  for(i = 0; i < P; i++)
1297  {
1298  swap( maparray[i*Q],
1299  maparray[i*Q+1]);
1300  swap( signarray[i*Q],
1301  signarray[i*Q+1]);
1302  }
1303  }
1304  else
1305  {
1306  for(i = 0; i < Q; i++)
1307  {
1308  for(j = 0; j < P/2; j++)
1309  {
1310  swap( maparray[i + j*Q] ,
1311  maparray[i+P*Q - Q -j*Q] );
1312  swap( signarray[i + j*Q] ,
1313  signarray[i+P*Q - Q -j*Q] );
1314  }
1315  }
1316  }
1317  }
1318  }
1319  }
1320 
1321 
1322 
1323  /**
1324  * Expansions in each of the three dimensions must be of type
1325  * LibUtilities#eModified_A or LibUtilities#eGLL_Lagrange.
1326  *
1327  * @param localVertexId ID of vertex (0..7)
1328  * @returns Position of vertex in local numbering scheme.
1329  */
1330  int StdHexExp::v_GetVertexMap(const int localVertexId, bool useCoeffPacking)
1331  {
1332  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1333  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1334  "BasisType is not a boundary interior form");
1335  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1336  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1337  "BasisType is not a boundary interior form");
1338  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
1339  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1340  "BasisType is not a boundary interior form");
1341 
1342  ASSERTL1((localVertexId>=0)&&(localVertexId<8),
1343  "local vertex id must be between 0 and 7");
1344 
1345  int p = 0;
1346  int q = 0;
1347  int r = 0;
1348 
1349  // Retrieve the number of modes in each dimension.
1350  int nummodes [3] = {m_base[0]->GetNumModes(),
1351  m_base[1]->GetNumModes(),
1352  m_base[2]->GetNumModes()};
1353 
1354  if(useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1355  {
1356  if(localVertexId > 3)
1357  {
1358  if( GetBasisType(2) == LibUtilities::eGLL_Lagrange)
1359  {
1360  r = nummodes[2]-1;
1361  }
1362  else
1363  {
1364  r = 1;
1365  }
1366  }
1367 
1368  switch(localVertexId % 4)
1369  {
1370  case 0:
1371  break;
1372  case 1:
1373  {
1374  if( GetBasisType(0) == LibUtilities::eGLL_Lagrange)
1375  {
1376  p = nummodes[0]-1;
1377  }
1378  else
1379  {
1380  p = 1;
1381  }
1382  }
1383  break;
1384  case 2:
1385  {
1386  if( GetBasisType(1) == LibUtilities::eGLL_Lagrange)
1387  {
1388  q = nummodes[1]-1;
1389  }
1390  else
1391  {
1392  q = 1;
1393  }
1394  }
1395  break;
1396  case 3:
1397  {
1398  if( GetBasisType(1) == LibUtilities::eGLL_Lagrange)
1399  {
1400  p = nummodes[0]-1;
1401  q = nummodes[1]-1;
1402  }
1403  else
1404  {
1405  p = 1;
1406  q = 1;
1407  }
1408  }
1409  break;
1410  }
1411  }
1412  else
1413  {
1414  // Right face (vertices 1,2,5,6)
1415  if( (localVertexId % 4) % 3 > 0 )
1416  {
1417  if( GetBasisType(0) == LibUtilities::eGLL_Lagrange)
1418  {
1419  p = nummodes[0]-1;
1420  }
1421  else
1422  {
1423  p = 1;
1424  }
1425  }
1426 
1427  // Back face (vertices 2,3,6,7)
1428  if( localVertexId % 4 > 1 )
1429  {
1430  if( GetBasisType(1) == LibUtilities::eGLL_Lagrange)
1431  {
1432  q = nummodes[1]-1;
1433  }
1434  else
1435  {
1436  q = 1;
1437  }
1438  }
1439 
1440  // Top face (vertices 4,5,6,7)
1441  if( localVertexId > 3)
1442  {
1443  if( GetBasisType(2) == LibUtilities::eGLL_Lagrange)
1444  {
1445  r = nummodes[2]-1;
1446  }
1447  else
1448  {
1449  r = 1;
1450  }
1451  }
1452  }
1453  // Compute the local number.
1454  return r*nummodes[0]*nummodes[1] + q*nummodes[0] + p;
1455  }
1456 
1457 
1458  /**
1459  * @param eid The edge to compute the numbering for.
1460  * @param edgeOrient Orientation of the edge.
1461  * @param maparray Storage for computed mapping array.
1462  * @param signarray ?
1463  */
1464  void StdHexExp::v_GetEdgeInteriorMap(const int eid,
1465  const Orientation edgeOrient,
1466  Array<OneD, unsigned int> &maparray,
1467  Array<OneD, int> &signarray)
1468  {
1469  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1470  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1471  "BasisType is not a boundary interior form");
1472  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1473  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1474  "BasisType is not a boundary interior form");
1475  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
1476  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1477  "BasisType is not a boundary interior form");
1478 
1479  ASSERTL1((eid>=0)&&(eid<12),
1480  "local edge id must be between 0 and 11");
1481 
1482  int nEdgeIntCoeffs = GetEdgeNcoeffs(eid)-2;
1483 
1484  if(maparray.num_elements()!=nEdgeIntCoeffs)
1485  {
1486  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1487  }
1488 
1489  if(signarray.num_elements() != nEdgeIntCoeffs)
1490  {
1491  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1492  }
1493  else
1494  {
1495  fill( signarray.get() , signarray.get()+nEdgeIntCoeffs, 1 );
1496  }
1497 
1498  int nummodes [3] = {m_base[0]->GetNumModes(),
1499  m_base[1]->GetNumModes(),
1500  m_base[2]->GetNumModes()};
1501 
1502  const LibUtilities::BasisType bType [3] = {GetBasisType(0),
1503  GetBasisType(1),
1504  GetBasisType(2)};
1505 
1506  bool reverseOrdering = false;
1507  bool signChange = false;
1508 
1509  int IdxRange [3][2] = {{0,0},{0,0},{0,0}};
1510 
1511  switch(eid)
1512  {
1513  case 0:
1514  case 1:
1515  case 2:
1516  case 3:
1517  {
1518  IdxRange[2][0] = 0;
1519  IdxRange[2][1] = 1;
1520  }
1521  break;
1522  case 8:
1523  case 9:
1524  case 10:
1525  case 11:
1526  {
1527  if( bType[2] == LibUtilities::eGLL_Lagrange)
1528  {
1529  IdxRange[2][0] = nummodes[2] - 1;
1530  IdxRange[2][1] = nummodes[2];
1531  }
1532  else
1533  {
1534  IdxRange[2][0] = 1;
1535  IdxRange[2][1] = 2;
1536  }
1537  }
1538  break;
1539  case 4:
1540  case 5:
1541  case 6:
1542  case 7:
1543  {
1544  if( bType[2] == LibUtilities::eGLL_Lagrange)
1545  {
1546  IdxRange[2][0] = 1;
1547  IdxRange[2][1] = nummodes[2] - 1;
1548 
1549  if(edgeOrient==eBackwards)
1550  {
1551  reverseOrdering = true;
1552  }
1553  }
1554  else
1555  {
1556  IdxRange[2][0] = 2;
1557  IdxRange[2][1] = nummodes[2];
1558 
1559  if(edgeOrient==eBackwards)
1560  {
1561  signChange = true;
1562  }
1563  }
1564  }
1565  break;
1566  }
1567 
1568  switch(eid)
1569  {
1570  case 0:
1571  case 4:
1572  case 5:
1573  case 8:
1574  {
1575  IdxRange[1][0] = 0;
1576  IdxRange[1][1] = 1;
1577  }
1578  break;
1579  case 2:
1580  case 6:
1581  case 7:
1582  case 10:
1583  {
1584  if( bType[1] == LibUtilities::eGLL_Lagrange)
1585  {
1586  IdxRange[1][0] = nummodes[1] - 1;
1587  IdxRange[1][1] = nummodes[1];
1588  }
1589  else
1590  {
1591  IdxRange[1][0] = 1;
1592  IdxRange[1][1] = 2;
1593  }
1594  }
1595  break;
1596  case 1:
1597  case 9:
1598  {
1599  if( bType[1] == LibUtilities::eGLL_Lagrange)
1600  {
1601  IdxRange[1][0] = 1;
1602  IdxRange[1][1] = nummodes[1] - 1;
1603 
1604  if(edgeOrient==eBackwards)
1605  {
1606  reverseOrdering = true;
1607  }
1608  }
1609  else
1610  {
1611  IdxRange[1][0] = 2;
1612  IdxRange[1][1] = nummodes[1];
1613 
1614  if(edgeOrient==eBackwards)
1615  {
1616  signChange = true;
1617  }
1618  }
1619  }
1620  break;
1621  case 3:
1622  case 11:
1623  {
1624  if( bType[1] == LibUtilities::eGLL_Lagrange)
1625  {
1626  IdxRange[1][0] = 1;
1627  IdxRange[1][1] = nummodes[1] - 1;
1628 
1629  if(edgeOrient==eForwards)
1630  {
1631  reverseOrdering = true;
1632  }
1633  }
1634  else
1635  {
1636  IdxRange[1][0] = 2;
1637  IdxRange[1][1] = nummodes[1];
1638 
1639  if(edgeOrient==eForwards)
1640  {
1641  signChange = true;
1642  }
1643  }
1644  }
1645  break;
1646  }
1647 
1648  switch(eid)
1649  {
1650  case 3:
1651  case 4:
1652  case 7:
1653  case 11:
1654  {
1655  IdxRange[0][0] = 0;
1656  IdxRange[0][1] = 1;
1657  }
1658  break;
1659  case 1:
1660  case 5:
1661  case 6:
1662  case 9:
1663  {
1664  if( bType[0] == LibUtilities::eGLL_Lagrange)
1665  {
1666  IdxRange[0][0] = nummodes[0] - 1;
1667  IdxRange[0][1] = nummodes[0];
1668  }
1669  else
1670  {
1671  IdxRange[0][0] = 1;
1672  IdxRange[0][1] = 2;
1673  }
1674  }
1675  break;
1676  case 0:
1677  case 8:
1678  {
1679  if( bType[0] == LibUtilities::eGLL_Lagrange)
1680  {
1681  IdxRange[0][0] = 1;
1682  IdxRange[0][1] = nummodes[0] - 1;
1683 
1684  if(edgeOrient==eBackwards)
1685  {
1686  reverseOrdering = true;
1687  }
1688  }
1689  else
1690  {
1691  IdxRange[0][0] = 2;
1692  IdxRange[0][1] = nummodes[0];
1693 
1694  if(edgeOrient==eBackwards)
1695  {
1696  signChange = true;
1697  }
1698  }
1699  }
1700  break;
1701  case 2:
1702  case 10:
1703  {
1704  if( bType[0] == LibUtilities::eGLL_Lagrange)
1705  {
1706  IdxRange[0][0] = 1;
1707  IdxRange[0][1] = nummodes[0] - 1;
1708 
1709  if(edgeOrient==eForwards)
1710  {
1711  reverseOrdering = true;
1712  }
1713  }
1714  else
1715  {
1716  IdxRange[0][0] = 2;
1717  IdxRange[0][1] = nummodes[0];
1718 
1719  if(edgeOrient==eForwards)
1720  {
1721  signChange = true;
1722  }
1723  }
1724  }
1725  break;
1726  }
1727 
1728  int p,q,r;
1729  int cnt = 0;
1730 
1731  for(r = IdxRange[2][0]; r < IdxRange[2][1]; r++)
1732  {
1733  for(q = IdxRange[1][0]; q < IdxRange[1][1]; q++)
1734  {
1735  for(p = IdxRange[0][0]; p < IdxRange[0][1]; p++)
1736  {
1737  maparray[cnt++]
1738  = r*nummodes[0]*nummodes[1] + q*nummodes[0] + p;
1739  }
1740  }
1741  }
1742 
1743  if( reverseOrdering )
1744  {
1745  reverse( maparray.get() , maparray.get()+nEdgeIntCoeffs );
1746  }
1747 
1748  if( signChange )
1749  {
1750  for(p = 1; p < nEdgeIntCoeffs; p+=2)
1751  {
1752  signarray[p] = -1;
1753  }
1754  }
1755  }
1756 
1757 
1758  /**
1759  * Generate mapping describing which elemental modes lie on the
1760  * interior of a given face. Accounts for face orientation.
1761  */
1762  void StdHexExp::v_GetFaceInteriorMap(const int fid,
1763  const Orientation faceOrient,
1764  Array<OneD, unsigned int> &maparray,
1765  Array<OneD, int>& signarray)
1766  {
1767  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1768  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1769  "BasisType is not a boundary interior form");
1770  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1771  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1772  "BasisType is not a boundary interior form");
1773  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
1774  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1775  "BasisType is not a boundary interior form");
1776 
1777  ASSERTL1((fid>=0)&&(fid<6),
1778  "local face id must be between 0 and 5");
1779 
1780  int nFaceIntCoeffs = GetFaceIntNcoeffs(fid);
1781 
1782  if(maparray.num_elements()!=nFaceIntCoeffs)
1783  {
1784  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1785  }
1786 
1787  if(signarray.num_elements() != nFaceIntCoeffs)
1788  {
1789  signarray = Array<OneD, int>(nFaceIntCoeffs,1);
1790  }
1791  else
1792  {
1793  fill( signarray.get() , signarray.get()+nFaceIntCoeffs, 1 );
1794  }
1795 
1796  int nummodes [3] = {m_base[0]->GetNumModes(),
1797  m_base[1]->GetNumModes(),
1798  m_base[2]->GetNumModes()};
1799 
1800  const LibUtilities::BasisType bType [3] = {GetBasisType(0),
1801  GetBasisType(1),
1802  GetBasisType(2)};
1803 
1804  int nummodesA = 0;
1805  int nummodesB = 0;
1806 
1807  // Determine the number of modes in face directions A & B based
1808  // on the face index given.
1809  switch(fid)
1810  {
1811  case 0:
1812  case 5:
1813  {
1814  nummodesA = nummodes[0];
1815  nummodesB = nummodes[1];
1816  }
1817  break;
1818  case 1:
1819  case 3:
1820  {
1821  nummodesA = nummodes[0];
1822  nummodesB = nummodes[2];
1823  }
1824  break;
1825  case 2:
1826  case 4:
1827  {
1828  nummodesA = nummodes[1];
1829  nummodesB = nummodes[2];
1830  }
1831  }
1832 
1833  int i,j;
1834  Array<OneD, int> arrayindx(nFaceIntCoeffs);
1835 
1836  // Create a mapping array to account for transposition of the
1837  // coordinates due to face orientation.
1838  for(i = 0; i < (nummodesB-2); i++)
1839  {
1840  for(j = 0; j < (nummodesA-2); j++)
1841  {
1842  if( faceOrient < eDir1FwdDir2_Dir2FwdDir1 )
1843  {
1844  arrayindx[i*(nummodesA-2)+j] = i*(nummodesA-2)+j;
1845  }
1846  else
1847  {
1848  arrayindx[i*(nummodesA-2)+j] = j*(nummodesB-2)+i;
1849  }
1850  }
1851  }
1852 
1853  int IdxRange [3][2];
1854  int Incr[3];
1855 
1856  Array<OneD, int> sign0(nummodes[0], 1);
1857  Array<OneD, int> sign1(nummodes[1], 1);
1858  Array<OneD, int> sign2(nummodes[2], 1);
1859 
1860  // Set the upper and lower bounds, and increment for the faces
1861  // involving the first coordinate direction.
1862  switch(fid)
1863  {
1864  case 0: // bottom face
1865  {
1866  IdxRange[2][0] = 0;
1867  IdxRange[2][1] = 1;
1868  Incr[2] = 1;
1869  }
1870  break;
1871  case 5: // top face
1872  {
1873  if( bType[2] == LibUtilities::eGLL_Lagrange)
1874  {
1875  IdxRange[2][0] = nummodes[2] - 1;
1876  IdxRange[2][1] = nummodes[2];
1877  Incr[2] = 1;
1878  }
1879  else
1880  {
1881  IdxRange[2][0] = 1;
1882  IdxRange[2][1] = 2;
1883  Incr[2] = 1;
1884  }
1885 
1886  }
1887  break;
1888  default: // all other faces
1889  {
1890  if( bType[2] == LibUtilities::eGLL_Lagrange)
1891  {
1892  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 2 )
1893  {
1894  IdxRange[2][0] = nummodes[2] - 2;
1895  IdxRange[2][1] = 0;
1896  Incr[2] = -1;
1897 
1898  }
1899  else
1900  {
1901  IdxRange[2][0] = 1;
1902  IdxRange[2][1] = nummodes[2] - 1;
1903  Incr[2] = 1;
1904  }
1905  }
1906  else
1907  {
1908  IdxRange[2][0] = 2;
1909  IdxRange[2][1] = nummodes[2];
1910  Incr[2] = 1;
1911 
1912  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 2 )
1913  {
1914  for(i = 3; i < nummodes[2]; i+=2)
1915  {
1916  sign2[i] = -1;
1917  }
1918  }
1919  }
1920  }
1921  }
1922 
1923  // Set the upper and lower bounds, and increment for the faces
1924  // involving the second coordinate direction.
1925  switch(fid)
1926  {
1927  case 1:
1928  {
1929  IdxRange[1][0] = 0;
1930  IdxRange[1][1] = 1;
1931  Incr[1] = 1;
1932  }
1933  break;
1934  case 3:
1935  {
1936  if( bType[1] == LibUtilities::eGLL_Lagrange)
1937  {
1938  IdxRange[1][0] = nummodes[1] - 1;
1939  IdxRange[1][1] = nummodes[1];
1940  Incr[1] = 1;
1941  }
1942  else
1943  {
1944  IdxRange[1][0] = 1;
1945  IdxRange[1][1] = 2;
1946  Incr[1] = 1;
1947  }
1948  }
1949  break;
1950  case 0:
1951  case 5:
1952  {
1953  if( bType[1] == LibUtilities::eGLL_Lagrange)
1954  {
1955  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 2 )
1956  {
1957  IdxRange[1][0] = nummodes[1] - 2;
1958  IdxRange[1][1] = 0;
1959  Incr[1] = -1;
1960 
1961  }
1962  else
1963  {
1964  IdxRange[1][0] = 1;
1965  IdxRange[1][1] = nummodes[1] - 1;
1966  Incr[1] = 1;
1967  }
1968  }
1969  else
1970  {
1971  IdxRange[1][0] = 2;
1972  IdxRange[1][1] = nummodes[1];
1973  Incr[1] = 1;
1974 
1975  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 2 )
1976  {
1977  for(i = 3; i < nummodes[1]; i+=2)
1978  {
1979  sign1[i] = -1;
1980  }
1981  }
1982  }
1983  }
1984  break;
1985  default: // case2: case4:
1986  {
1987  if( bType[1] == LibUtilities::eGLL_Lagrange)
1988  {
1989  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1 )
1990  {
1991  IdxRange[1][0] = nummodes[1] - 2;
1992  IdxRange[1][1] = 0;
1993  Incr[1] = -1;
1994 
1995  }
1996  else
1997  {
1998  IdxRange[1][0] = 1;
1999  IdxRange[1][1] = nummodes[1] - 1;
2000  Incr[1] = 1;
2001  }
2002  }
2003  else
2004  {
2005  IdxRange[1][0] = 2;
2006  IdxRange[1][1] = nummodes[1];
2007  Incr[1] = 1;
2008 
2009  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1 )
2010  {
2011  for(i = 3; i < nummodes[1]; i+=2)
2012  {
2013  sign1[i] = -1;
2014  }
2015  }
2016  }
2017  }
2018  }
2019 
2020  switch(fid)
2021  {
2022  case 4:
2023  {
2024  IdxRange[0][0] = 0;
2025  IdxRange[0][1] = 1;
2026  Incr[0] = 1;
2027  }
2028  break;
2029  case 2:
2030  {
2031  if( bType[0] == LibUtilities::eGLL_Lagrange)
2032  {
2033  IdxRange[0][0] = nummodes[0] - 1;
2034  IdxRange[0][1] = nummodes[0];
2035  Incr[0] = 1;
2036  }
2037  else
2038  {
2039  IdxRange[0][0] = 1;
2040  IdxRange[0][1] = 2;
2041  Incr[0] = 1;
2042  }
2043  }
2044  break;
2045  default:
2046  {
2047  if( bType[0] == LibUtilities::eGLL_Lagrange)
2048  {
2049  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1 )
2050  {
2051  IdxRange[0][0] = nummodes[0] - 2;
2052  IdxRange[0][1] = 0;
2053  Incr[0] = -1;
2054 
2055  }
2056  else
2057  {
2058  IdxRange[0][0] = 1;
2059  IdxRange[0][1] = nummodes[0] - 1;
2060  Incr[0] = 1;
2061  }
2062  }
2063  else
2064  {
2065  IdxRange[0][0] = 2;
2066  IdxRange[0][1] = nummodes[0];
2067  Incr[0] = 1;
2068 
2069  if( ((int) (faceOrient-eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1 )
2070  {
2071  for(i = 3; i < nummodes[0]; i+=2)
2072  {
2073  sign0[i] = -1;
2074  }
2075  }
2076  }
2077  }
2078  }
2079 
2080  int p,q,r;
2081  int cnt = 0;
2082 
2083  for(r = IdxRange[2][0]; r != IdxRange[2][1]; r+=Incr[2])
2084  {
2085  for(q = IdxRange[1][0]; q != IdxRange[1][1]; q+=Incr[1])
2086  {
2087  for(p = IdxRange[0][0]; p != IdxRange[0][1]; p+=Incr[0])
2088  {
2089  maparray [ arrayindx[cnt ] ]
2090  = r*nummodes[0]*nummodes[1] + q*nummodes[0] + p;
2091  signarray[ arrayindx[cnt++] ]
2092  = sign0[p] * sign1[q] * sign2[r];
2093  }
2094  }
2095  }
2096  }
2097 
2098 
2099  /**
2100  * @param outarray Storage area for computed map.
2101  */
2102  void StdHexExp::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
2103  {
2104  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
2105  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
2106  "BasisType is not a boundary interior form");
2107  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
2108  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
2109  "BasisType is not a boundary interior form");
2110  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
2111  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
2112  "BasisType is not a boundary interior form");
2113 
2114  int i;
2115  int nummodes [3] = {m_base[0]->GetNumModes(),
2116  m_base[1]->GetNumModes(),
2117  m_base[2]->GetNumModes()};
2118 
2119  int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
2120 
2121  if(outarray.num_elements() != nIntCoeffs)
2122  {
2123  outarray = Array<OneD, unsigned int>(nIntCoeffs);
2124  }
2125 
2126  const LibUtilities::BasisType Btype [3] = {GetBasisType(0),
2127  GetBasisType(1),
2128  GetBasisType(2)};
2129 
2130  int p,q,r;
2131  int cnt = 0;
2132 
2133  int IntIdx [3][2];
2134 
2135  for(i = 0; i < 3; i++)
2136  {
2137  if( Btype[i] == LibUtilities::eModified_A)
2138  {
2139  IntIdx[i][0] = 2;
2140  IntIdx[i][1] = nummodes[i];
2141  }
2142  else
2143  {
2144  IntIdx[i][0] = 1;
2145  IntIdx[i][1] = nummodes[i]-1;
2146  }
2147  }
2148 
2149  for(r = IntIdx[2][0]; r < IntIdx[2][1]; r++)
2150  {
2151  for( q = IntIdx[1][0]; q < IntIdx[1][1]; q++)
2152  {
2153  for( p = IntIdx[0][0]; p < IntIdx[0][1]; p++)
2154  {
2155  outarray[cnt++] = r*nummodes[0]*nummodes[1] +
2156  q*nummodes[0] + p;
2157  }
2158  }
2159  }
2160  }
2161 
2162 
2163  /**
2164  * @param outarray Storage for computed map.
2165  */
2166  void StdHexExp::v_GetBoundaryMap(Array<OneD, unsigned int>& outarray)
2167  {
2168  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
2169  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
2170  "BasisType is not a boundary interior form");
2171  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
2172  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
2173  "BasisType is not a boundary interior form");
2174  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
2175  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
2176  "BasisType is not a boundary interior form");
2177 
2178  int i;
2179  int nummodes [3] = {m_base[0]->GetNumModes(),
2180  m_base[1]->GetNumModes(),
2181  m_base[2]->GetNumModes()};
2182 
2183  int nBndCoeffs = NumBndryCoeffs();
2184 
2185  if(outarray.num_elements()!=nBndCoeffs)
2186  {
2187  outarray = Array<OneD, unsigned int>(nBndCoeffs);
2188  }
2189 
2190  const LibUtilities::BasisType Btype [3] = {GetBasisType(0),
2191  GetBasisType(1),
2192  GetBasisType(2)};
2193 
2194  int p,q,r;
2195  int cnt = 0;
2196 
2197  int BndIdx [3][2];
2198  int IntIdx [3][2];
2199 
2200  for(i = 0; i < 3; i++)
2201  {
2202  BndIdx[i][0] = 0;
2203 
2204  if( Btype[i] == LibUtilities::eModified_A)
2205  {
2206  BndIdx[i][1] = 1;
2207  IntIdx[i][0] = 2;
2208  IntIdx[i][1] = nummodes[i];
2209  }
2210  else
2211  {
2212  BndIdx[i][1] = nummodes[i]-1;
2213  IntIdx[i][0] = 1;
2214  IntIdx[i][1] = nummodes[i]-1;
2215  }
2216  }
2217 
2218 
2219  for(i = 0; i < 2; i++)
2220  {
2221  r = BndIdx[2][i];
2222  for( q = 0; q < nummodes[1]; q++)
2223  {
2224  for( p = 0; p < nummodes[0]; p++)
2225  {
2226  outarray[cnt++] = r*nummodes[0]*nummodes[1]+q*nummodes[0] + p;
2227  }
2228  }
2229  }
2230 
2231  for(r = IntIdx[2][0]; r < IntIdx[2][1]; r++)
2232  {
2233  for( i = 0; i < 2; i++)
2234  {
2235  q = BndIdx[1][i];
2236  for( p = 0; p < nummodes[0]; p++)
2237  {
2238  outarray[cnt++] = r*nummodes[0]*nummodes[1] +
2239  q*nummodes[0] + p;
2240  }
2241  }
2242 
2243  for( q = IntIdx[1][0]; q < IntIdx[1][1]; q++)
2244  {
2245  for( i = 0; i < 2; i++)
2246  {
2247  p = BndIdx[0][i];
2248  outarray[cnt++] = r*nummodes[0]*nummodes[1] +
2249  q*nummodes[0] + p;
2250  }
2251  }
2252  }
2253 
2254  sort(outarray.get(), outarray.get() + nBndCoeffs);
2255  }
2256 
2257  DNekMatSharedPtr StdHexExp::v_GenMatrix(const StdMatrixKey &mkey)
2258  {
2259  return StdExpansion::CreateGeneralMatrix(mkey);
2260  }
2261 
2262 
2263  DNekMatSharedPtr StdHexExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
2264  {
2265  return StdExpansion::CreateGeneralMatrix(mkey);
2266  }
2267 
2268 
2269  void StdHexExp::v_MassMatrixOp(
2270  const Array<OneD, const NekDouble> &inarray,
2271  Array<OneD,NekDouble> &outarray,
2272  const StdMatrixKey &mkey)
2273  {
2274  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
2275  }
2276 
2277 
2278  void StdHexExp::v_LaplacianMatrixOp(
2279  const Array<OneD, const NekDouble> &inarray,
2280  Array<OneD,NekDouble> &outarray,
2281  const StdMatrixKey &mkey)
2282  {
2283  StdHexExp::v_LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
2284  }
2285 
2286 
2287  void StdHexExp::v_LaplacianMatrixOp(const int k1, const int k2,
2288  const Array<OneD, const NekDouble> &inarray,
2289  Array<OneD,NekDouble> &outarray,
2290  const StdMatrixKey &mkey)
2291  {
2292  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,
2293  mkey);
2294  }
2295 
2296 
2297  void StdHexExp::v_WeakDerivMatrixOp(const int i,
2298  const Array<OneD, const NekDouble> &inarray,
2299  Array<OneD,NekDouble> &outarray,
2300  const StdMatrixKey &mkey)
2301  {
2302  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,
2303  mkey);
2304  }
2305 
2306  void StdHexExp::v_HelmholtzMatrixOp(
2307  const Array<OneD, const NekDouble> &inarray,
2308  Array<OneD,NekDouble> &outarray,
2309  const StdMatrixKey &mkey)
2310  {
2311  StdHexExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
2312  }
2313 
2314 
2315  void StdHexExp::v_GeneralMatrixOp_MatOp(
2316  const Array<OneD, const NekDouble> &inarray,
2317  Array<OneD,NekDouble> &outarray,
2318  const StdMatrixKey &mkey)
2319  {
2320  DNekMatSharedPtr mat = m_stdMatrixManager[mkey];
2321 
2322  if(inarray.get() == outarray.get())
2323  {
2324  Array<OneD,NekDouble> tmp(m_ncoeffs);
2325  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
2326 
2327  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
2328  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
2329  }
2330  else
2331  {
2332  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
2333  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
2334  }
2335  }
2336 
2337 
2338  void StdHexExp::v_MultiplyByStdQuadratureMetric(const Array<OneD, const NekDouble>& inarray,
2339  Array<OneD, NekDouble> &outarray)
2340  {
2341  int i;
2342  int nquad0 = m_base[0]->GetNumPoints();
2343  int nquad1 = m_base[1]->GetNumPoints();
2344  int nquad2 = m_base[2]->GetNumPoints();
2345  int nq01 = nquad0*nquad1;
2346  int nq12 = nquad1*nquad2;
2347 
2348  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
2349  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
2350  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
2351 
2352  for(i = 0; i < nq12; ++i)
2353  {
2354  Vmath::Vmul(nquad0, inarray.get()+i*nquad0, 1,
2355  w0.get(), 1, outarray.get()+i*nquad0,1);
2356  }
2357 
2358  for(i = 0; i < nq12; ++i)
2359  {
2360  Vmath::Smul(nquad0, w1[i%nquad1], outarray.get()+i*nquad0, 1,
2361  outarray.get()+i*nquad0, 1);
2362  }
2363 
2364  for(i = 0; i < nquad2; ++i)
2365  {
2366  Vmath::Smul(nq01, w2[i], outarray.get()+i*nq01, 1,
2367  outarray.get()+i*nq01, 1);
2368  }
2369  }
2370 
2371  void StdHexExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
2372  const StdMatrixKey &mkey)
2373  {
2374  // Generate an orthonogal expansion
2375  int qa = m_base[0]->GetNumPoints();
2376  int qb = m_base[1]->GetNumPoints();
2377  int qc = m_base[2]->GetNumPoints();
2378  int nmodes_a = m_base[0]->GetNumModes();
2379  int nmodes_b = m_base[1]->GetNumModes();
2380  int nmodes_c = m_base[2]->GetNumModes();
2381  // Declare orthogonal basis.
2382  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
2383  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
2384  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
2385 
2386  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
2387  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A,nmodes_b,pb);
2388  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_A,nmodes_c,pc);
2389  StdHexExp OrthoExp(Ba,Bb,Bc);
2390 
2391  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2392  int i,j,k,cnt=0;
2393 
2394  // project onto modal space.
2395  OrthoExp.FwdTrans(array,orthocoeffs);
2396 
2397  if(mkey.ConstFactorExists(eFactorSVVPowerKerDiffCoeff))
2398  {
2399  // Rodrigo's power kernel
2400  NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
2401  NekDouble SvvDiffCoeff =
2402  mkey.GetConstFactor(eFactorSVVPowerKerDiffCoeff)*
2403  mkey.GetConstFactor(eFactorSVVDiffCoeff);
2404 
2405  for(int i = 0; i < nmodes_a; ++i)
2406  {
2407  for(int j = 0; j < nmodes_b; ++j)
2408  {
2409  NekDouble fac1 = std::max(
2410  pow((1.0*i)/(nmodes_a-1),cutoff*nmodes_a),
2411  pow((1.0*j)/(nmodes_b-1),cutoff*nmodes_b));
2412 
2413  for(int k = 0; k < nmodes_c; ++k)
2414  {
2415  NekDouble fac = std::max(fac1,
2416  pow((1.0*k)/(nmodes_c-1),cutoff*nmodes_c));
2417 
2418  orthocoeffs[cnt]
2419  *= SvvDiffCoeff * fac;
2420  cnt++;
2421  }
2422  }
2423  }
2424  }
2425  else if(mkey.ConstFactorExists(eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
2426  {
2427  NekDouble SvvDiffCoeff =
2428  mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff)*
2429  mkey.GetConstFactor(eFactorSVVDiffCoeff);
2430 
2431  int max_abc = max(nmodes_a-kSVVDGFiltermodesmin,
2432  nmodes_b-kSVVDGFiltermodesmin);
2433  max_abc = max(max_abc, nmodes_c-kSVVDGFiltermodesmin);
2434  // clamp max_abc
2435  max_abc = max(max_abc,0);
2436  max_abc = min(max_abc,kSVVDGFiltermodesmax-kSVVDGFiltermodesmin);
2437 
2438  for(int i = 0; i < nmodes_a; ++i)
2439  {
2440  for(int j = 0; j < nmodes_b; ++j)
2441  {
2442  int maxij = max(i,j);
2443 
2444  for(int k = 0; k < nmodes_c; ++k)
2445  {
2446  int maxijk = max(maxij,k);
2447  maxijk = min(maxijk,kSVVDGFiltermodesmax-1);
2448 
2449  orthocoeffs[cnt] *= SvvDiffCoeff *
2450  kSVVDGFilter[max_abc][maxijk];
2451  cnt++;
2452  }
2453  }
2454  }
2455  }
2456  else
2457  {
2458 
2459  int cutoff = (int) (mkey.GetConstFactor(eFactorSVVCutoffRatio)*min(nmodes_a,nmodes_b));
2460  NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
2461  // Filter just trilinear space
2462  int nmodes = max(nmodes_a,nmodes_b);
2463  nmodes = max(nmodes,nmodes_c);
2464 
2465  Array<OneD, NekDouble> fac(nmodes,1.0);
2466  for(j = cutoff; j < nmodes; ++j)
2467  {
2468  fac[j] = fabs((j-nmodes)/((NekDouble) (j-cutoff+1.0)));
2469  fac[j] *= fac[j]; //added this line to conform with equation
2470  }
2471 
2472  for(i = 0; i < nmodes_a; ++i)
2473  {
2474  for(j = 0; j < nmodes_b; ++j)
2475  {
2476  for(k = 0; k < nmodes_c; ++k)
2477  {
2478  if((i >= cutoff)||(j >= cutoff)||(k >= cutoff))
2479  {
2480  orthocoeffs[i*nmodes_a*nmodes_b +
2481  j*nmodes_c + k] *=
2482  (SvvDiffCoeff*exp(-(fac[i]+fac[j]+fac[k])));
2483  }
2484  else
2485  {
2486  orthocoeffs[i*nmodes_a*nmodes_b + j*nmodes_c + k] *= 0.0;
2487  }
2488  }
2489  }
2490  }
2491  }
2492 
2493  // backward transform to physical space
2494  OrthoExp.BwdTrans(orthocoeffs,array);
2495  }
2496 
2497  void StdHexExp::v_ExponentialFilter(
2498  Array<OneD, NekDouble> &array,
2499  const NekDouble alpha,
2500  const NekDouble exponent,
2501  const NekDouble cutoff)
2502  {
2503  // Generate an orthogonal expansion
2504  int qa = m_base[0]->GetNumPoints();
2505  int qb = m_base[1]->GetNumPoints();
2506  int qc = m_base[2]->GetNumPoints();
2507  int nmodesA = m_base[0]->GetNumModes();
2508  int nmodesB = m_base[1]->GetNumModes();
2509  int nmodesC = m_base[2]->GetNumModes();
2510  int P = nmodesA - 1;
2511  int Q = nmodesB - 1;
2512  int R = nmodesC - 1;
2513 
2514  // Declare orthogonal basis.
2515  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
2516  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
2517  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
2518 
2519  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, nmodesA, pa);
2520  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, nmodesB, pb);
2521  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_A, nmodesC, pc);
2522  StdHexExp OrthoExp(Ba,Bb,Bc);
2523 
2524  // Cutoff
2525  int Pcut = cutoff*P;
2526  int Qcut = cutoff*Q;
2527  int Rcut = cutoff*R;
2528 
2529  // Project onto orthogonal space.
2530  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2531  OrthoExp.FwdTrans(array,orthocoeffs);
2532 
2533  //
2534  NekDouble fac, fac1, fac2, fac3;
2535  int index = 0;
2536  for(int i = 0; i < nmodesA; ++i)
2537  {
2538  for(int j = 0; j < nmodesB; ++j)
2539  {
2540  for(int k = 0; k < nmodesC; ++k, ++index)
2541  {
2542  //to filter out only the "high-modes"
2543  if(i > Pcut || j > Qcut || k > Rcut)
2544  {
2545  fac1 = (NekDouble) (i - Pcut)/( (NekDouble)(P - Pcut) );
2546  fac2 = (NekDouble) (j - Qcut)/( (NekDouble)(Q - Qcut) );
2547  fac3 = (NekDouble) (k - Rcut)/( (NekDouble)(R - Rcut) );
2548  fac = max( max(fac1, fac2), fac3);
2549  fac = pow(fac, exponent);
2550  orthocoeffs[index] *= exp(-alpha*fac);
2551  }
2552  }
2553  }
2554  }
2555 
2556  // backward transform to physical space
2557  OrthoExp.BwdTrans(orthocoeffs,array);
2558  }
2559 
2560  }
2561 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:274
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:60
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::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::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::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
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::StdExpansion
The base class for all shapes.
Definition: StdExpansion.h:69
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::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:140
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::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::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:714
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::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::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::CreateGeneralMatrix
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
Definition: StdExpansion.cpp:324
Nektar::StdRegions::StdExpansion::GetFaceIntNcoeffs
int GetFaceIntNcoeffs(const int i) const
Definition: StdExpansion.h:358
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::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::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:469
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::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::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::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
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::m_stdMatrixManager
LibUtilities::NekManager< StdMatrixKey, DNekMat, StdMatrixKey::opLess > m_stdMatrixManager
Definition: StdExpansion.h:1432
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::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
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1429
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
Class representing a hexehedral element in reference space.
Definition: StdHexExp.h:48
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::StdHexExp::v_ExponentialFilter
virtual void v_ExponentialFilter(Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff)
Definition: StdHexExp.cpp:2497
Nektar::StdRegions::StdHexExp::v_GetNfaces
virtual int v_GetNfaces() const
Definition: StdHexExp.cpp:678
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::StdHexExp::v_NumBndryCoeffs
virtual int v_NumBndryCoeffs() const
Definition: StdHexExp.cpp:690
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:2306
Nektar::StdRegions::StdHexExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2371
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_FillMode
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:618
Nektar::StdRegions::StdHexExp::v_GetTotalEdgeIntNcoeffs
virtual int v_GetTotalEdgeIntNcoeffs() const
Definition: StdHexExp.cpp:749
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
Nektar::StdRegions::StdHexExp::v_GetEdgeNcoeffs
virtual int v_GetEdgeNcoeffs(const int i) const
Definition: StdHexExp.cpp:731
Nektar::StdRegions::StdHexExp::v_GetBoundaryMap
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdHexExp.cpp:2166
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::StdRegions::StdHexExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2257
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:1464
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::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:2315
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:1762
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::StdRegions::StdHexExp::v_GetTotalFaceIntNcoeffs
virtual int v_GetTotalFaceIntNcoeffs() const
Definition: StdHexExp.cpp:791
Nektar::StdRegions::StdHexExp::v_GetVertexMap
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdHexExp.cpp:1330
Nektar::StdRegions::StdHexExp::v_GetNverts
virtual int v_GetNverts() const
Definition: StdHexExp.cpp:666
Nektar::StdRegions::StdHexExp::v_GetFaceIntNcoeffs
virtual int v_GetFaceIntNcoeffs(const int i) const
Definition: StdHexExp.cpp:773
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:2269
Nektar::StdRegions::StdHexExp::v_DetFaceBasisKey
virtual const LibUtilities::BasisKey v_DetFaceBasisKey(const int i, const int k) const
Definition: StdHexExp.cpp:848
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
Nektar::StdRegions::StdHexExp::v_GetNedges
virtual int v_GetNedges() const
Definition: StdHexExp.cpp:672
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:895
Nektar::StdRegions::StdHexExp::v_GetFaceNumModes
virtual void v_GetFaceNumModes(const int fid, const Orientation faceOrient, int &numModes0, int &numModes1)
Definition: StdHexExp.cpp:921
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:2278
Nektar::StdRegions::StdHexExp::v_IsBoundaryInteriorExpansion
virtual bool v_IsBoundaryInteriorExpansion()
Definition: StdHexExp.cpp:75
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:969
Nektar::StdRegions::StdHexExp::v_GetFacePointsKey
virtual LibUtilities::PointsKey v_GetFacePointsKey(const int i, const int j) const
Definition: StdHexExp.cpp:819
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::StdHexExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2263
Nektar::StdRegions::StdHexExp::v_MultiplyByStdQuadratureMetric
virtual void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:2338
Nektar::StdRegions::StdHexExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdHexExp.cpp:2102
Nektar::StdRegions::StdHexExp::v_NumDGBndryCoeffs
virtual int v_NumDGBndryCoeffs() const
Definition: StdHexExp.cpp:711
Nektar::StdRegions::StdHexExp::v_CalcNumberOfCoefficients
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdHexExp.cpp:839
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_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_DetShapeType
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdHexExp.cpp:684
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:2297
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::StdHexExp::v_GetFaceNcoeffs
virtual int v_GetFaceNcoeffs(const int i) const
Definition: StdHexExp.cpp:755
Nektar::StdRegions::StdHexExp::v_GetEdgeBasisType
virtual LibUtilities::BasisType v_GetEdgeBasisType(const int i) const
Definition: StdHexExp.cpp:877
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
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::StdHexExp::v_GetFaceNumPoints
virtual int v_GetFaceNumPoints(const int i) const
Definition: StdHexExp.cpp:798
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:121
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:130
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
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
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
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::LibUtilities::eModified_B
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::LibUtilities::eOrtho_A
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:45
Nektar::LibUtilities::eModified_C
@ eModified_C
Principle Modified Functions .
Definition: BasisType.h:50
Nektar::LibUtilities::eGLL_Lagrange
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:54
Nektar::LibUtilities::eOrtho_C
@ eOrtho_C
Principle Orthogonal Functions .
Definition: BasisType.h:47
Nektar::LibUtilities::eOrtho_B
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:46
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
Nektar::StdRegions::kSVVDGFiltermodesmin
const int kSVVDGFiltermodesmin
Definition: StdRegions.hpp:384
Nektar::StdRegions::kSVVDGFiltermodesmax
const int kSVVDGFiltermodesmax
Definition: StdRegions.hpp:385
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::kSVVDGFilter
const NekDouble kSVVDGFilter[9][11]
Definition: StdRegions.hpp:387
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:787
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
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
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
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
main.P
P
Definition: main.py:133
Generated on Mon Aug 23 2021 06:57:50 for Nektar++ by doxygen 1.9.1