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