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