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