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 // 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 nummodesA,
938  int nummodesB)
939  {
940  int i,j;
941  const int nummodes0 = m_base[0]->GetNumModes();
942  const int nummodes1 = m_base[1]->GetNumModes();
943  const int nummodes2 = m_base[2]->GetNumModes();
944 
945  ASSERTL1(GetEdgeBasisType(0) == GetEdgeBasisType(1) &&
946  GetEdgeBasisType(0) == GetEdgeBasisType(2),
947  "Method only implemented if BasisType is indentical in "
948  "all directions");
949  ASSERTL1(GetEdgeBasisType(0) == LibUtilities::eModified_A ||
950  GetEdgeBasisType(0) == LibUtilities::eGLL_Lagrange,
951  "Method only implemented for Modified_A or GLL_Lagrange BasisType");
952 
953  if (nummodesA == -1)
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 
975 
976  bool modified = GetEdgeBasisType(0) == LibUtilities::eModified_A;
977  int nFaceCoeffs = nummodesA*nummodesB;
978 
979  if(maparray.num_elements() != nFaceCoeffs)
980  {
981  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
982  }
983 
984  if(signarray.num_elements() != nFaceCoeffs)
985  {
986  signarray = Array<OneD, int>(nFaceCoeffs,1);
987  }
988  else
989  {
990  fill( signarray.get() , signarray.get()+nFaceCoeffs, 1 );
991  }
992 
993  Array<OneD, int> arrayindx(nFaceCoeffs);
994 
995  for(i = 0; i < nummodesB; i++)
996  {
997  for(j = 0; j < nummodesA; j++)
998  {
999  if( faceOrient < 9 )
1000  {
1001  arrayindx[i*nummodesA+j] = i*nummodesA+j;
1002  }
1003  else
1004  {
1005  arrayindx[i*nummodesA+j] = j*nummodesB+i;
1006  }
1007  }
1008  }
1009 
1010  int offset = 0;
1011  int jump1 = 1;
1012  int jump2 = 1;
1013 
1014  switch(fid)
1015  {
1016  case 5:
1017  {
1018  if (modified)
1019  {
1020  offset = nummodes0*nummodes1;
1021  }
1022  else
1023  {
1024  offset = (nummodes2-1)*nummodes0*nummodes1;
1025  jump1 = nummodes0;
1026  }
1027  }
1028  case 0:
1029  {
1030  jump1 = nummodes0;
1031  }
1032  break;
1033  case 3:
1034  {
1035  if (modified)
1036  {
1037  offset = nummodes0;
1038  }
1039  else
1040  {
1041  offset = nummodes0*(nummodes1-1);
1042  jump1 = nummodes0*nummodes1;
1043  }
1044  }
1045  case 1:
1046  {
1047  jump1 = nummodes0*nummodes1;
1048  }
1049  break;
1050  case 2:
1051  {
1052  if (modified)
1053  {
1054  offset = 1;
1055  }
1056  else
1057  {
1058  offset = nummodes0-1;
1059  jump1 = nummodes0*nummodes1;
1060  jump2 = nummodes0;
1061 
1062  }
1063  }
1064  case 4:
1065  {
1066  jump1 = nummodes0*nummodes1;
1067  jump2 = nummodes0;
1068  }
1069  break;
1070  default:
1071  ASSERTL0(false,"fid must be between 0 and 5");
1072  }
1073 
1074  for(i = 0; i < nummodesB; i++)
1075  {
1076  for(j = 0; j < nummodesA; j++)
1077  {
1078  maparray[ arrayindx[i*nummodesA+j] ]
1079  = i*jump1 + j*jump2 + offset;
1080  }
1081  }
1082 
1083 
1084 
1085  if( (faceOrient==6) || (faceOrient==8) ||
1086  (faceOrient==11) || (faceOrient==12) )
1087  {
1088  if(faceOrient<9)
1089  {
1090  if (modified)
1091  {
1092  for(i = 3; i < nummodesB; i+=2)
1093  {
1094  for(j = 0; j < nummodesA; j++)
1095  {
1096  signarray[ arrayindx[i*nummodesA+j] ] *= -1;
1097  }
1098  }
1099 
1100  for(i = 0; i < nummodesA; i++)
1101  {
1102  swap( maparray[i] , maparray[i+nummodesA] );
1103  swap( signarray[i] , signarray[i+nummodesA] );
1104  }
1105 
1106  }
1107  else
1108  {
1109  for(i = 0; i < nummodesA; i++)
1110  {
1111  for(j = 0; j < nummodesB/2; j++)
1112  {
1113  swap( maparray[i + j*nummodesA],
1114  maparray[i+nummodesA*nummodesB
1115  -nummodesA -j*nummodesA] );
1116  swap( signarray[i + j*nummodesA],
1117  signarray[i+nummodesA*nummodesB
1118  -nummodesA -j*nummodesA]);
1119  }
1120  }
1121  }
1122  }
1123  else
1124  {
1125  if (modified)
1126  {
1127  for(i = 0; i < nummodesB; i++)
1128  {
1129  for(j = 3; j < nummodesA; j+=2)
1130  {
1131  signarray[ arrayindx[i*nummodesA+j] ] *= -1;
1132  }
1133  }
1134 
1135  for(i = 0; i < nummodesB; i++)
1136  {
1137  swap( maparray[i] , maparray[i+nummodesB] );
1138  swap( signarray[i] , signarray[i+nummodesB] );
1139  }
1140 
1141  }
1142  else
1143  {
1144  for(i = 0; i < nummodesA; i++)
1145  {
1146  for(j = 0; j < nummodesB/2; j++)
1147  {
1148  swap( maparray[i*nummodesB + j],
1149  maparray[i*nummodesB + nummodesB -1 -j]);
1150  swap( signarray[i*nummodesB + j],
1151  signarray[i*nummodesB + nummodesB -1 -j]);
1152  }
1153  }
1154  }
1155  }
1156  }
1157 
1158  if( (faceOrient==7) || (faceOrient==8) ||
1159  (faceOrient==10) || (faceOrient==12) )
1160  {
1161  if(faceOrient<9)
1162  {
1163  if (modified)
1164  {
1165  for(i = 0; i < nummodesB; i++)
1166  {
1167  for(j = 3; j < nummodesA; j+=2)
1168  {
1169  signarray[ arrayindx[i*nummodesA+j] ] *= -1;
1170  }
1171  }
1172 
1173  for(i = 0; i < nummodesB; i++)
1174  {
1175  swap( maparray[i*nummodesA],
1176  maparray[i*nummodesA+1]);
1177  swap( signarray[i*nummodesA],
1178  signarray[i*nummodesA+1]);
1179  }
1180  }
1181  else
1182  {
1183  for(i = 0; i < nummodesB; i++)
1184  {
1185  for(j = 0; j < nummodesA/2; j++)
1186  {
1187  swap( maparray[i*nummodesA + j],
1188  maparray[i*nummodesA + nummodesA -1 -j]);
1189  swap( signarray[i*nummodesA + j],
1190  signarray[i*nummodesA + nummodesA -1 -j]);
1191  }
1192  }
1193  }
1194 
1195 
1196 
1197  }
1198  else
1199  {
1200  if (modified)
1201  {
1202  for(i = 3; i < nummodesB; i+=2)
1203  {
1204  for(j = 0; j < nummodesA; j++)
1205  {
1206  signarray[ arrayindx[i*nummodesA+j] ] *= -1;
1207  }
1208  }
1209 
1210  for(i = 0; i < nummodesA; i++)
1211  {
1212  swap( maparray[i*nummodesB],
1213  maparray[i*nummodesB+1]);
1214  swap( signarray[i*nummodesB],
1215  signarray[i*nummodesB+1]);
1216  }
1217  }
1218  else
1219  {
1220  for(i = 0; i < nummodesB; i++)
1221  {
1222  for(j = 0; j < nummodesA/2; j++)
1223  {
1224  swap( maparray[i + j*nummodesB] ,
1225  maparray[i+nummodesA*nummodesB -
1226  nummodesB -j*nummodesB] );
1227  swap( signarray[i + j*nummodesB] ,
1228  signarray[i+nummodesA*nummodesB -
1229  nummodesB -j*nummodesB] );
1230 
1231  }
1232  }
1233 
1234  }
1235  }
1236  }
1237  }
1238 
1239 
1240 
1241  /**
1242  * Expansions in each of the three dimensions must be of type
1243  * LibUtilities#eModified_A or LibUtilities#eGLL_Lagrange.
1244  *
1245  * @param localVertexId ID of vertex (0..7)
1246  * @returns Position of vertex in local numbering scheme.
1247  */
1248  int StdHexExp::v_GetVertexMap(const int localVertexId, bool useCoeffPacking)
1249  {
1250  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1251  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1252  "BasisType is not a boundary interior form");
1253  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1254  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1255  "BasisType is not a boundary interior form");
1256  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
1257  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1258  "BasisType is not a boundary interior form");
1259 
1260  ASSERTL1((localVertexId>=0)&&(localVertexId<8),
1261  "local vertex id must be between 0 and 7");
1262 
1263  int p = 0;
1264  int q = 0;
1265  int r = 0;
1266 
1267  // Retrieve the number of modes in each dimension.
1268  int nummodes [3] = {m_base[0]->GetNumModes(),
1269  m_base[1]->GetNumModes(),
1270  m_base[2]->GetNumModes()};
1271 
1272  if(useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1273  {
1274  if(localVertexId > 3)
1275  {
1276  if( GetBasisType(2) == LibUtilities::eGLL_Lagrange)
1277  {
1278  r = nummodes[2]-1;
1279  }
1280  else
1281  {
1282  r = 1;
1283  }
1284  }
1285 
1286  switch(localVertexId % 4)
1287  {
1288  case 0:
1289  break;
1290  case 1:
1291  {
1292  if( GetBasisType(0) == LibUtilities::eGLL_Lagrange)
1293  {
1294  p = nummodes[0]-1;
1295  }
1296  else
1297  {
1298  p = 1;
1299  }
1300  }
1301  break;
1302  case 2:
1303  {
1304  if( GetBasisType(1) == LibUtilities::eGLL_Lagrange)
1305  {
1306  q = nummodes[1]-1;
1307  }
1308  else
1309  {
1310  q = 1;
1311  }
1312  }
1313  break;
1314  case 3:
1315  {
1316  if( GetBasisType(1) == LibUtilities::eGLL_Lagrange)
1317  {
1318  p = nummodes[0]-1;
1319  q = nummodes[1]-1;
1320  }
1321  else
1322  {
1323  p = 1;
1324  q = 1;
1325  }
1326  }
1327  break;
1328  }
1329  }
1330  else
1331  {
1332  // Right face (vertices 1,2,5,6)
1333  if( (localVertexId % 4) % 3 > 0 )
1334  {
1335  if( GetBasisType(0) == LibUtilities::eGLL_Lagrange)
1336  {
1337  p = nummodes[0]-1;
1338  }
1339  else
1340  {
1341  p = 1;
1342  }
1343  }
1344 
1345  // Back face (vertices 2,3,6,7)
1346  if( localVertexId % 4 > 1 )
1347  {
1348  if( GetBasisType(1) == LibUtilities::eGLL_Lagrange)
1349  {
1350  q = nummodes[1]-1;
1351  }
1352  else
1353  {
1354  q = 1;
1355  }
1356  }
1357 
1358  // Top face (vertices 4,5,6,7)
1359  if( localVertexId > 3)
1360  {
1361  if( GetBasisType(2) == LibUtilities::eGLL_Lagrange)
1362  {
1363  r = nummodes[2]-1;
1364  }
1365  else
1366  {
1367  r = 1;
1368  }
1369  }
1370  }
1371  // Compute the local number.
1372  return r*nummodes[0]*nummodes[1] + q*nummodes[0] + p;
1373  }
1374 
1375 
1376  /**
1377  * @param eid The edge to compute the numbering for.
1378  * @param edgeOrient Orientation of the edge.
1379  * @param maparray Storage for computed mapping array.
1380  * @param signarray ?
1381  */
1382  void StdHexExp::v_GetEdgeInteriorMap(const int eid,
1383  const Orientation edgeOrient,
1384  Array<OneD, unsigned int> &maparray,
1385  Array<OneD, int> &signarray)
1386  {
1387  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1388  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1389  "BasisType is not a boundary interior form");
1390  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1391  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1392  "BasisType is not a boundary interior form");
1393  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
1394  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1395  "BasisType is not a boundary interior form");
1396 
1397  ASSERTL1((eid>=0)&&(eid<12),
1398  "local edge id must be between 0 and 11");
1399 
1400  int nEdgeIntCoeffs = GetEdgeNcoeffs(eid)-2;
1401 
1402  if(maparray.num_elements()!=nEdgeIntCoeffs)
1403  {
1404  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1405  }
1406 
1407  if(signarray.num_elements() != nEdgeIntCoeffs)
1408  {
1409  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
1410  }
1411  else
1412  {
1413  fill( signarray.get() , signarray.get()+nEdgeIntCoeffs, 1 );
1414  }
1415 
1416  int nummodes [3] = {m_base[0]->GetNumModes(),
1417  m_base[1]->GetNumModes(),
1418  m_base[2]->GetNumModes()};
1419 
1420  const LibUtilities::BasisType bType [3] = {GetBasisType(0),
1421  GetBasisType(1),
1422  GetBasisType(2)};
1423 
1424  bool reverseOrdering = false;
1425  bool signChange = false;
1426 
1427  int IdxRange [3][2] = {{0,0},{0,0},{0,0}};
1428 
1429  switch(eid)
1430  {
1431  case 0:
1432  case 1:
1433  case 2:
1434  case 3:
1435  {
1436  IdxRange[2][0] = 0;
1437  IdxRange[2][1] = 1;
1438  }
1439  break;
1440  case 8:
1441  case 9:
1442  case 10:
1443  case 11:
1444  {
1445  if( bType[2] == LibUtilities::eGLL_Lagrange)
1446  {
1447  IdxRange[2][0] = nummodes[2] - 1;
1448  IdxRange[2][1] = nummodes[2];
1449  }
1450  else
1451  {
1452  IdxRange[2][0] = 1;
1453  IdxRange[2][1] = 2;
1454  }
1455  }
1456  break;
1457  case 4:
1458  case 5:
1459  case 6:
1460  case 7:
1461  {
1462  if( bType[2] == LibUtilities::eGLL_Lagrange)
1463  {
1464  IdxRange[2][0] = 1;
1465  IdxRange[2][1] = nummodes[2] - 1;
1466 
1467  if(edgeOrient==eBackwards)
1468  {
1469  reverseOrdering = true;
1470  }
1471  }
1472  else
1473  {
1474  IdxRange[2][0] = 2;
1475  IdxRange[2][1] = nummodes[2];
1476 
1477  if(edgeOrient==eBackwards)
1478  {
1479  signChange = true;
1480  }
1481  }
1482  }
1483  break;
1484  }
1485 
1486  switch(eid)
1487  {
1488  case 0:
1489  case 4:
1490  case 5:
1491  case 8:
1492  {
1493  IdxRange[1][0] = 0;
1494  IdxRange[1][1] = 1;
1495  }
1496  break;
1497  case 2:
1498  case 6:
1499  case 7:
1500  case 10:
1501  {
1502  if( bType[1] == LibUtilities::eGLL_Lagrange)
1503  {
1504  IdxRange[1][0] = nummodes[1] - 1;
1505  IdxRange[1][1] = nummodes[1];
1506  }
1507  else
1508  {
1509  IdxRange[1][0] = 1;
1510  IdxRange[1][1] = 2;
1511  }
1512  }
1513  break;
1514  case 1:
1515  case 9:
1516  {
1517  if( bType[1] == LibUtilities::eGLL_Lagrange)
1518  {
1519  IdxRange[1][0] = 1;
1520  IdxRange[1][1] = nummodes[1] - 1;
1521 
1522  if(edgeOrient==eBackwards)
1523  {
1524  reverseOrdering = true;
1525  }
1526  }
1527  else
1528  {
1529  IdxRange[1][0] = 2;
1530  IdxRange[1][1] = nummodes[1];
1531 
1532  if(edgeOrient==eBackwards)
1533  {
1534  signChange = true;
1535  }
1536  }
1537  }
1538  break;
1539  case 3:
1540  case 11:
1541  {
1542  if( bType[1] == LibUtilities::eGLL_Lagrange)
1543  {
1544  IdxRange[1][0] = 1;
1545  IdxRange[1][1] = nummodes[1] - 1;
1546 
1547  if(edgeOrient==eForwards)
1548  {
1549  reverseOrdering = true;
1550  }
1551  }
1552  else
1553  {
1554  IdxRange[1][0] = 2;
1555  IdxRange[1][1] = nummodes[1];
1556 
1557  if(edgeOrient==eForwards)
1558  {
1559  signChange = true;
1560  }
1561  }
1562  }
1563  break;
1564  }
1565 
1566  switch(eid)
1567  {
1568  case 3:
1569  case 4:
1570  case 7:
1571  case 11:
1572  {
1573  IdxRange[0][0] = 0;
1574  IdxRange[0][1] = 1;
1575  }
1576  break;
1577  case 1:
1578  case 5:
1579  case 6:
1580  case 9:
1581  {
1582  if( bType[0] == LibUtilities::eGLL_Lagrange)
1583  {
1584  IdxRange[0][0] = nummodes[0] - 1;
1585  IdxRange[0][1] = nummodes[0];
1586  }
1587  else
1588  {
1589  IdxRange[0][0] = 1;
1590  IdxRange[0][1] = 2;
1591  }
1592  }
1593  break;
1594  case 0:
1595  case 8:
1596  {
1597  if( bType[0] == LibUtilities::eGLL_Lagrange)
1598  {
1599  IdxRange[0][0] = 1;
1600  IdxRange[0][1] = nummodes[0] - 1;
1601 
1602  if(edgeOrient==eBackwards)
1603  {
1604  reverseOrdering = true;
1605  }
1606  }
1607  else
1608  {
1609  IdxRange[0][0] = 2;
1610  IdxRange[0][1] = nummodes[0];
1611 
1612  if(edgeOrient==eBackwards)
1613  {
1614  signChange = true;
1615  }
1616  }
1617  }
1618  break;
1619  case 2:
1620  case 10:
1621  {
1622  if( bType[0] == LibUtilities::eGLL_Lagrange)
1623  {
1624  IdxRange[0][0] = 1;
1625  IdxRange[0][1] = nummodes[0] - 1;
1626 
1627  if(edgeOrient==eForwards)
1628  {
1629  reverseOrdering = true;
1630  }
1631  }
1632  else
1633  {
1634  IdxRange[0][0] = 2;
1635  IdxRange[0][1] = nummodes[0];
1636 
1637  if(edgeOrient==eForwards)
1638  {
1639  signChange = true;
1640  }
1641  }
1642  }
1643  break;
1644  }
1645 
1646  int p,q,r;
1647  int cnt = 0;
1648 
1649  for(r = IdxRange[2][0]; r < IdxRange[2][1]; r++)
1650  {
1651  for(q = IdxRange[1][0]; q < IdxRange[1][1]; q++)
1652  {
1653  for(p = IdxRange[0][0]; p < IdxRange[0][1]; p++)
1654  {
1655  maparray[cnt++]
1656  = r*nummodes[0]*nummodes[1] + q*nummodes[0] + p;
1657  }
1658  }
1659  }
1660 
1661  if( reverseOrdering )
1662  {
1663  reverse( maparray.get() , maparray.get()+nEdgeIntCoeffs );
1664  }
1665 
1666  if( signChange )
1667  {
1668  for(p = 1; p < nEdgeIntCoeffs; p+=2)
1669  {
1670  signarray[p] = -1;
1671  }
1672  }
1673  }
1674 
1675 
1676  /**
1677  * Generate mapping describing which elemental modes lie on the
1678  * interior of a given face. Accounts for face orientation.
1679  */
1680  void StdHexExp::v_GetFaceInteriorMap(const int fid,
1681  const Orientation faceOrient,
1682  Array<OneD, unsigned int> &maparray,
1683  Array<OneD, int>& signarray)
1684  {
1685  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1686  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1687  "BasisType is not a boundary interior form");
1688  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1689  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1690  "BasisType is not a boundary interior form");
1691  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
1692  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1693  "BasisType is not a boundary interior form");
1694 
1695  ASSERTL1((fid>=0)&&(fid<6),
1696  "local face id must be between 0 and 5");
1697 
1698  int nFaceIntCoeffs = GetFaceIntNcoeffs(fid);
1699 
1700  if(maparray.num_elements()!=nFaceIntCoeffs)
1701  {
1702  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1703  }
1704 
1705  if(signarray.num_elements() != nFaceIntCoeffs)
1706  {
1707  signarray = Array<OneD, int>(nFaceIntCoeffs,1);
1708  }
1709  else
1710  {
1711  fill( signarray.get() , signarray.get()+nFaceIntCoeffs, 1 );
1712  }
1713 
1714  int nummodes [3] = {m_base[0]->GetNumModes(),
1715  m_base[1]->GetNumModes(),
1716  m_base[2]->GetNumModes()};
1717 
1718  const LibUtilities::BasisType bType [3] = {GetBasisType(0),
1719  GetBasisType(1),
1720  GetBasisType(2)};
1721 
1722  int nummodesA = 0;
1723  int nummodesB = 0;
1724 
1725  // Determine the number of modes in face directions A & B based
1726  // on the face index given.
1727  switch(fid)
1728  {
1729  case 0:
1730  case 5:
1731  {
1732  nummodesA = nummodes[0];
1733  nummodesB = nummodes[1];
1734  }
1735  break;
1736  case 1:
1737  case 3:
1738  {
1739  nummodesA = nummodes[0];
1740  nummodesB = nummodes[2];
1741  }
1742  break;
1743  case 2:
1744  case 4:
1745  {
1746  nummodesA = nummodes[1];
1747  nummodesB = nummodes[2];
1748  }
1749  }
1750 
1751  int i,j;
1752  Array<OneD, int> arrayindx(nFaceIntCoeffs);
1753 
1754  // Create a mapping array to account for transposition of the
1755  // coordinates due to face orientation.
1756  for(i = 0; i < (nummodesB-2); i++)
1757  {
1758  for(j = 0; j < (nummodesA-2); j++)
1759  {
1760  if( faceOrient < 9 )
1761  {
1762  arrayindx[i*(nummodesA-2)+j] = i*(nummodesA-2)+j;
1763  }
1764  else
1765  {
1766  arrayindx[i*(nummodesA-2)+j] = j*(nummodesB-2)+i;
1767  }
1768  }
1769  }
1770 
1771  int IdxRange [3][2];
1772  int Incr[3];
1773 
1774  Array<OneD, int> sign0(nummodes[0], 1);
1775  Array<OneD, int> sign1(nummodes[1], 1);
1776  Array<OneD, int> sign2(nummodes[2], 1);
1777 
1778  // Set the upper and lower bounds, and increment for the faces
1779  // involving the first coordinate direction.
1780  switch(fid)
1781  {
1782  case 0: // bottom face
1783  {
1784  IdxRange[2][0] = 0;
1785  IdxRange[2][1] = 1;
1786  Incr[2] = 1;
1787  }
1788  break;
1789  case 5: // top face
1790  {
1791  if( bType[2] == LibUtilities::eGLL_Lagrange)
1792  {
1793  IdxRange[2][0] = nummodes[2] - 1;
1794  IdxRange[2][1] = nummodes[2];
1795  Incr[2] = 1;
1796  }
1797  else
1798  {
1799  IdxRange[2][0] = 1;
1800  IdxRange[2][1] = 2;
1801  Incr[2] = 1;
1802  }
1803 
1804  }
1805  break;
1806  default: // all other faces
1807  {
1808  if( bType[2] == LibUtilities::eGLL_Lagrange)
1809  {
1810  if( (((int) faceOrient)-5) % 2 )
1811  {
1812  IdxRange[2][0] = nummodes[2] - 2;
1813  IdxRange[2][1] = 0;
1814  Incr[2] = -1;
1815 
1816  }
1817  else
1818  {
1819  IdxRange[2][0] = 1;
1820  IdxRange[2][1] = nummodes[2] - 1;
1821  Incr[2] = 1;
1822  }
1823  }
1824  else
1825  {
1826  IdxRange[2][0] = 2;
1827  IdxRange[2][1] = nummodes[2];
1828  Incr[2] = 1;
1829 
1830  if( (((int) faceOrient)-5) % 2 )
1831  {
1832  for(i = 3; i < nummodes[2]; i+=2)
1833  {
1834  sign2[i] = -1;
1835  }
1836  }
1837  }
1838  }
1839  }
1840 
1841  // Set the upper and lower bounds, and increment for the faces
1842  // involving the second coordinate direction.
1843  switch(fid)
1844  {
1845  case 1:
1846  {
1847  IdxRange[1][0] = 0;
1848  IdxRange[1][1] = 1;
1849  Incr[1] = 1;
1850  }
1851  break;
1852  case 3:
1853  {
1854  if( bType[1] == LibUtilities::eGLL_Lagrange)
1855  {
1856  IdxRange[1][0] = nummodes[1] - 1;
1857  IdxRange[1][1] = nummodes[1];
1858  Incr[1] = 1;
1859  }
1860  else
1861  {
1862  IdxRange[1][0] = 1;
1863  IdxRange[1][1] = 2;
1864  Incr[1] = 1;
1865  }
1866  }
1867  break;
1868  case 0:
1869  case 5:
1870  {
1871  if( bType[1] == LibUtilities::eGLL_Lagrange)
1872  {
1873  if( (((int) faceOrient)-5) % 2 )
1874  {
1875  IdxRange[1][0] = nummodes[1] - 2;
1876  IdxRange[1][1] = 0;
1877  Incr[1] = -1;
1878 
1879  }
1880  else
1881  {
1882  IdxRange[1][0] = 1;
1883  IdxRange[1][1] = nummodes[1] - 1;
1884  Incr[1] = 1;
1885  }
1886  }
1887  else
1888  {
1889  IdxRange[1][0] = 2;
1890  IdxRange[1][1] = nummodes[1];
1891  Incr[1] = 1;
1892 
1893  if( (((int) faceOrient)-5) % 2 )
1894  {
1895  for(i = 3; i < nummodes[1]; i+=2)
1896  {
1897  sign1[i] = -1;
1898  }
1899  }
1900  }
1901  }
1902  break;
1903  default: // case2: case4:
1904  {
1905  if( bType[1] == LibUtilities::eGLL_Lagrange)
1906  {
1907  if( (((int) faceOrient)-5) % 4 > 1 )
1908  {
1909  IdxRange[1][0] = nummodes[1] - 2;
1910  IdxRange[1][1] = 0;
1911  Incr[1] = -1;
1912 
1913  }
1914  else
1915  {
1916  IdxRange[1][0] = 1;
1917  IdxRange[1][1] = nummodes[1] - 1;
1918  Incr[1] = 1;
1919  }
1920  }
1921  else
1922  {
1923  IdxRange[1][0] = 2;
1924  IdxRange[1][1] = nummodes[1];
1925  Incr[1] = 1;
1926 
1927  if( (((int) faceOrient)-5) % 4 > 1 )
1928  {
1929  for(i = 3; i < nummodes[1]; i+=2)
1930  {
1931  sign1[i] = -1;
1932  }
1933  }
1934  }
1935  }
1936  }
1937 
1938  switch(fid)
1939  {
1940  case 4:
1941  {
1942  IdxRange[0][0] = 0;
1943  IdxRange[0][1] = 1;
1944  Incr[0] = 1;
1945  }
1946  break;
1947  case 2:
1948  {
1949  if( bType[0] == LibUtilities::eGLL_Lagrange)
1950  {
1951  IdxRange[0][0] = nummodes[0] - 1;
1952  IdxRange[0][1] = nummodes[0];
1953  Incr[0] = 1;
1954  }
1955  else
1956  {
1957  IdxRange[0][0] = 1;
1958  IdxRange[0][1] = 2;
1959  Incr[0] = 1;
1960  }
1961  }
1962  break;
1963  default:
1964  {
1965  if( bType[0] == LibUtilities::eGLL_Lagrange)
1966  {
1967  if( (((int) faceOrient)-5) % 4 > 1 )
1968  {
1969  IdxRange[0][0] = nummodes[0] - 2;
1970  IdxRange[0][1] = 0;
1971  Incr[0] = -1;
1972 
1973  }
1974  else
1975  {
1976  IdxRange[0][0] = 1;
1977  IdxRange[0][1] = nummodes[0] - 1;
1978  Incr[0] = 1;
1979  }
1980  }
1981  else
1982  {
1983  IdxRange[0][0] = 2;
1984  IdxRange[0][1] = nummodes[0];
1985  Incr[0] = 1;
1986 
1987  if( (((int) faceOrient)-5) % 4 > 1 )
1988  {
1989  for(i = 3; i < nummodes[0]; i+=2)
1990  {
1991  sign0[i] = -1;
1992  }
1993  }
1994  }
1995  }
1996  }
1997 
1998  int p,q,r;
1999  int cnt = 0;
2000 
2001  for(r = IdxRange[2][0]; r != IdxRange[2][1]; r+=Incr[2])
2002  {
2003  for(q = IdxRange[1][0]; q != IdxRange[1][1]; q+=Incr[1])
2004  {
2005  for(p = IdxRange[0][0]; p != IdxRange[0][1]; p+=Incr[0])
2006  {
2007  maparray [ arrayindx[cnt ] ]
2008  = r*nummodes[0]*nummodes[1] + q*nummodes[0] + p;
2009  signarray[ arrayindx[cnt++] ]
2010  = sign0[p] * sign1[q] * sign2[r];
2011  }
2012  }
2013  }
2014  }
2015 
2016 
2017  /**
2018  * @param outarray Storage area for computed map.
2019  */
2020  void StdHexExp::v_GetInteriorMap(Array<OneD, unsigned int>& outarray)
2021  {
2022  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
2023  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
2024  "BasisType is not a boundary interior form");
2025  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
2026  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
2027  "BasisType is not a boundary interior form");
2028  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
2029  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
2030  "BasisType is not a boundary interior form");
2031 
2032  int i;
2033  int nummodes [3] = {m_base[0]->GetNumModes(),
2034  m_base[1]->GetNumModes(),
2035  m_base[2]->GetNumModes()};
2036 
2037  int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
2038 
2039  if(outarray.num_elements() != nIntCoeffs)
2040  {
2041  outarray = Array<OneD, unsigned int>(nIntCoeffs);
2042  }
2043 
2044  const LibUtilities::BasisType Btype [3] = {GetBasisType(0),
2045  GetBasisType(1),
2046  GetBasisType(2)};
2047 
2048  int p,q,r;
2049  int cnt = 0;
2050 
2051  int IntIdx [3][2];
2052 
2053  for(i = 0; i < 3; i++)
2054  {
2055  if( Btype[i] == LibUtilities::eModified_A)
2056  {
2057  IntIdx[i][0] = 2;
2058  IntIdx[i][1] = nummodes[i];
2059  }
2060  else
2061  {
2062  IntIdx[i][0] = 1;
2063  IntIdx[i][1] = nummodes[i]-1;
2064  }
2065  }
2066 
2067  for(r = IntIdx[2][0]; r < IntIdx[2][1]; r++)
2068  {
2069  for( q = IntIdx[1][0]; q < IntIdx[1][1]; q++)
2070  {
2071  for( p = IntIdx[0][0]; p < IntIdx[0][1]; p++)
2072  {
2073  outarray[cnt++] = r*nummodes[0]*nummodes[1] +
2074  q*nummodes[0] + p;
2075  }
2076  }
2077  }
2078  }
2079 
2080 
2081  /**
2082  * @param outarray Storage for computed map.
2083  */
2084  void StdHexExp::v_GetBoundaryMap(Array<OneD, unsigned int>& outarray)
2085  {
2086  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
2087  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
2088  "BasisType is not a boundary interior form");
2089  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
2090  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
2091  "BasisType is not a boundary interior form");
2092  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
2093  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
2094  "BasisType is not a boundary interior form");
2095 
2096  int i;
2097  int nummodes [3] = {m_base[0]->GetNumModes(),
2098  m_base[1]->GetNumModes(),
2099  m_base[2]->GetNumModes()};
2100 
2101  int nBndCoeffs = NumBndryCoeffs();
2102 
2103  if(outarray.num_elements()!=nBndCoeffs)
2104  {
2105  outarray = Array<OneD, unsigned int>(nBndCoeffs);
2106  }
2107 
2108  const LibUtilities::BasisType Btype [3] = {GetBasisType(0),
2109  GetBasisType(1),
2110  GetBasisType(2)};
2111 
2112  int p,q,r;
2113  int cnt = 0;
2114 
2115  int BndIdx [3][2];
2116  int IntIdx [3][2];
2117 
2118  for(i = 0; i < 3; i++)
2119  {
2120  BndIdx[i][0] = 0;
2121 
2122  if( Btype[i] == LibUtilities::eModified_A)
2123  {
2124  BndIdx[i][1] = 1;
2125  IntIdx[i][0] = 2;
2126  IntIdx[i][1] = nummodes[i];
2127  }
2128  else
2129  {
2130  BndIdx[i][1] = nummodes[i]-1;
2131  IntIdx[i][0] = 1;
2132  IntIdx[i][1] = nummodes[i]-1;
2133  }
2134  }
2135 
2136 
2137  for(i = 0; i < 2; i++)
2138  {
2139  r = BndIdx[2][i];
2140  for( q = 0; q < nummodes[1]; q++)
2141  {
2142  for( p = 0; p < nummodes[0]; p++)
2143  {
2144  outarray[cnt++] = r*nummodes[0]*nummodes[1]+q*nummodes[0] + p;
2145  }
2146  }
2147  }
2148 
2149  for(r = IntIdx[2][0]; r < IntIdx[2][1]; r++)
2150  {
2151  for( i = 0; i < 2; i++)
2152  {
2153  q = BndIdx[1][i];
2154  for( p = 0; p < nummodes[0]; p++)
2155  {
2156  outarray[cnt++] = r*nummodes[0]*nummodes[1] +
2157  q*nummodes[0] + p;
2158  }
2159  }
2160 
2161  for( q = IntIdx[1][0]; q < IntIdx[1][1]; q++)
2162  {
2163  for( i = 0; i < 2; i++)
2164  {
2165  p = BndIdx[0][i];
2166  outarray[cnt++] = r*nummodes[0]*nummodes[1] +
2167  q*nummodes[0] + p;
2168  }
2169  }
2170  }
2171 
2172  sort(outarray.get(), outarray.get() + nBndCoeffs);
2173  }
2174 
2175  DNekMatSharedPtr StdHexExp::v_GenMatrix(const StdMatrixKey &mkey)
2176  {
2177  return StdExpansion::CreateGeneralMatrix(mkey);
2178  }
2179 
2180 
2181  DNekMatSharedPtr StdHexExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
2182  {
2183  return StdExpansion::CreateGeneralMatrix(mkey);
2184  }
2185 
2186 
2187  void StdHexExp::v_MassMatrixOp(
2188  const Array<OneD, const NekDouble> &inarray,
2189  Array<OneD,NekDouble> &outarray,
2190  const StdMatrixKey &mkey)
2191  {
2192  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
2193  }
2194 
2195 
2196  void StdHexExp::v_LaplacianMatrixOp(
2197  const Array<OneD, const NekDouble> &inarray,
2198  Array<OneD,NekDouble> &outarray,
2199  const StdMatrixKey &mkey)
2200  {
2201  StdHexExp::v_LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
2202  }
2203 
2204 
2205  void StdHexExp::v_LaplacianMatrixOp(const int k1, const int k2,
2206  const Array<OneD, const NekDouble> &inarray,
2207  Array<OneD,NekDouble> &outarray,
2208  const StdMatrixKey &mkey)
2209  {
2210  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,
2211  mkey);
2212  }
2213 
2214 
2215  void StdHexExp::v_WeakDerivMatrixOp(const int i,
2216  const Array<OneD, const NekDouble> &inarray,
2217  Array<OneD,NekDouble> &outarray,
2218  const StdMatrixKey &mkey)
2219  {
2220  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,
2221  mkey);
2222  }
2223 
2224  void StdHexExp::v_HelmholtzMatrixOp(
2225  const Array<OneD, const NekDouble> &inarray,
2226  Array<OneD,NekDouble> &outarray,
2227  const StdMatrixKey &mkey)
2228  {
2229  StdHexExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
2230  }
2231 
2232 
2233  void StdHexExp::v_GeneralMatrixOp_MatOp(
2234  const Array<OneD, const NekDouble> &inarray,
2235  Array<OneD,NekDouble> &outarray,
2236  const StdMatrixKey &mkey)
2237  {
2238  DNekMatSharedPtr mat = m_stdMatrixManager[mkey];
2239 
2240  if(inarray.get() == outarray.get())
2241  {
2242  Array<OneD,NekDouble> tmp(m_ncoeffs);
2243  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
2244 
2245  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
2246  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
2247  }
2248  else
2249  {
2250  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
2251  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
2252  }
2253  }
2254 
2255 
2256  void StdHexExp::v_MultiplyByStdQuadratureMetric(const Array<OneD, const NekDouble>& inarray,
2257  Array<OneD, NekDouble> &outarray)
2258  {
2259  int i;
2260  int nquad0 = m_base[0]->GetNumPoints();
2261  int nquad1 = m_base[1]->GetNumPoints();
2262  int nquad2 = m_base[2]->GetNumPoints();
2263  int nq01 = nquad0*nquad1;
2264  int nq12 = nquad1*nquad2;
2265 
2266  const Array<OneD, const NekDouble>& w0 = m_base[0]->GetW();
2267  const Array<OneD, const NekDouble>& w1 = m_base[1]->GetW();
2268  const Array<OneD, const NekDouble>& w2 = m_base[2]->GetW();
2269 
2270  for(i = 0; i < nq12; ++i)
2271  {
2272  Vmath::Vmul(nquad0, inarray.get()+i*nquad0, 1,
2273  w0.get(), 1, outarray.get()+i*nquad0,1);
2274  }
2275 
2276  for(i = 0; i < nq12; ++i)
2277  {
2278  Vmath::Smul(nquad0, w1[i%nquad1], outarray.get()+i*nquad0, 1,
2279  outarray.get()+i*nquad0, 1);
2280  }
2281 
2282  for(i = 0; i < nquad2; ++i)
2283  {
2284  Vmath::Smul(nq01, w2[i], outarray.get()+i*nq01, 1,
2285  outarray.get()+i*nq01, 1);
2286  }
2287  }
2288 
2289  void StdHexExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
2290  const StdMatrixKey &mkey)
2291  {
2292  // Generate an orthonogal expansion
2293  int qa = m_base[0]->GetNumPoints();
2294  int qb = m_base[1]->GetNumPoints();
2295  int qc = m_base[2]->GetNumPoints();
2296  int nmodes_a = m_base[0]->GetNumModes();
2297  int nmodes_b = m_base[1]->GetNumModes();
2298  int nmodes_c = m_base[2]->GetNumModes();
2299  // Declare orthogonal basis.
2300  LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
2301  LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
2302  LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
2303 
2304  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A,nmodes_a,pa);
2305  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A,nmodes_b,pb);
2306  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_A,nmodes_c,pc);
2307  StdHexExp OrthoExp(Ba,Bb,Bc);
2308 
2309  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2310  int i,j,k;
2311 
2312  int cutoff = (int) (mkey.GetConstFactor(eFactorSVVCutoffRatio)*min(nmodes_a,nmodes_b));
2313  NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
2314 
2315  // project onto modal space.
2316  OrthoExp.FwdTrans(array,orthocoeffs);
2317 
2318 
2319  // Filter just trilinear space
2320  int nmodes = max(nmodes_a,nmodes_b);
2321  nmodes = max(nmodes,nmodes_c);
2322 
2323  Array<OneD, NekDouble> fac(nmodes,1.0);
2324  for(j = cutoff; j < nmodes; ++j)
2325  {
2326  fac[j] = fabs((j-nmodes)/((NekDouble) (j-cutoff+1.0)));
2327  fac[j] *= fac[j]; //added this line to conform with equation
2328  }
2329 
2330  for(i = 0; i < nmodes_a; ++i)
2331  {
2332  for(j = 0; j < nmodes_b; ++j)
2333  {
2334  for(k = 0; k < nmodes_c; ++k)
2335  {
2336  if((i >= cutoff)||(j >= cutoff)||(k >= cutoff))
2337  {
2338  orthocoeffs[i*nmodes_a*nmodes_b + j*nmodes_c + k] *= (SvvDiffCoeff*exp( -(fac[i]+fac[j]+fac[k]) ));
2339  }
2340  else
2341  {
2342  orthocoeffs[i*nmodes_a*nmodes_b + j*nmodes_c + k] *= 0.0;
2343  }
2344  }
2345  }
2346  }
2347 
2348  // backward transform to physical space
2349  OrthoExp.BwdTrans(orthocoeffs,array);
2350  }
2351  }
2352 }
2353 
Nektar::StdRegions::StdHexExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdHexExp.cpp:2020
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:454
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:2175
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:135
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:2224
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:1401
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:902
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:613
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:1423
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:2084
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:2187
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:2215
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:2233
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:1014
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:684
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:397
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:2256
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:1382
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:344
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:2196
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:743
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:43
Nektar::StdRegions::StdHexExp::v_BwdTrans_SumFac
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:212
Nektar::StdRegions::StdHexExp::v_GetVertexMap
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdHexExp.cpp:1248
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:2289
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:187
Nektar::StdRegions::StdExpansion::MassMatrixOp_MatFree
void MassMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.cpp:627
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:509
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:1359
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:1680
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:165
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1356
Nektar::StdRegions::StdHexExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2181
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1016
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:1394
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:1823
Nektar::StdRegions::StdHexExp::v_NumDGBndryCoeffs
virtual int v_NumDGBndryCoeffs() const
Definition: StdHexExp.cpp:718
Generated on Tue Dec 8 2015 23:41:30 for Nektar++ by   doxygen 1.8.9.1