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