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