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