Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HexExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File HexExp.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: Methods for Hex expansion in local regoins
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 
37 #include <LocalRegions/HexExp.h>
38 #include <LibUtilities/Foundations/Interp.h>
39 #include <LibUtilities/Foundations/InterpCoeff.h>
40 #include <SpatialDomains/HexGeom.h>
41 
42 namespace Nektar
43 {
44  namespace LocalRegions
45  {
46  /**
47  * @class HexExp
48  * Defines a hexahedral local expansion.
49  */
50 
51  /**
52  * \brief Constructor using BasisKey class for quadrature points and
53  * order definition
54  *
55  * @param Ba Basis key for first coordinate.
56  * @param Bb Basis key for second coordinate.
57  * @param Bc Basis key for third coordinate.
58  */
59  HexExp::HexExp(const LibUtilities::BasisKey &Ba,
60  const LibUtilities::BasisKey &Bb,
61  const LibUtilities::BasisKey &Bc,
62  const SpatialDomains::HexGeomSharedPtr &geom):
63  StdExpansion (Ba.GetNumModes()*Bb.GetNumModes()*Bc.GetNumModes(),3,Ba,Bb,Bc),
64  StdExpansion3D(Ba.GetNumModes()*Bb.GetNumModes()*Bc.GetNumModes(),Ba,Bb,Bc),
65  StdRegions::StdHexExp(Ba,Bb,Bc),
66  Expansion (geom),
67  Expansion3D (geom),
68  m_matrixManager(
69  boost::bind(&HexExp::CreateMatrix, this, _1),
70  std::string("HexExpMatrix")),
71  m_staticCondMatrixManager(
72  boost::bind(&HexExp::CreateStaticCondMatrix, this, _1),
73  std::string("HexExpStaticCondMatrix"))
74  {
75  }
76 
77 
78  /**
79  * \brief Copy Constructor
80  *
81  * @param T HexExp to copy.
82  */
83  HexExp::HexExp(const HexExp &T):
84  StdExpansion(T),
85  StdExpansion3D(T),
86  StdRegions::StdHexExp(T),
87  Expansion(T),
88  Expansion3D(T),
89  m_matrixManager(T.m_matrixManager),
90  m_staticCondMatrixManager(T.m_staticCondMatrixManager)
91  {
92  }
93 
94  /**
95  * \brief Destructor
96  */
97  HexExp::~HexExp()
98  {
99  }
100 
101 
102  //-----------------------------
103  // Integration Methods
104  //-----------------------------
105  /**
106  * \brief Integrate the physical point list \a inarray over region
107  *
108  * @param inarray definition of function to be returned at
109  * quadrature points of expansion.
110  * @returns \f$\int^1_{-1}\int^1_{-1} \int^1_{-1}
111  * u(\eta_1, \eta_2, \eta_3) J[i,j,k] d \eta_1 d \eta_2 d \eta_3 \f$
112  * where \f$inarray[i,j,k] = u(\eta_{1i},\eta_{2j},\eta_{3k}) \f$
113  * and \f$ J[i,j,k] \f$ is the Jacobian evaluated at the quadrature
114  * point.
115  */
116  NekDouble HexExp::v_Integral(
117  const Array<OneD, const NekDouble> &inarray)
118  {
119  int nquad0 = m_base[0]->GetNumPoints();
120  int nquad1 = m_base[1]->GetNumPoints();
121  int nquad2 = m_base[2]->GetNumPoints();
122  Array<OneD, const NekDouble> jac = m_metricinfo->GetJac(GetPointsKeys());
123  NekDouble returnVal;
124  Array<OneD,NekDouble> tmp(nquad0*nquad1*nquad2);
125 
126  // multiply inarray with Jacobian
127 
128  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
129  {
130  Vmath::Vmul(nquad0*nquad1*nquad2,&jac[0],1,
131  (NekDouble*)&inarray[0],1,&tmp[0],1);
132  }
133  else
134  {
135  Vmath::Smul(nquad0*nquad1*nquad2,(NekDouble) jac[0],
136  (NekDouble*)&inarray[0],1,&tmp[0],1);
137  }
138 
139  // call StdHexExp version;
140  returnVal = StdHexExp::v_Integral(tmp);
141 
142  return returnVal;
143  }
144 
145 
146  //-----------------------------
147  // Differentiation Methods
148  //-----------------------------
149  /**
150  * \brief Calculate the derivative of the physical points
151  *
152  * For Hexahedral region can use the Tensor_Deriv function defined
153  * under StdExpansion.
154  * @param inarray Input array
155  * @param out_d0 Derivative of \a inarray in first direction.
156  * @param out_d1 Derivative of \a inarray in second direction.
157  * @param out_d2 Derivative of \a inarray in third direction.
158  */
159  void HexExp::v_PhysDeriv(
160  const Array<OneD, const NekDouble> & inarray,
161  Array<OneD,NekDouble> &out_d0,
162  Array<OneD,NekDouble> &out_d1,
163  Array<OneD,NekDouble> &out_d2)
164  {
165  int nquad0 = m_base[0]->GetNumPoints();
166  int nquad1 = m_base[1]->GetNumPoints();
167  int nquad2 = m_base[2]->GetNumPoints();
168  int ntot = nquad0 * nquad1 * nquad2;
169 
170  Array<TwoD, const NekDouble> df =
171  m_metricinfo->GetDerivFactors(GetPointsKeys());
172  Array<OneD,NekDouble> Diff0 = Array<OneD,NekDouble>(ntot);
173  Array<OneD,NekDouble> Diff1 = Array<OneD,NekDouble>(ntot);
174  Array<OneD,NekDouble> Diff2 = Array<OneD,NekDouble>(ntot);
175 
176  StdHexExp::v_PhysDeriv(inarray, Diff0, Diff1, Diff2);
177 
178  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
179  {
180  if(out_d0.num_elements())
181  {
182  Vmath::Vmul (ntot,&df[0][0],1,&Diff0[0],1, &out_d0[0], 1);
183  Vmath::Vvtvp(ntot,&df[1][0],1,&Diff1[0],1, &out_d0[0], 1,
184  &out_d0[0],1);
185  Vmath::Vvtvp(ntot,&df[2][0],1,&Diff2[0],1, &out_d0[0], 1,
186  &out_d0[0],1);
187  }
188 
189  if(out_d1.num_elements())
190  {
191  Vmath::Vmul (ntot,&df[3][0],1,&Diff0[0],1, &out_d1[0], 1);
192  Vmath::Vvtvp(ntot,&df[4][0],1,&Diff1[0],1, &out_d1[0], 1,
193  &out_d1[0],1);
194  Vmath::Vvtvp(ntot,&df[5][0],1,&Diff2[0],1, &out_d1[0], 1,
195  &out_d1[0],1);
196  }
197 
198  if(out_d2.num_elements())
199  {
200  Vmath::Vmul (ntot,&df[6][0],1,&Diff0[0],1, &out_d2[0], 1);
201  Vmath::Vvtvp(ntot,&df[7][0],1,&Diff1[0],1, &out_d2[0], 1,
202  &out_d2[0],1);
203  Vmath::Vvtvp(ntot,&df[8][0],1,&Diff2[0],1, &out_d2[0], 1,
204  &out_d2[0],1);
205  }
206  }
207  else // regular geometry
208  {
209  if(out_d0.num_elements())
210  {
211  Vmath::Smul (ntot,df[0][0],&Diff0[0],1, &out_d0[0], 1);
212  Blas::Daxpy (ntot,df[1][0],&Diff1[0],1, &out_d0[0], 1);
213  Blas::Daxpy (ntot,df[2][0],&Diff2[0],1, &out_d0[0], 1);
214  }
215 
216  if(out_d1.num_elements())
217  {
218  Vmath::Smul (ntot,df[3][0],&Diff0[0],1, &out_d1[0], 1);
219  Blas::Daxpy (ntot,df[4][0],&Diff1[0],1, &out_d1[0], 1);
220  Blas::Daxpy (ntot,df[5][0],&Diff2[0],1, &out_d1[0], 1);
221  }
222 
223  if(out_d2.num_elements())
224  {
225  Vmath::Smul (ntot,df[6][0],&Diff0[0],1, &out_d2[0], 1);
226  Blas::Daxpy (ntot,df[7][0],&Diff1[0],1, &out_d2[0], 1);
227  Blas::Daxpy (ntot,df[8][0],&Diff2[0],1, &out_d2[0], 1);
228  }
229  }
230  }
231 
232 
233  /**
234  * \brief Calculate the derivative of the physical points in a single
235  * direction.
236  *
237  * @param dir Direction in which to compute derivative.
238  * Valid values are 0, 1, 2.
239  * @param inarray Input array.
240  * @param outarray Output array.
241  */
242  void HexExp::v_PhysDeriv(
243  const int dir,
244  const Array<OneD, const NekDouble>& inarray,
245  Array<OneD, NekDouble>& outarray)
246  {
247  switch(dir)
248  {
249  case 0:
250  {
251  PhysDeriv(inarray, outarray, NullNekDouble1DArray,
252  NullNekDouble1DArray);
253  }
254  break;
255  case 1:
256  {
257  PhysDeriv(inarray, NullNekDouble1DArray, outarray,
258  NullNekDouble1DArray);
259  }
260  break;
261  case 2:
262  {
263  PhysDeriv(inarray, NullNekDouble1DArray,
264  NullNekDouble1DArray, outarray);
265  }
266  break;
267  default:
268  {
269  ASSERTL1(false,"input dir is out of range");
270  }
271  break;
272  }
273  }
274 
275 
276  //-----------------------------
277  // Transforms
278  //-----------------------------
279 
280  /**
281  * \brief Forward transform from physical quadrature space stored in \a
282  * inarray and evaluate the expansion coefficients and store in
283  * \a (this)->_coeffs
284  *
285  * @param inarray Input array
286  * @param outarray Output array
287  */
288  void HexExp::v_FwdTrans(
289  const Array<OneD, const NekDouble> & inarray,
290  Array<OneD,NekDouble> &outarray)
291  {
292  if( m_base[0]->Collocation() && m_base[1]->Collocation()
293  && m_base[2]->Collocation())
294  {
295  Vmath::Vcopy(GetNcoeffs(),&inarray[0],1,&outarray[0],1);
296  }
297  else
298  {
299  IProductWRTBase(inarray,outarray);
300 
301  // get Mass matrix inverse
302  MatrixKey masskey(StdRegions::eInvMass,
303  DetShapeType(),*this);
304  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
305 
306  // copy inarray in case inarray == outarray
307  DNekVec in (m_ncoeffs,outarray);
308  DNekVec out(m_ncoeffs,outarray,eWrapper);
309 
310  out = (*matsys)*in;
311  }
312  }
313 
314 
315  //-----------------------------
316  // Inner product functions
317  //-----------------------------
318 
319  /**
320  * \brief Calculate the inner product of inarray with respect to the
321  * elements basis.
322  *
323  * @param inarray Input array of physical space data.
324  * @param outarray Output array of data.
325  */
326  void HexExp::v_IProductWRTBase(
327  const Array<OneD, const NekDouble> &inarray,
328  Array<OneD, NekDouble> &outarray)
329  {
330  HexExp::v_IProductWRTBase_SumFac(inarray, outarray);
331  }
332 
333  /**
334  * \brief Calculate the inner product of inarray with respect to the
335  * given basis B = base0 * base1 * base2.
336  *
337  * \f$ \begin{array}{rcl} 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}
340  * (\xi_{3k}) w_i w_j w_k u(\xi_{1,i} \xi_{2,j} \xi_{3,k})
341  * J_{i,j,k}\\ & = & \sum_{i=0}^{nq_0} \psi_p^a(\xi_{1,i})
342  * \sum_{j=0}^{nq_1} \psi_{q}^a(\xi_{2,j}) \sum_{k=0}^{nq_2}
343  * \psi_{r}^a u(\xi_{1i},\xi_{2j},\xi_{3k})
344  * J_{i,j,k} \end{array} \f$ \n
345  * where
346  * \f$ \phi_{pqr} (\xi_1 , \xi_2 , \xi_3)
347  * = \psi_p^a ( \xi_1) \psi_{q}^a (\xi_2) \psi_{r}^a (\xi_3) \f$ \n
348  * which can be implemented as \n
349  * \f$f_{r} (\xi_{3k})
350  * = \sum_{k=0}^{nq_3} \psi_{r}^a u(\xi_{1i},\xi_{2j},\xi_{3k})
351  * J_{i,j,k} = {\bf B_3 U} \f$ \n
352  * \f$ g_{q} (\xi_{3k}) = \sum_{j=0}^{nq_1} \psi_{q}^a (\xi_{2j})
353  * f_{r} (\xi_{3k}) = {\bf B_2 F} \f$ \n
354  * \f$ (\phi_{pqr}, u)_{\delta}
355  * = \sum_{k=0}^{nq_0} \psi_{p}^a (\xi_{3k}) g_{q} (\xi_{3k})
356  * = {\bf B_1 G} \f$
357  *
358  * @param base0 Basis to integrate wrt in first dimension.
359  * @param base1 Basis to integrate wrt in second dimension.
360  * @param base2 Basis to integrate wrt in third dimension.
361  * @param inarray Input array.
362  * @param outarray Output array.
363  * @param coll_check (not used)
364  */
365  void HexExp::v_IProductWRTBase_SumFac(
366  const Array<OneD, const NekDouble> &inarray,
367  Array<OneD, NekDouble> &outarray)
368  {
369  int nquad0 = m_base[0]->GetNumPoints();
370  int nquad1 = m_base[1]->GetNumPoints();
371  int nquad2 = m_base[2]->GetNumPoints();
372  int order0 = m_base[0]->GetNumModes();
373  int order1 = m_base[1]->GetNumModes();
374 
375  Array<OneD, NekDouble> tmp(inarray.num_elements());
376  Array<OneD, NekDouble> wsp(nquad0*nquad1*(nquad2+order0) +
377  order0*order1*nquad2);
378 
379  MultiplyByQuadratureMetric(inarray, tmp);
380  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
381  m_base[1]->GetBdata(),
382  m_base[2]->GetBdata(),
383  tmp,outarray,wsp,
384  true,true,true);
385  }
386 
387  void HexExp::v_IProductWRTDerivBase(
388  const int dir,
389  const Array<OneD, const NekDouble>& inarray,
390  Array<OneD, NekDouble> & outarray)
391  {
392  HexExp::IProductWRTDerivBase_SumFac(dir,inarray,outarray);
393  }
394 
395 
396  /**
397  * @brief Calculates the inner product \f$ I_{pqr} = (u,
398  * \partial_{x_i} \phi_{pqr}) \f$.
399  *
400  * The derivative of the basis functions is performed using the chain
401  * rule in order to incorporate the geometric factors. Assuming that
402  * the basis functions are a tensor product
403  * \f$\phi_{pqr}(\xi_1,\xi_2,\xi_3) =
404  * \phi_1(\xi_1)\phi_2(\xi_2)\phi_3(\xi_3)\f$, in the hexahedral
405  * element, this is straightforward and yields the result
406  *
407  * \f[
408  * I_{pqr} = \sum_{k=1}^3 \left(u, \frac{\partial u}{\partial \xi_k}
409  * \frac{\partial \xi_k}{\partial x_i}\right)
410  * \f]
411  *
412  * @param dir Direction in which to take the derivative.
413  * @param inarray The function \f$ u \f$.
414  * @param outarray Value of the inner product.
415  */
416  void HexExp::IProductWRTDerivBase_SumFac(
417  const int dir,
418  const Array<OneD, const NekDouble>& inarray,
419  Array<OneD, NekDouble> & outarray)
420  {
421  ASSERTL1((dir==0)||(dir==1)||(dir==2),"Invalid direction.");
422 
423  const int nq0 = m_base[0]->GetNumPoints();
424  const int nq1 = m_base[1]->GetNumPoints();
425  const int nq2 = m_base[2]->GetNumPoints();
426  const int nq = nq0*nq1*nq2;
427  const int nm0 = m_base[0]->GetNumModes();
428  const int nm1 = m_base[1]->GetNumModes();
429 
430  const Array<TwoD, const NekDouble>& df =
431  m_metricinfo->GetDerivFactors(GetPointsKeys());
432 
433  Array<OneD, NekDouble> alloc(4*nq + m_ncoeffs + nm0*nq2*(nq1+nm1));
434  Array<OneD, NekDouble> tmp1 (alloc); // Quad metric
435  Array<OneD, NekDouble> tmp2 (alloc + nq); // Dir1 metric
436  Array<OneD, NekDouble> tmp3 (alloc + 2*nq); // Dir2 metric
437  Array<OneD, NekDouble> tmp4 (alloc + 3*nq); // Dir3 metric
438  Array<OneD, NekDouble> tmp5 (alloc + 4*nq); // iprod tmp
439  Array<OneD, NekDouble> wsp (tmp5 + m_ncoeffs); // Wsp
440 
441  MultiplyByQuadratureMetric(inarray, tmp1);
442 
443  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
444  {
445  Vmath::Vmul(nq,&df[3*dir][0], 1,tmp1.get(),1,tmp2.get(),1);
446  Vmath::Vmul(nq,&df[3*dir+1][0],1,tmp1.get(),1,tmp3.get(),1);
447  Vmath::Vmul(nq,&df[3*dir+2][0],1,tmp1.get(),1,tmp4.get(),1);
448  }
449  else
450  {
451  Vmath::Smul(nq, df[3*dir][0], tmp1.get(),1,tmp2.get(), 1);
452  Vmath::Smul(nq, df[3*dir+1][0],tmp1.get(),1,tmp3.get(), 1);
453  Vmath::Smul(nq, df[3*dir+2][0],tmp1.get(),1,tmp4.get(), 1);
454  }
455 
456  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
457  m_base[1]->GetBdata(),
458  m_base[2]->GetBdata(),
459  tmp2,outarray,wsp,
460  false,true,true);
461 
462  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
463  m_base[1]->GetDbdata(),
464  m_base[2]->GetBdata(),
465  tmp3,tmp5,wsp,
466  true,false,true);
467  Vmath::Vadd(m_ncoeffs, tmp5, 1, outarray, 1, outarray, 1);
468 
469  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
470  m_base[1]->GetBdata(),
471  m_base[2]->GetDbdata(),
472  tmp4,tmp5,wsp,
473  true,true,false);
474  Vmath::Vadd(m_ncoeffs, tmp5, 1, outarray, 1, outarray, 1);
475  }
476 
477 
478  void HexExp::IProductWRTDerivBase_MatOp(
479  const int dir,
480  const Array<OneD, const NekDouble>& inarray,
481  Array<OneD, NekDouble> &outarray)
482  {
483  int nq = GetTotPoints();
484  StdRegions::MatrixType mtype = StdRegions::eIProductWRTDerivBase0;
485 
486  switch(dir)
487  {
488  case 0:
489  {
490  mtype = StdRegions::eIProductWRTDerivBase0;
491  }
492  break;
493  case 1:
494  {
495  mtype = StdRegions::eIProductWRTDerivBase1;
496  }
497  break;
498  case 2:
499  {
500  mtype = StdRegions::eIProductWRTDerivBase2;
501  }
502  break;
503  default:
504  {
505  ASSERTL1(false,"input dir is out of range");
506  }
507  break;
508  }
509 
510  MatrixKey iprodmatkey(mtype,DetShapeType(),*this);
511  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
512 
513  Blas::Dgemv('N',m_ncoeffs,nq,iprodmat->Scale(),(iprodmat->GetOwnedMatrix())->GetPtr().get(),
514  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
515  }
516 
517 
518  //-----------------------------
519  // Evaluation functions
520  //-----------------------------
521 
522 
523  /**
524  * Given the local cartesian coordinate \a Lcoord evaluate the
525  * value of physvals at this point by calling through to the
526  * StdExpansion method
527  */
528  NekDouble HexExp::v_StdPhysEvaluate(
529  const Array<OneD, const NekDouble> &Lcoord,
530  const Array<OneD, const NekDouble> &physvals)
531  {
532  // Evaluate point in local coordinates.
533  return StdHexExp::v_PhysEvaluate(Lcoord,physvals);
534  }
535 
536  NekDouble HexExp::v_PhysEvaluate(
537  const Array<OneD, const NekDouble> &coord,
538  const Array<OneD, const NekDouble> & physvals)
539  {
540  Array<OneD,NekDouble> Lcoord = Array<OneD,NekDouble>(3);
541 
542  ASSERTL0(m_geom,"m_geom not defined");
543  m_geom->GetLocCoords(coord,Lcoord);
544  return StdHexExp::v_PhysEvaluate(Lcoord, physvals);
545  }
546 
547  /**
548  * \brief Retrieves the physical coordinates of a given set of
549  * reference coordinates.
550  *
551  * @param Lcoords Local coordinates in reference space.
552  * @param coords Corresponding coordinates in physical space.
553  */
554  void HexExp::v_GetCoord(
555  const Array<OneD, const NekDouble> &Lcoords,
556  Array<OneD,NekDouble> &coords)
557  {
558  int i;
559 
560  ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[0] <= 1.0 &&
561  Lcoords[1] >= -1.0 && Lcoords[1] <= 1.0 &&
562  Lcoords[2] >= -1.0 && Lcoords[2] <= 1.0,
563  "Local coordinates are not in region [-1,1]");
564 
565  m_geom->FillGeom();
566 
567  for(i = 0; i < m_geom->GetCoordim(); ++i)
568  {
569  coords[i] = m_geom->GetCoord(i,Lcoords);
570  }
571  }
572 
573  void HexExp::v_GetCoords(
574  Array<OneD, NekDouble> &coords_0,
575  Array<OneD, NekDouble> &coords_1,
576  Array<OneD, NekDouble> &coords_2)
577  {
578  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
579  }
580 
581  //-----------------------------
582  // Helper functions
583  //-----------------------------
584 
585  /// Return the region shape using the enum-list of ShapeType
586  LibUtilities::ShapeType HexExp::v_DetShapeType() const
587  {
588  return LibUtilities::eHexahedron;
589  }
590 
591 
592  void HexExp::v_ExtractDataToCoeffs(
593  const NekDouble *data,
594  const std::vector<unsigned int > &nummodes,
595  const int mode_offset,
596  NekDouble * coeffs)
597  {
598  int data_order0 = nummodes[mode_offset];
599  int fillorder0 = min(m_base[0]->GetNumModes(),data_order0);
600  int data_order1 = nummodes[mode_offset+1];
601  int order1 = m_base[1]->GetNumModes();
602  int fillorder1 = min(order1,data_order1);
603  int data_order2 = nummodes[mode_offset+2];
604  int order2 = m_base[2]->GetNumModes();
605  int fillorder2 = min(order2,data_order2);
606 
607  switch(m_base[0]->GetBasisType())
608  {
609  case LibUtilities::eModified_A:
610  {
611  int i,j;
612  int cnt = 0;
613  int cnt1 = 0;
614 
615  ASSERTL1(m_base[1]->GetBasisType() ==
616  LibUtilities::eModified_A,
617  "Extraction routine not set up for this basis");
618  ASSERTL1(m_base[2]->GetBasisType() ==
619  LibUtilities::eModified_A,
620  "Extraction routine not set up for this basis");
621 
622  Vmath::Zero(m_ncoeffs,coeffs,1);
623  for(j = 0; j < fillorder0; ++j)
624  {
625  for(i = 0; i < fillorder1; ++i)
626  {
627  Vmath::Vcopy(fillorder2, &data[cnt], 1,
628  &coeffs[cnt1], 1);
629  cnt += data_order2;
630  cnt1 += order2;
631  }
632 
633  // count out data for j iteration
634  for(i = fillorder1; i < data_order1; ++i)
635  {
636  cnt += data_order2;
637  }
638 
639  for(i = fillorder1; i < order1; ++i)
640  {
641  cnt1 += order2;
642  }
643  }
644  }
645  break;
646  default:
647  ASSERTL0(false, "basis is either not set up or not "
648  "hierarchicial");
649  }
650  }
651 
652  StdRegions::Orientation HexExp::v_GetFaceOrient(int face)
653  {
654  return GetGeom3D()->GetFaceOrient(face);
655  }
656 
657 
658  bool HexExp::v_GetFaceDGForwards(const int i) const
659  {
660  StdRegions::Orientation fo = GetGeom3D()->GetFaceOrient(i);
661 
662  return fo == StdRegions::eDir1FwdDir1_Dir2FwdDir2 ||
663  fo == StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
664  fo == StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
665  fo == StdRegions::eDir1FwdDir2_Dir2BwdDir1;
666  }
667 
668 
669  void HexExp::v_GetTracePhysVals(
670  const int face,
671  const StdRegions::StdExpansionSharedPtr &FaceExp,
672  const Array<OneD, const NekDouble> &inarray,
673  Array<OneD, NekDouble> &outarray,
674  StdRegions::Orientation orient)
675  {
676  v_GetFacePhysVals(face,FaceExp,inarray,outarray,orient);
677  }
678 
679 
680  ///Returns the physical values at the quadrature points of a face
681  void HexExp::v_GetFacePhysVals(
682  const int face,
683  const StdRegions::StdExpansionSharedPtr &FaceExp,
684  const Array<OneD, const NekDouble> &inarray,
685  Array<OneD, NekDouble> &outarray,
686  StdRegions::Orientation orient)
687  {
688  int nquad0 = m_base[0]->GetNumPoints();
689  int nquad1 = m_base[1]->GetNumPoints();
690  int nquad2 = m_base[2]->GetNumPoints();
691  Array<OneD, NekDouble> o_tmp(nquad0*nquad1*nquad2);
692 
693  if (orient == StdRegions::eNoOrientation)
694  {
695  orient = GetFaceOrient(face);
696  }
697 
698  switch(face)
699  {
700  case 0:
701  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
702  {
703  //Directions A and B positive
704  Vmath::Vcopy(nquad0*nquad1,&(inarray[0]),1,&(outarray[0]),1);
705  }
706  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
707  {
708  //Direction A negative and B positive
709  for (int j=0; j<nquad1; j++)
710  {
711  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+j*nquad0,-1,&(outarray[0])+(j*nquad0),1);
712  }
713  }
714  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
715  {
716  //Direction A positive and B negative
717  for (int j=0; j<nquad1; j++)
718  {
719  Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1-j),1,&(outarray[0])+(j*nquad0),1);
720  }
721  }
722  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
723  {
724  //Direction A negative and B negative
725  for(int j=0; j<nquad1; j++)
726  {
727  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1-j*nquad0),-1,&(outarray[0])+(j*nquad0),1);
728  }
729  }
730  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
731  {
732  //Transposed, Direction A and B positive
733  for (int i=0; i<nquad0; i++)
734  {
735  Vmath::Vcopy(nquad1,&(inarray[0])+i,nquad0,&(outarray[0])+(i*nquad1),1);
736  }
737  }
738  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
739  {
740  //Transposed, Direction A negative and B positive
741  for (int i=0; i<nquad0; i++)
742  {
743  Vmath::Vcopy(nquad1,&(inarray[0])+i+nquad0*(nquad1-1),-nquad0,&(outarray[0])+(i*nquad1),1);
744  }
745  }
746  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
747  {
748  //Transposed, Direction A positive and B negative
749  for (int i=0; i<nquad0; i++)
750  {
751  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0-1-i),nquad0,&(outarray[0])+(i*nquad1),1);
752  }
753  }
754  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
755  {
756  //Transposed, Direction A and B negative
757  for (int i=0; i<nquad0; i++)
758  {
759  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0*nquad1-1-i),-nquad0,&(outarray[0])+(i*nquad1),1);
760  }
761  }
762  o_tmp = outarray;
763  //interpolate
764  LibUtilities::Interp2D(m_base[0]->GetPointsKey(),
765  m_base[1]->GetPointsKey(), o_tmp,
766  FaceExp->GetBasis(0)->GetPointsKey(),
767  FaceExp->GetBasis(1)->GetPointsKey(),
768  outarray);
769  break;
770  case 1:
771  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
772  {
773  //Direction A and B positive
774  for (int k=0; k<nquad2; k++)
775  {
776  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1*k),
777  1,&(outarray[0])+(k*nquad0),1);
778  }
779  }
780  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
781  {
782  //Direction A negative and B positive
783  for (int k=0; k<nquad2; k++)
784  {
785  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+(nquad0*nquad1*k),
786  -1,&(outarray[0])+(k*nquad0),1);
787  }
788  }
789  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
790  {
791  //Direction A positive and B negative
792  for (int k=0; k<nquad2; k++)
793  {
794  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1*(nquad2-1-k)),
795  1,&(outarray[0])+(k*nquad0),1);
796  }
797  }
798  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
799  {
800  //Direction A negative and B negative
801  for(int k=0; k<nquad2; k++)
802  {
803  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+(nquad0*nquad1*(nquad2-1-k)),
804  -1,&(outarray[0])+(k*nquad0),1);
805  }
806  }
807  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
808  {
809  //Transposed, Direction A and B positive
810  for (int i=0; i<nquad0; i++)
811  {
812  Vmath::Vcopy(nquad2,&(inarray[0])+i,nquad0*nquad1,
813  &(outarray[0])+(i*nquad2),1);
814  }
815  }
816  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
817  {
818  //Transposed, Direction A negative and B positive
819  for (int i=0; i<nquad0; i++)
820  {
821  Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*nquad1*(nquad2-1)+i,
822  -nquad0*nquad1,&(outarray[0])+(i*nquad2),1);
823  }
824  }
825  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
826  {
827  //Transposed, Direction A positive and B negative
828  for (int i=0; i<nquad0; i++)
829  {
830  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0-1-i),nquad0*nquad1,
831  &(outarray[0])+(i*nquad2),1);
832  }
833  }
834  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
835  {
836  //Transposed, Direction A and B negative
837  for (int i=0; i<nquad0; i++)
838  {
839  Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*nquad1*nquad2+(nquad0-1-i),
840  -nquad0*nquad1,&(outarray[0])+(i*nquad2),1);
841  }
842  }
843  o_tmp = outarray;
844  //interpolate
845  LibUtilities::Interp2D(m_base[0]->GetPointsKey(),
846  m_base[2]->GetPointsKey(), o_tmp,
847  FaceExp->GetBasis(0)->GetPointsKey(),
848  FaceExp->GetBasis(1)->GetPointsKey(),
849  outarray);
850  break;
851  case 2:
852  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
853  {
854  //Directions A and B positive
855  Vmath::Vcopy(nquad0*nquad1,&(inarray[0])+(nquad0-1),
856  nquad0,&(outarray[0]),1);
857  }
858  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
859  {
860  //Direction A negative and B positive
861  for (int k=0; k<nquad2; k++)
862  {
863  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1)+(k*nquad0*nquad1),
864  -nquad0,&(outarray[0])+(k*nquad0),1);
865  }
866  }
867  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
868  {
869  //Direction A positive and B negative
870  for (int k=0; k<nquad2; k++)
871  {
872  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+(nquad0*nquad1*(nquad2-1-k)),
873  nquad0,&(outarray[0])+(k*nquad0),1);
874  }
875  }
876  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
877  {
878  //Direction A negative and B negative
879  for (int k=0; k<nquad2; k++)
880  {
881  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1)+(nquad0*nquad1*(nquad2-1-k)),
882  -nquad0,&(outarray[0])+(k*nquad0),1);
883  }
884  }
885  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
886  {
887  //Transposed, Direction A and B positive
888  for (int j=0; j<nquad1; j++)
889  {
890  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0-1)+(j*nquad0),
891  nquad0*nquad1,&(outarray[0])+(j*nquad2),1);
892  }
893  }
894  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
895  {
896  //Transposed, Direction A negative and B positive
897  for (int j=0; j<nquad0; j++)
898  {
899  Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*nquad1*(nquad2-1)+nquad0+j*nquad0,
900  -nquad0*nquad1,&(outarray[0])+(j*nquad2),1);
901  }
902  }
903  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
904  {
905  //Transposed, Direction A positive and B negative
906  for (int j=0; j<nquad0; j++)
907  {
908  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*nquad1-1-j*nquad0),
909  nquad0*nquad1,&(outarray[0])+(j*nquad2),1);
910  }
911  }
912  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
913  {
914  //Transposed, Direction A and B negative
915  for (int j=0; j<nquad0; j++)
916  {
917  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*nquad1*nquad2-1-j*nquad0),
918  -nquad0*nquad1,&(outarray[0])+(j*nquad2),1);
919  }
920  }
921  o_tmp = outarray;
922  //interpolate
923  LibUtilities::Interp2D(m_base[1]->GetPointsKey(),
924  m_base[2]->GetPointsKey(), o_tmp,
925  FaceExp->GetBasis(0)->GetPointsKey(),
926  FaceExp->GetBasis(1)->GetPointsKey(),
927  outarray);
928  break;
929  case 3:
930  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
931  {
932  //Directions A and B positive
933  for (int k=0; k<nquad2; k++)
934  {
935  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*(nquad1-1))+(k*nquad0*nquad1),
936  1,&(outarray[0])+(k*nquad0),1);
937  }
938  }
939  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
940  {
941  //Direction A negative and B positive
942  for (int k=0; k<nquad2; k++)
943  {
944  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1)+(k*nquad0*nquad1),
945  -1,&(outarray[0])+(k*nquad0),1);
946  }
947  }
948  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
949  {
950  //Direction A positive and B negative
951  for (int k=0; k<nquad2; k++)
952  {
953  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*(nquad1-1))+(nquad0*nquad1*(nquad2-1-k)),
954  1,&(outarray[0])+(k*nquad0),1);
955  }
956  }
957  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
958  {
959  //Direction A negative and B negative
960  for (int k=0; k<nquad2; k++)
961  {
962  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1)+(nquad0*nquad1*(nquad2-1-k)),
963  -1,&(outarray[0])+(k*nquad0),1);
964  }
965  }
966  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
967  {
968  //Transposed, Direction A and B positive
969  for (int i=0; i<nquad0; i++)
970  {
971  Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*(nquad1-1)+i,nquad0*nquad1,
972  &(outarray[0])+(i*nquad2),1);
973  }
974  }
975  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
976  {
977  //Transposed, Direction A negative and B positive
978  for (int i=0; i<nquad0; i++)
979  {
980  Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*(nquad1*nquad2-1)+i,-nquad0*nquad1,
981  &(outarray[0])+(i*nquad2),1);
982  }
983  }
984  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
985  {
986  //Transposed, Direction A positive and B negative
987  for (int i=0; i<nquad0; i++)
988  {
989  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*nquad1-1-i),nquad0*nquad1,
990  &(outarray[0])+(i*nquad2),1);
991  }
992  }
993  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
994  {
995  //Transposed, Direction A and B negative
996  for (int i=0; i<nquad0; i++)
997  {
998  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*nquad1*nquad2-1-i),-nquad0*nquad1,
999  &(outarray[0])+(i*nquad2),1);
1000  }
1001  }
1002  o_tmp = outarray;
1003  //interpolate
1004  LibUtilities::Interp2D(m_base[0]->GetPointsKey(),
1005  m_base[2]->GetPointsKey(), o_tmp,
1006  FaceExp->GetBasis(0)->GetPointsKey(),
1007  FaceExp->GetBasis(1)->GetPointsKey(),
1008  outarray);
1009  break;
1010  case 4:
1011  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
1012  {
1013  //Directions A and B positive
1014  Vmath::Vcopy(nquad0*nquad1,&(inarray[0]),nquad0,&(outarray[0]),1);
1015  }
1016  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
1017  {
1018  //Direction A negative and B positive
1019  for (int k=0; k<nquad2; k++)
1020  {
1021  Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1)+(k*nquad0*nquad1),
1022  -nquad0,&(outarray[0])+(k*nquad0),1);
1023  }
1024  }
1025  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
1026  {
1027  //Direction A positive and B negative
1028  for (int k=0; k<nquad2; k++)
1029  {
1030  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1*(nquad2-1-k)),
1031  nquad0,&(outarray[0])+(k*nquad0),1);
1032  }
1033  }
1034  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
1035  {
1036  //Direction A negative and B negative
1037  for (int k=0; k<nquad2; k++)
1038  {
1039  Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1)+(nquad0*nquad1*(nquad2-1-k)),
1040  -nquad0,&(outarray[0])+(k*nquad0),1);
1041  }
1042  }
1043  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
1044  {
1045  //Transposed, Direction A and B positive
1046  for (int j=0; j<nquad0; j++)
1047  {
1048  Vmath::Vcopy(nquad2,&(inarray[0])+j*nquad0,nquad0*nquad1,
1049  &(outarray[0])+(j*nquad2),1);
1050  }
1051  }
1052  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
1053  {
1054  //Transposed, Direction A negative and B positive
1055  for (int j=0; j<nquad0; j++)
1056  {
1057  Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*nquad1*(nquad2-1)+j*nquad0,
1058  -nquad0*nquad1,&(outarray[0])+(j*nquad2),1);
1059  }
1060  }
1061  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
1062  {
1063  //Transposed, Direction A positive and B negative
1064  for (int j=0; j<nquad0; j++)
1065  {
1066  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*(nquad1-1)-j*nquad0),
1067  nquad0*nquad1,&(outarray[0])+(j*nquad2),1);
1068  }
1069  }
1070  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
1071  {
1072  //Transposed, Direction A and B negative
1073  for (int j=0; j<nquad0; j++)
1074  {
1075  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*(nquad1*nquad2-1)-j*nquad0),
1076  -nquad0*nquad1,&(outarray[0])+(j*nquad2),1);
1077  }
1078  }
1079  o_tmp = outarray;
1080  //interpolate
1081  LibUtilities::Interp2D(m_base[1]->GetPointsKey(),
1082  m_base[2]->GetPointsKey(), o_tmp,
1083  FaceExp->GetBasis(0)->GetPointsKey(),
1084  FaceExp->GetBasis(1)->GetPointsKey(),
1085  outarray);
1086  break;
1087  case 5:
1088  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
1089  {
1090  //Directions A and B positive
1091  Vmath::Vcopy(nquad0*nquad1,&(inarray[0])+nquad0*nquad1*(nquad2-1),1,&(outarray[0]),1);
1092  }
1093  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
1094  {
1095  //Direction A negative and B positive
1096  for (int j=0; j<nquad1; j++)
1097  {
1098  Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*nquad1*(nquad2-1)+(nquad0-1+j*nquad0),
1099  -1,&(outarray[0])+(j*nquad0),1);
1100  }
1101  }
1102  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
1103  {
1104  //Direction A positive and B negative
1105  for (int j=0; j<nquad1; j++)
1106  {
1107  Vmath::Vcopy(nquad0,&(inarray[0])+((nquad0*nquad1*nquad2-1)-(nquad0-1)-j*nquad0),
1108  1,&(outarray[0])+(j*nquad0),1);
1109  }
1110  }
1111  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
1112  {
1113  //Direction A negative and B negative
1114  for (int j=0; j<nquad1; j++)
1115  {
1116  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1*nquad2-1-j*nquad0),
1117  -1,&(outarray[0])+(j*nquad0),1);
1118  }
1119  }
1120  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
1121  {
1122  //Transposed, Direction A and B positive
1123  for (int i=0; i<nquad0; i++)
1124  {
1125  Vmath::Vcopy(nquad1,&(inarray[0])+nquad0*nquad1*(nquad2-1)+i,nquad0,
1126  &(outarray[0])+(i*nquad1),1);
1127  }
1128  }
1129  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
1130  {
1131  //Transposed, Direction A negative and B positive
1132  for (int i=0; i<nquad0; i++)
1133  {
1134  Vmath::Vcopy(nquad1,&(inarray[0])+nquad0*(nquad1*nquad2-1)+i,-nquad0,
1135  &(outarray[0])+(i*nquad1),1);
1136  }
1137  }
1138  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
1139  {
1140  //Transposed, Direction A positive and B negative
1141  for (int i=0; i<nquad0; i++)
1142  {
1143  Vmath::Vcopy(nquad1,&(inarray[0])+nquad0*nquad1*(nquad2-1)+(nquad0-1-i),
1144  nquad0,&(outarray[0])+(i*nquad1),1);
1145  }
1146  }
1147  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
1148  {
1149  //Transposed, Direction A and B negative
1150  for (int i=0; i<nquad0; i++)
1151  {
1152  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0*nquad1*nquad2-1-i),-nquad0,
1153  &(outarray[0])+(i*nquad1),1);
1154  }
1155  }
1156  o_tmp = outarray;
1157  //interpolate
1158  LibUtilities::Interp2D(m_base[0]->GetPointsKey(),
1159  m_base[1]->GetPointsKey(), o_tmp,
1160  FaceExp->GetBasis(0)->GetPointsKey(),
1161  FaceExp->GetBasis(1)->GetPointsKey(),
1162  outarray);
1163  break;
1164  default:
1165  ASSERTL0(false,"face value (> 5) is out of range");
1166  break;
1167  }
1168  }
1169 
1170 
1171  void HexExp::v_ComputeFaceNormal(const int face)
1172  {
1173  int i;
1174  const SpatialDomains::GeomFactorsSharedPtr & geomFactors =
1175  GetGeom()->GetMetricInfo();
1176  SpatialDomains::GeomType type = geomFactors->GetGtype();
1177  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1178  const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys);
1179  const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys);
1180 
1181  int nqe0 = m_base[0]->GetNumPoints();
1182  int nqe1 = m_base[1]->GetNumPoints();
1183  int nqe2 = m_base[2]->GetNumPoints();
1184  int nqe01 = nqe0*nqe1;
1185  int nqe02 = nqe0*nqe2;
1186  int nqe12 = nqe1*nqe2;
1187 
1188  int nqe;
1189  if (face == 0 || face == 5)
1190  {
1191  nqe = nqe01;
1192  }
1193  else if (face == 1 || face == 3)
1194  {
1195  nqe = nqe02;
1196  }
1197  else
1198  {
1199  nqe = nqe12;
1200  }
1201 
1202  int vCoordDim = GetCoordim();
1203 
1204  m_faceNormals[face] = Array<OneD, Array<OneD, NekDouble> >(vCoordDim);
1205  Array<OneD, Array<OneD, NekDouble> > &normal = m_faceNormals[face];
1206  for (i = 0; i < vCoordDim; ++i)
1207  {
1208  normal[i] = Array<OneD, NekDouble>(nqe);
1209  }
1210  // Regular geometry case
1211  if((type == SpatialDomains::eRegular)||(type == SpatialDomains::eMovingRegular))
1212  {
1213  NekDouble fac;
1214  // Set up normals
1215  switch(face)
1216  {
1217  case 0:
1218  for(i = 0; i < vCoordDim; ++i)
1219  {
1220  Vmath::Fill(nqe,-df[3*i+2][0],normal[i],1);
1221  }
1222  break;
1223  case 1:
1224  for(i = 0; i < vCoordDim; ++i)
1225  {
1226  Vmath::Fill(nqe,-df[3*i+1][0],normal[i],1);
1227  }
1228  break;
1229  case 2:
1230  for(i = 0; i < vCoordDim; ++i)
1231  {
1232  Vmath::Fill(nqe,df[3*i][0],normal[i],1);
1233  }
1234  break;
1235  case 3:
1236  for(i = 0; i < vCoordDim; ++i)
1237  {
1238  Vmath::Fill(nqe,df[3*i+1][0],normal[i],1);
1239  }
1240  break;
1241  case 4:
1242  for(i = 0; i < vCoordDim; ++i)
1243  {
1244  Vmath::Fill(nqe,-df[3*i][0],normal[i],1);
1245  }
1246  break;
1247  case 5:
1248  for(i = 0; i < vCoordDim; ++i)
1249  {
1250  Vmath::Fill(nqe,df[3*i+2][0],normal[i],1);
1251  }
1252  break;
1253  default:
1254  ASSERTL0(false,"face is out of range (edge < 5)");
1255  }
1256 
1257  // normalise
1258  fac = 0.0;
1259  for(i =0 ; i < vCoordDim; ++i)
1260  {
1261  fac += normal[i][0]*normal[i][0];
1262  }
1263  fac = 1.0/sqrt(fac);
1264  for (i = 0; i < vCoordDim; ++i)
1265  {
1266  Vmath::Smul(nqe,fac,normal[i],1,normal[i],1);
1267  }
1268 
1269  }
1270  else // Set up deformed normals
1271  {
1272  int j, k;
1273 
1274  Array<OneD,NekDouble> work(nqe,0.0);
1275 
1276  // Extract Jacobian along face and recover local
1277  // derivates (dx/dr) for polynomial interpolation by
1278  // multiplying m_gmat by jacobian
1279  switch(face)
1280  {
1281  case 0:
1282  for(j = 0; j < nqe; ++j)
1283  {
1284  normal[0][j] = -df[2][j]*jac[j];
1285  normal[1][j] = -df[5][j]*jac[j];
1286  normal[2][j] = -df[8][j]*jac[j];
1287  }
1288  break;
1289  case 1:
1290  for (j = 0; j < nqe0; ++j)
1291  {
1292  for(k = 0; k < nqe2; ++k)
1293  {
1294  int idx = j + nqe01*k;
1295  normal[0][j+k*nqe0] = -df[1][idx]*jac[idx];
1296  normal[1][j+k*nqe0] = -df[4][idx]*jac[idx];
1297  normal[2][j+k*nqe0] = -df[7][idx]*jac[idx];
1298  }
1299  }
1300  break;
1301  case 2:
1302  for (j = 0; j < nqe1; ++j)
1303  {
1304  for(k = 0; k < nqe2; ++k)
1305  {
1306  int idx = nqe0-1+nqe0*j+nqe01*k;
1307  normal[0][j+k*nqe0] = df[0][idx]*jac[idx];
1308  normal[1][j+k*nqe0] = df[3][idx]*jac[idx];
1309  normal[2][j+k*nqe0] = df[6][idx]*jac[idx];
1310  }
1311  }
1312  break;
1313  case 3:
1314  for (j = 0; j < nqe0; ++j)
1315  {
1316  for(k = 0; k < nqe2; ++k)
1317  {
1318  int idx = nqe0*(nqe1-1)+j+nqe01*k;
1319  normal[0][j+k*nqe0] = df[1][idx]*jac[idx];
1320  normal[1][j+k*nqe0] = df[4][idx]*jac[idx];
1321  normal[2][j+k*nqe0] = df[7][idx]*jac[idx];
1322  }
1323  }
1324  break;
1325  case 4:
1326  for (j = 0; j < nqe0; ++j)
1327  {
1328  for(k = 0; k < nqe2; ++k)
1329  {
1330  int idx = j*nqe0+nqe01*k;
1331  normal[0][j+k*nqe0] = -df[0][idx]*jac[idx];
1332  normal[1][j+k*nqe0] = -df[3][idx]*jac[idx];
1333  normal[2][j+k*nqe0] = -df[6][idx]*jac[idx];
1334  }
1335  }
1336  break;
1337  case 5:
1338  for (j = 0; j < nqe01; ++j)
1339  {
1340  int idx = j+nqe01*(nqe2-1);
1341  normal[0][j] = df[2][idx]*jac[idx];
1342  normal[1][j] = df[5][idx]*jac[idx];
1343  normal[2][j] = df[8][idx]*jac[idx];
1344  }
1345  break;
1346  default:
1347  ASSERTL0(false,"face is out of range (face < 5)");
1348  }
1349 
1350  Vmath::Sdiv(nqe,1.0,&jac[0],1,&work[0],1);
1351 
1352  // interpolate
1353  for(i = 0; i < GetCoordim(); ++i)
1354  {
1355  Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
1356  }
1357 
1358  //normalise normal vectors
1359  Vmath::Zero(nqe,work,1);
1360  for(i = 0; i < GetCoordim(); ++i)
1361  {
1362  Vmath::Vvtvp(nqe,normal[i],1, normal[i],1,work,1,work,1);
1363  }
1364 
1365  Vmath::Vsqrt(nqe,work,1,work,1);
1366  Vmath::Sdiv(nqe,1.0,work,1,work,1);
1367 
1368  for(i = 0; i < GetCoordim(); ++i)
1369  {
1370  Vmath::Vmul(nqe,normal[i],1,work,1,normal[i],1);
1371  }
1372  }
1373  }
1374 
1375  //-----------------------------
1376  // Operator creation functions
1377  //-----------------------------
1378  void HexExp::v_MassMatrixOp(
1379  const Array<OneD, const NekDouble> &inarray,
1380  Array<OneD,NekDouble> &outarray,
1381  const StdRegions::StdMatrixKey &mkey)
1382  {
1383  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
1384  }
1385 
1386  void HexExp::v_LaplacianMatrixOp(
1387  const Array<OneD, const NekDouble> &inarray,
1388  Array<OneD,NekDouble> &outarray,
1389  const StdRegions::StdMatrixKey &mkey)
1390  {
1391  HexExp::v_LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1392  }
1393 
1394  void HexExp::v_LaplacianMatrixOp(
1395  const int k1,
1396  const int k2,
1397  const Array<OneD, const NekDouble> &inarray,
1398  Array<OneD,NekDouble> &outarray,
1399  const StdRegions::StdMatrixKey &mkey)
1400  {
1401  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,
1402  mkey);
1403  }
1404 
1405  void HexExp::v_WeakDerivMatrixOp(
1406  const int i,
1407  const Array<OneD, const NekDouble> &inarray,
1408  Array<OneD,NekDouble> &outarray,
1409  const StdRegions::StdMatrixKey &mkey)
1410  {
1411  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
1412  }
1413 
1414  void HexExp::v_WeakDirectionalDerivMatrixOp(
1415  const Array<OneD, const NekDouble> &inarray,
1416  Array<OneD,NekDouble> &outarray,
1417  const StdRegions::StdMatrixKey &mkey)
1418  {
1419  StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray,
1420  outarray,mkey);
1421  }
1422 
1423  void HexExp::v_MassLevelCurvatureMatrixOp(
1424  const Array<OneD, const NekDouble> &inarray,
1425  Array<OneD,NekDouble> &outarray,
1426  const StdRegions::StdMatrixKey &mkey)
1427  {
1428  StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray,
1429  outarray,mkey);
1430  }
1431 
1432  void HexExp::v_HelmholtzMatrixOp(
1433  const Array<OneD, const NekDouble> &inarray,
1434  Array<OneD,NekDouble> &outarray,
1435  const StdRegions::StdMatrixKey &mkey)
1436  {
1437  HexExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1438  }
1439 
1440 
1441  void HexExp::v_GeneralMatrixOp_MatOp(
1442  const Array<OneD, const NekDouble> &inarray,
1443  Array<OneD,NekDouble> &outarray,
1444  const StdRegions::StdMatrixKey &mkey)
1445  {
1446  //int nConsts = mkey.GetNconstants();
1447  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1448 
1449 // switch(nConsts)
1450 // {
1451 // case 0:
1452 // {
1453 // mat = GetLocMatrix(mkey.GetMatrixType());
1454 // }
1455 // break;
1456 // case 1:
1457 // {
1458 // mat = GetLocMatrix(mkey.GetMatrixType(),mkey.GetConstant(0));
1459 // }
1460 // break;
1461 // case 2:
1462 // {
1463 // mat = GetLocMatrix(mkey.GetMatrixType(),mkey.GetConstant(0),mkey.GetConstant(1));
1464 // }
1465 // break;
1466 //
1467 // default:
1468 // {
1469 // NEKERROR(ErrorUtil::efatal, "Unknown number of constants");
1470 // }
1471 // break;
1472 // }
1473 
1474  if(inarray.get() == outarray.get())
1475  {
1476  Array<OneD,NekDouble> tmp(m_ncoeffs);
1477  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1478 
1479  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1480  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1481  }
1482  else
1483  {
1484  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1485  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1486  }
1487  }
1488 
1489  /**
1490  * This function is used to compute exactly the advective numerical flux
1491  * on the interface of two elements with different expansions, hence an
1492  * appropriate number of Gauss points has to be used. The number of
1493  * Gauss points has to be equal to the number used by the highest
1494  * polynomial degree of the two adjacent elements
1495  *
1496  * @param numMin Is the reduced polynomial order
1497  * @param inarray Input array of coefficients
1498  * @param dumpVar Output array of reduced coefficients.
1499  */
1500  void HexExp::v_ReduceOrderCoeffs(
1501  int numMin,
1502  const Array<OneD, const NekDouble> &inarray,
1503  Array<OneD, NekDouble> &outarray)
1504  {
1505  int n_coeffs = inarray.num_elements();
1506  int nmodes0 = m_base[0]->GetNumModes();
1507  int nmodes1 = m_base[1]->GetNumModes();
1508  int nmodes2 = m_base[2]->GetNumModes();
1509  int numMax = nmodes0;
1510 
1511  Array<OneD, NekDouble> coeff (n_coeffs);
1512  Array<OneD, NekDouble> coeff_tmp1(nmodes0*nmodes1, 0.0);
1513  Array<OneD, NekDouble> coeff_tmp2(n_coeffs, 0.0);
1514  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
1515 
1516  Vmath::Vcopy(n_coeffs,inarray,1,coeff_tmp2,1);
1517 
1518  const LibUtilities::PointsKey Pkey0(
1519  nmodes0, LibUtilities::eGaussLobattoLegendre);
1520  const LibUtilities::PointsKey Pkey1(
1521  nmodes1, LibUtilities::eGaussLobattoLegendre);
1522  const LibUtilities::PointsKey Pkey2(
1523  nmodes2, LibUtilities::eGaussLobattoLegendre);
1524 
1525  LibUtilities::BasisKey b0(
1526  m_base[0]->GetBasisType(), nmodes0, Pkey0);
1527  LibUtilities::BasisKey b1(
1528  m_base[1]->GetBasisType(), nmodes1, Pkey1);
1529  LibUtilities::BasisKey b2(
1530  m_base[2]->GetBasisType(), nmodes2, Pkey2);
1531  LibUtilities::BasisKey bortho0(
1532  LibUtilities::eOrtho_A, nmodes0, Pkey0);
1533  LibUtilities::BasisKey bortho1(
1534  LibUtilities::eOrtho_A, nmodes1, Pkey1);
1535  LibUtilities::BasisKey bortho2(
1536  LibUtilities::eOrtho_A, nmodes2, Pkey2);
1537 
1538  LibUtilities::InterpCoeff3D(
1539  b0, b1, b2, coeff_tmp2,
1540  bortho0, bortho1, bortho2, coeff);
1541 
1542  Vmath::Zero(n_coeffs, coeff_tmp2, 1);
1543 
1544  int cnt = 0, cnt2 = 0;
1545 
1546  for (int u = 0; u < numMin+1; ++u)
1547  {
1548  for (int i = 0; i < numMin; ++i)
1549  {
1550  Vmath::Vcopy(numMin,
1551  tmp = coeff+cnt+cnt2,1,
1552  tmp2 = coeff_tmp1+cnt,1);
1553 
1554  cnt = i*numMax;
1555  }
1556 
1557  Vmath::Vcopy(nmodes0*nmodes1,
1558  tmp3 = coeff_tmp1,1,
1559  tmp4 = coeff_tmp2+cnt2,1);
1560 
1561  cnt2 = u*nmodes0*nmodes1;
1562  }
1563 
1564  LibUtilities::InterpCoeff3D(
1565  bortho0, bortho1, bortho2, coeff_tmp2,
1566  b0, b1, b2, outarray);
1567  }
1568 
1569  //-----------------------------
1570  // Matrix creation functions
1571  //-----------------------------
1572  DNekMatSharedPtr HexExp::v_GenMatrix(
1573  const StdRegions::StdMatrixKey &mkey)
1574  {
1575  DNekMatSharedPtr returnval;
1576 
1577  switch(mkey.GetMatrixType())
1578  {
1579  case StdRegions::eHybridDGHelmholtz:
1580  case StdRegions::eHybridDGLamToU:
1581  case StdRegions::eHybridDGLamToQ0:
1582  case StdRegions::eHybridDGLamToQ1:
1583  case StdRegions::eHybridDGLamToQ2:
1584  case StdRegions::eHybridDGHelmBndLam:
1585  case StdRegions::eInvLaplacianWithUnityMean:
1586  returnval = Expansion3D::v_GenMatrix(mkey);
1587  break;
1588  default:
1589  returnval = StdHexExp::v_GenMatrix(mkey);
1590  }
1591 
1592  return returnval;
1593  }
1594 
1595 
1596  DNekMatSharedPtr HexExp::v_CreateStdMatrix(
1597  const StdRegions::StdMatrixKey &mkey)
1598  {
1599  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1600  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1601  LibUtilities::BasisKey bkey2 = m_base[2]->GetBasisKey();
1602 
1603  StdRegions::StdHexExpSharedPtr tmp = MemoryManager<StdHexExp>
1604  ::AllocateSharedPtr(bkey0,bkey1,bkey2);
1605 
1606  return tmp->GetStdMatrix(mkey);
1607  }
1608 
1609 
1610  DNekScalMatSharedPtr HexExp::CreateMatrix(const MatrixKey &mkey)
1611  {
1612  DNekScalMatSharedPtr returnval;
1613  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1614 
1615  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,
1616  "Geometric information is not set up");
1617 
1618  switch(mkey.GetMatrixType())
1619  {
1620  case StdRegions::eMass:
1621  {
1622  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1623  mkey.GetNVarCoeff())
1624  {
1625  NekDouble one = 1.0;
1626  DNekMatSharedPtr mat = GenMatrix(mkey);
1627  returnval = MemoryManager<DNekScalMat>
1628  ::AllocateSharedPtr(one,mat);
1629  }
1630  else
1631  {
1632  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1633  DNekMatSharedPtr mat
1634  = GetStdMatrix(mkey);
1635  returnval = MemoryManager<DNekScalMat>
1636  ::AllocateSharedPtr(jac,mat);
1637  }
1638  }
1639  break;
1640  case StdRegions::eInvMass:
1641  {
1642  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1643  {
1644  NekDouble one = 1.0;
1645  StdRegions::StdMatrixKey masskey(StdRegions::eMass,
1646  DetShapeType(), *this);
1647  DNekMatSharedPtr mat = GenMatrix(masskey);
1648  mat->Invert();
1649 
1650  returnval = MemoryManager<DNekScalMat>
1651  ::AllocateSharedPtr(one,mat);
1652  }
1653  else
1654  {
1655  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1656  DNekMatSharedPtr mat
1657  = GetStdMatrix(mkey);
1658  returnval = MemoryManager<DNekScalMat>
1659  ::AllocateSharedPtr(fac,mat);
1660  }
1661  }
1662  break;
1663  case StdRegions::eWeakDeriv0:
1664  case StdRegions::eWeakDeriv1:
1665  case StdRegions::eWeakDeriv2:
1666  {
1667  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1668  mkey.GetNVarCoeff())
1669  {
1670  NekDouble one = 1.0;
1671  DNekMatSharedPtr mat = GenMatrix(mkey);
1672 
1673  returnval = MemoryManager<DNekScalMat>
1674  ::AllocateSharedPtr(one,mat);
1675  }
1676  else
1677  {
1678  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1679  Array<TwoD, const NekDouble> df
1680  = m_metricinfo->GetDerivFactors(ptsKeys);
1681  int dir = 0;
1682 
1683  switch(mkey.GetMatrixType())
1684  {
1685  case StdRegions::eWeakDeriv0:
1686  dir = 0;
1687  break;
1688  case StdRegions::eWeakDeriv1:
1689  dir = 1;
1690  break;
1691  case StdRegions::eWeakDeriv2:
1692  dir = 2;
1693  break;
1694  default:
1695  break;
1696  }
1697 
1698  MatrixKey deriv0key(StdRegions::eWeakDeriv0,
1699  mkey.GetShapeType(), *this);
1700  MatrixKey deriv1key(StdRegions::eWeakDeriv1,
1701  mkey.GetShapeType(), *this);
1702  MatrixKey deriv2key(StdRegions::eWeakDeriv2,
1703  mkey.GetShapeType(), *this);
1704 
1705  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
1706  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
1707  DNekMat &deriv2 = *GetStdMatrix(deriv2key);
1708 
1709  int rows = deriv0.GetRows();
1710  int cols = deriv1.GetColumns();
1711 
1712  DNekMatSharedPtr WeakDeriv = MemoryManager<DNekMat>
1713  ::AllocateSharedPtr(rows,cols);
1714 
1715  (*WeakDeriv) = df[3*dir ][0]*deriv0
1716  + df[3*dir+1][0]*deriv1
1717  + df[3*dir+2][0]*deriv2;
1718 
1719  returnval = MemoryManager<DNekScalMat>
1720  ::AllocateSharedPtr(jac,WeakDeriv);
1721  }
1722  }
1723  break;
1724  case StdRegions::eLaplacian:
1725  {
1726  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1727  mkey.GetNVarCoeff()||
1728  mkey.ConstFactorExists(
1729  StdRegions::eFactorSVVCutoffRatio))
1730  {
1731  NekDouble one = 1.0;
1732  DNekMatSharedPtr mat = GenMatrix(mkey);
1733 
1734  returnval = MemoryManager<DNekScalMat>
1735  ::AllocateSharedPtr(one,mat);
1736  }
1737  else
1738  {
1739  MatrixKey lap00key(StdRegions::eLaplacian00,
1740  mkey.GetShapeType(), *this);
1741  MatrixKey lap01key(StdRegions::eLaplacian01,
1742  mkey.GetShapeType(), *this);
1743  MatrixKey lap02key(StdRegions::eLaplacian02,
1744  mkey.GetShapeType(), *this);
1745  MatrixKey lap11key(StdRegions::eLaplacian11,
1746  mkey.GetShapeType(), *this);
1747  MatrixKey lap12key(StdRegions::eLaplacian12,
1748  mkey.GetShapeType(), *this);
1749  MatrixKey lap22key(StdRegions::eLaplacian22,
1750  mkey.GetShapeType(), *this);
1751 
1752  DNekMat &lap00 = *GetStdMatrix(lap00key);
1753  DNekMat &lap01 = *GetStdMatrix(lap01key);
1754  DNekMat &lap02 = *GetStdMatrix(lap02key);
1755  DNekMat &lap11 = *GetStdMatrix(lap11key);
1756  DNekMat &lap12 = *GetStdMatrix(lap12key);
1757  DNekMat &lap22 = *GetStdMatrix(lap22key);
1758 
1759  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1760  Array<TwoD, const NekDouble> gmat
1761  = m_metricinfo->GetGmat(ptsKeys);
1762 
1763  int rows = lap00.GetRows();
1764  int cols = lap00.GetColumns();
1765 
1766  DNekMatSharedPtr lap = MemoryManager<DNekMat>
1767  ::AllocateSharedPtr(rows,cols);
1768 
1769  (*lap) = gmat[0][0]*lap00
1770  + gmat[4][0]*lap11
1771  + gmat[8][0]*lap22
1772  + gmat[3][0]*(lap01 + Transpose(lap01))
1773  + gmat[6][0]*(lap02 + Transpose(lap02))
1774  + gmat[7][0]*(lap12 + Transpose(lap12));
1775 
1776  returnval = MemoryManager<DNekScalMat>
1777  ::AllocateSharedPtr(jac,lap);
1778  }
1779  }
1780  break;
1781  case StdRegions::eHelmholtz:
1782  {
1783  NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
1784  MatrixKey masskey(StdRegions::eMass,
1785  mkey.GetShapeType(), *this);
1786  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1787  MatrixKey lapkey(StdRegions::eLaplacian,
1788  mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1789  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1790 
1791  int rows = LapMat.GetRows();
1792  int cols = LapMat.GetColumns();
1793 
1794  DNekMatSharedPtr helm = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1795 
1796  NekDouble one = 1.0;
1797  (*helm) = LapMat + lambda*MassMat;
1798 
1799  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
1800  }
1801  break;
1802  case StdRegions::eIProductWRTBase:
1803  {
1804  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1805  {
1806  NekDouble one = 1.0;
1807  DNekMatSharedPtr mat = GenMatrix(mkey);
1808  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1809  }
1810  else
1811  {
1812  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1813  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1814  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1815  }
1816  }
1817  break;
1818  case StdRegions::eHybridDGHelmholtz:
1819  case StdRegions::eHybridDGLamToU:
1820  case StdRegions::eHybridDGLamToQ0:
1821  case StdRegions::eHybridDGLamToQ1:
1822  case StdRegions::eHybridDGLamToQ2:
1823  case StdRegions::eHybridDGHelmBndLam:
1824  case StdRegions::eInvLaplacianWithUnityMean:
1825  {
1826  NekDouble one = 1.0;
1827 
1828  DNekMatSharedPtr mat = GenMatrix(mkey);
1829  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1830  }
1831  break;
1832  case StdRegions::eInvHybridDGHelmholtz:
1833  {
1834  NekDouble one = 1.0;
1835 
1836 // StdRegions::StdMatrixKey hkey(StdRegions::eHybridDGHelmholtz,
1837 // DetShapeType(),*this,
1838 // mkey.GetConstant(0),
1839 // mkey.GetConstant(1));
1840  MatrixKey hkey(StdRegions::eHybridDGHelmholtz, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1841  DNekMatSharedPtr mat = GenMatrix(hkey);
1842 
1843  mat->Invert();
1844  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1845  }
1846  break;
1847  case StdRegions::ePreconLinearSpace:
1848  {
1849  NekDouble one = 1.0;
1850  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1851  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1852  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1853  DNekMatSharedPtr R=BuildVertexMatrix(A);
1854 
1855  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1856  }
1857  break;
1858  case StdRegions::ePreconLinearSpaceMass:
1859  {
1860  NekDouble one = 1.0;
1861  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1862  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1863  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1864  DNekMatSharedPtr R=BuildVertexMatrix(A);
1865 
1866  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1867  }
1868  break;
1869  case StdRegions::ePreconR:
1870  {
1871  NekDouble one = 1.0;
1872  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this,mkey.GetConstFactors(), mkey.GetVarCoeffs());
1873  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1874  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1875 
1876  DNekScalMatSharedPtr Atmp;
1877  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1878 
1879  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1880  }
1881  break;
1882  case StdRegions::ePreconRT:
1883  {
1884  NekDouble one = 1.0;
1885  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this,mkey.GetConstFactors(), mkey.GetVarCoeffs());
1886  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1887  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1888 
1889  DNekScalMatSharedPtr Atmp;
1890  DNekMatSharedPtr RT=BuildTransformationMatrix(A,mkey.GetMatrixType());
1891 
1892  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,RT);
1893  }
1894  break;
1895  case StdRegions::ePreconRMass:
1896  {
1897  NekDouble one = 1.0;
1898  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1899  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1900  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1901 
1902  DNekScalMatSharedPtr Atmp;
1903  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1904 
1905  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1906  }
1907  break;
1908  case StdRegions::ePreconRTMass:
1909  {
1910  NekDouble one = 1.0;
1911  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1912  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1913  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1914 
1915  DNekScalMatSharedPtr Atmp;
1916  DNekMatSharedPtr RT=BuildTransformationMatrix(A,mkey.GetMatrixType());
1917 
1918  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,RT);
1919  }
1920  break;
1921  default:
1922  {
1923  NekDouble one = 1.0;
1924  DNekMatSharedPtr mat = GenMatrix(mkey);
1925 
1926  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1927  }
1928  break;
1929  }
1930 
1931  return returnval;
1932  }
1933 
1934 
1935  DNekScalBlkMatSharedPtr HexExp::CreateStaticCondMatrix(const MatrixKey &mkey)
1936  {
1937  DNekScalBlkMatSharedPtr returnval;
1938 
1939  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
1940 
1941  // set up block matrix system
1942  unsigned int nbdry = NumBndryCoeffs();
1943  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
1944  unsigned int exp_size[] = {nbdry,nint};
1945  unsigned int nblks = 2;
1946  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks,nblks,exp_size,exp_size); //Really need a constructor which takes Arrays
1947  NekDouble factor = 1.0;
1948 
1949  switch(mkey.GetMatrixType())
1950  {
1951  case StdRegions::eLaplacian:
1952  case StdRegions::eHelmholtz: // special case since Helmholtz not defined in StdRegions
1953 
1954  // use Deformed case for both regular and deformed geometries
1955  factor = 1.0;
1956  goto UseLocRegionsMatrix;
1957  break;
1958  default:
1959  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1960  mkey.GetNVarCoeff())
1961  {
1962  factor = 1.0;
1963  goto UseLocRegionsMatrix;
1964  }
1965  else
1966  {
1967  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1968  factor = mat->Scale();
1969  goto UseStdRegionsMatrix;
1970  }
1971  break;
1972  UseStdRegionsMatrix:
1973  {
1974  NekDouble invfactor = 1.0/factor;
1975  NekDouble one = 1.0;
1976  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
1977  DNekScalMatSharedPtr Atmp;
1978  DNekMatSharedPtr Asubmat;
1979 
1980  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
1981  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
1982  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
1983  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
1984  }
1985  break;
1986  UseLocRegionsMatrix:
1987  {
1988  int i,j;
1989  NekDouble invfactor = 1.0/factor;
1990  NekDouble one = 1.0;
1991  DNekScalMat &mat = *GetLocMatrix(mkey);
1992  DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
1993  DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
1994  DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
1995  DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
1996 
1997  Array<OneD,unsigned int> bmap(nbdry);
1998  Array<OneD,unsigned int> imap(nint);
1999  GetBoundaryMap(bmap);
2000  GetInteriorMap(imap);
2001 
2002  for(i = 0; i < nbdry; ++i)
2003  {
2004  for(j = 0; j < nbdry; ++j)
2005  {
2006  (*A)(i,j) = mat(bmap[i],bmap[j]);
2007  }
2008 
2009  for(j = 0; j < nint; ++j)
2010  {
2011  (*B)(i,j) = mat(bmap[i],imap[j]);
2012  }
2013  }
2014 
2015  for(i = 0; i < nint; ++i)
2016  {
2017  for(j = 0; j < nbdry; ++j)
2018  {
2019  (*C)(i,j) = mat(imap[i],bmap[j]);
2020  }
2021 
2022  for(j = 0; j < nint; ++j)
2023  {
2024  (*D)(i,j) = mat(imap[i],imap[j]);
2025  }
2026  }
2027 
2028  // Calculate static condensed system
2029  if(nint)
2030  {
2031  D->Invert();
2032  (*B) = (*B)*(*D);
2033  (*A) = (*A) - (*B)*(*C);
2034  }
2035 
2036  DNekScalMatSharedPtr Atmp;
2037 
2038  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
2039  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
2040  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
2041  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
2042 
2043  }
2044  }
2045  return returnval;
2046  }
2047 
2048 
2049  DNekScalMatSharedPtr HexExp::v_GetLocMatrix(const MatrixKey &mkey)
2050  {
2051  return m_matrixManager[mkey];
2052  }
2053 
2054 
2055  DNekScalBlkMatSharedPtr HexExp::v_GetLocStaticCondMatrix(
2056  const MatrixKey &mkey)
2057  {
2058  return m_staticCondMatrixManager[mkey];
2059  }
2060 
2061  void HexExp::v_DropLocStaticCondMatrix(const MatrixKey &mkey)
2062  {
2063  m_staticCondMatrixManager.DeleteObject(mkey);
2064  }
2065 
2066  void HexExp::v_LaplacianMatrixOp_MatFree_Kernel(
2067  const Array<OneD, const NekDouble> &inarray,
2068  Array<OneD, NekDouble> &outarray,
2069  Array<OneD, NekDouble> &wsp)
2070  {
2071  // This implementation is only valid when there are no
2072  // coefficients associated to the Laplacian operator
2073  if (m_metrics.count(MetricLaplacian00) == 0)
2074  {
2075  ComputeLaplacianMetric();
2076  }
2077 
2078  int nquad0 = m_base[0]->GetNumPoints();
2079  int nquad1 = m_base[1]->GetNumPoints();
2080  int nquad2 = m_base[2]->GetNumPoints();
2081  int nqtot = nquad0*nquad1*nquad2;
2082 
2083  ASSERTL1(wsp.num_elements() >= 6*nqtot,
2084  "Insufficient workspace size.");
2085 
2086  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
2087  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
2088  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata();
2089  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
2090  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
2091  const Array<OneD, const NekDouble>& dbase2 = m_base[2]->GetDbdata();
2092  const Array<OneD, const NekDouble>& metric00 = m_metrics[MetricLaplacian00];
2093  const Array<OneD, const NekDouble>& metric01 = m_metrics[MetricLaplacian01];
2094  const Array<OneD, const NekDouble>& metric02 = m_metrics[MetricLaplacian02];
2095  const Array<OneD, const NekDouble>& metric11 = m_metrics[MetricLaplacian11];
2096  const Array<OneD, const NekDouble>& metric12 = m_metrics[MetricLaplacian12];
2097  const Array<OneD, const NekDouble>& metric22 = m_metrics[MetricLaplacian22];
2098 
2099  // Allocate temporary storage
2100  Array<OneD,NekDouble> wsp0(wsp);
2101  Array<OneD,NekDouble> wsp1(wsp+1*nqtot);
2102  Array<OneD,NekDouble> wsp2(wsp+2*nqtot);
2103  Array<OneD,NekDouble> wsp3(wsp+3*nqtot);
2104  Array<OneD,NekDouble> wsp4(wsp+4*nqtot);
2105  Array<OneD,NekDouble> wsp5(wsp+5*nqtot);
2106 
2107  StdExpansion3D::PhysTensorDeriv(inarray,wsp0,wsp1,wsp2);
2108 
2109  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
2110  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
2111  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
2112  // especially for this purpose
2113  Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp0[0],1,&metric01[0],1,&wsp1[0],1,&wsp3[0],1);
2114  Vmath::Vvtvp (nqtot,&metric02[0],1,&wsp2[0],1,&wsp3[0],1,&wsp3[0],1);
2115  Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp0[0],1,&metric11[0],1,&wsp1[0],1,&wsp4[0],1);
2116  Vmath::Vvtvp (nqtot,&metric12[0],1,&wsp2[0],1,&wsp4[0],1,&wsp4[0],1);
2117  Vmath::Vvtvvtp(nqtot,&metric02[0],1,&wsp0[0],1,&metric12[0],1,&wsp1[0],1,&wsp5[0],1);
2118  Vmath::Vvtvp (nqtot,&metric22[0],1,&wsp2[0],1,&wsp5[0],1,&wsp5[0],1);
2119 
2120  // outarray = m = (D_xi1 * B)^T * k
2121  // wsp1 = n = (D_xi2 * B)^T * l
2122  IProductWRTBase_SumFacKernel(dbase0,base1,base2,wsp3,outarray,wsp0,false,true,true);
2123  IProductWRTBase_SumFacKernel(base0,dbase1,base2,wsp4,wsp2, wsp0,true,false,true);
2124  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
2125  IProductWRTBase_SumFacKernel(base0,base1,dbase2,wsp5,wsp2, wsp0,true,true,false);
2126  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
2127  }
2128 
2129 
2130  void HexExp::v_ComputeLaplacianMetric()
2131  {
2132  if (m_metrics.count(MetricQuadrature) == 0)
2133  {
2134  ComputeQuadratureMetric();
2135  }
2136 
2137  const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
2138  const unsigned int nqtot = GetTotPoints();
2139  const unsigned int dim = 3;
2140  const MetricType m[3][3] = { {MetricLaplacian00, MetricLaplacian01, MetricLaplacian02},
2141  {MetricLaplacian01, MetricLaplacian11, MetricLaplacian12},
2142  {MetricLaplacian02, MetricLaplacian12, MetricLaplacian22}
2143  };
2144 
2145  for (unsigned int i = 0; i < dim; ++i)
2146  {
2147  for (unsigned int j = i; j < dim; ++j)
2148  {
2149  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
2150  const Array<TwoD, const NekDouble> &gmat =
2151  m_metricinfo->GetGmat(GetPointsKeys());
2152  if (type == SpatialDomains::eDeformed)
2153  {
2154  Vmath::Vcopy(nqtot, &gmat[i*dim+j][0], 1,
2155  &m_metrics[m[i][j]][0], 1);
2156  }
2157  else
2158  {
2159  Vmath::Fill(nqtot, gmat[i*dim+j][0],
2160  &m_metrics[m[i][j]][0], 1);
2161  }
2162  MultiplyByQuadratureMetric(m_metrics[m[i][j]],
2163  m_metrics[m[i][j]]);
2164 
2165  }
2166  }
2167  }
2168 
2169  }//end of namespace
2170 }//end of namespace
Generated on Sun Mar 15 2015 00:11:32 for Nektar++ by   doxygen 1.8.1.2