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  bool HexExp::v_GetFaceDGForwards(const int i) const
653  {
654  StdRegions::Orientation fo = GetGeom3D()->GetForient(i);
655 
656  return fo == StdRegions::eDir1FwdDir1_Dir2FwdDir2 ||
657  fo == StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
658  fo == StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
659  fo == StdRegions::eDir1FwdDir2_Dir2BwdDir1;
660  }
661 
662 
663  void HexExp::v_GetTracePhysVals(
664  const int face,
665  const StdRegions::StdExpansionSharedPtr &FaceExp,
666  const Array<OneD, const NekDouble> &inarray,
667  Array<OneD, NekDouble> &outarray,
668  StdRegions::Orientation orient)
669  {
670  v_GetFacePhysVals(face,FaceExp,inarray,outarray,orient);
671  }
672 
673 
674  ///Returns the physical values at the quadrature points of a face
675  void HexExp::v_GetFacePhysVals(
676  const int face,
677  const StdRegions::StdExpansionSharedPtr &FaceExp,
678  const Array<OneD, const NekDouble> &inarray,
679  Array<OneD, NekDouble> &outarray,
680  StdRegions::Orientation orient)
681  {
682  int nquad0 = m_base[0]->GetNumPoints();
683  int nquad1 = m_base[1]->GetNumPoints();
684  int nquad2 = m_base[2]->GetNumPoints();
685  Array<OneD, NekDouble> o_tmp(GetFaceNumPoints(face));
686 
687  if (orient == StdRegions::eNoOrientation)
688  {
689  orient = GetForient(face);
690  }
691 
692  switch(face)
693  {
694  case 0:
695  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
696  {
697  //Directions A and B positive
698  Vmath::Vcopy(nquad0*nquad1,&(inarray[0]),1,&(o_tmp[0]),1);
699  }
700  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
701  {
702  //Direction A negative and B positive
703  for (int j=0; j<nquad1; j++)
704  {
705  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+j*nquad0,-1,&(o_tmp[0])+(j*nquad0),1);
706  }
707  }
708  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
709  {
710  //Direction A positive and B negative
711  for (int j=0; j<nquad1; j++)
712  {
713  Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1-j),1,&(o_tmp[0])+(j*nquad0),1);
714  }
715  }
716  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
717  {
718  //Direction A negative and B negative
719  for(int j=0; j<nquad1; j++)
720  {
721  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1-j*nquad0),-1,&(o_tmp[0])+(j*nquad0),1);
722  }
723  }
724  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
725  {
726  //Transposed, Direction A and B positive
727  for (int i=0; i<nquad0; i++)
728  {
729  Vmath::Vcopy(nquad1,&(inarray[0])+i,nquad0,&(o_tmp[0])+(i*nquad1),1);
730  }
731  }
732  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
733  {
734  //Transposed, Direction A negative and B positive
735  for (int i=0; i<nquad0; i++)
736  {
737  Vmath::Vcopy(nquad1,&(inarray[0])+i+nquad0*(nquad1-1),-nquad0,&(o_tmp[0])+(i*nquad1),1);
738  }
739  }
740  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
741  {
742  //Transposed, Direction A positive and B negative
743  for (int i=0; i<nquad0; i++)
744  {
745  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0-1-i),nquad0,&(o_tmp[0])+(i*nquad1),1);
746  }
747  }
748  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
749  {
750  //Transposed, Direction A and B negative
751  for (int i=0; i<nquad0; i++)
752  {
753  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0*nquad1-1-i),-nquad0,&(o_tmp[0])+(i*nquad1),1);
754  }
755  }
756  //interpolate
757  LibUtilities::Interp2D(m_base[0]->GetPointsKey(),
758  m_base[1]->GetPointsKey(), o_tmp,
759  FaceExp->GetBasis(0)->GetPointsKey(),
760  FaceExp->GetBasis(1)->GetPointsKey(),
761  outarray);
762  break;
763  case 1:
764  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
765  {
766  //Direction A and B positive
767  for (int k=0; k<nquad2; k++)
768  {
769  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1*k),
770  1,&(o_tmp[0])+(k*nquad0),1);
771  }
772  }
773  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
774  {
775  //Direction A negative and B positive
776  for (int k=0; k<nquad2; k++)
777  {
778  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+(nquad0*nquad1*k),
779  -1,&(o_tmp[0])+(k*nquad0),1);
780  }
781  }
782  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
783  {
784  //Direction A positive and B negative
785  for (int k=0; k<nquad2; k++)
786  {
787  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1*(nquad2-1-k)),
788  1,&(o_tmp[0])+(k*nquad0),1);
789  }
790  }
791  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
792  {
793  //Direction A negative and B negative
794  for(int k=0; k<nquad2; k++)
795  {
796  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+(nquad0*nquad1*(nquad2-1-k)),
797  -1,&(o_tmp[0])+(k*nquad0),1);
798  }
799  }
800  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
801  {
802  //Transposed, Direction A and B positive
803  for (int i=0; i<nquad0; i++)
804  {
805  Vmath::Vcopy(nquad2,&(inarray[0])+i,nquad0*nquad1,
806  &(o_tmp[0])+(i*nquad2),1);
807  }
808  }
809  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
810  {
811  //Transposed, Direction A negative and B positive
812  for (int i=0; i<nquad0; i++)
813  {
814  Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*nquad1*(nquad2-1)+i,
815  -nquad0*nquad1,&(o_tmp[0])+(i*nquad2),1);
816  }
817  }
818  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
819  {
820  //Transposed, Direction A positive and B negative
821  for (int i=0; i<nquad0; i++)
822  {
823  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0-1-i),nquad0*nquad1,
824  &(o_tmp[0])+(i*nquad2),1);
825  }
826  }
827  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
828  {
829  //Transposed, Direction A and B negative
830  for (int i=0; i<nquad0; i++)
831  {
832  Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*nquad1*nquad2+(nquad0-1-i),
833  -nquad0*nquad1,&(o_tmp[0])+(i*nquad2),1);
834  }
835  }
836  //interpolate
837  LibUtilities::Interp2D(m_base[0]->GetPointsKey(),
838  m_base[2]->GetPointsKey(), o_tmp,
839  FaceExp->GetBasis(0)->GetPointsKey(),
840  FaceExp->GetBasis(1)->GetPointsKey(),
841  outarray);
842  break;
843  case 2:
844  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
845  {
846  //Directions A and B positive
847  Vmath::Vcopy(nquad0*nquad1,&(inarray[0])+(nquad0-1),
848  nquad0,&(o_tmp[0]),1);
849  }
850  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
851  {
852  //Direction A negative and B positive
853  for (int k=0; k<nquad2; k++)
854  {
855  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1)+(k*nquad0*nquad1),
856  -nquad0,&(o_tmp[0])+(k*nquad0),1);
857  }
858  }
859  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
860  {
861  //Direction A positive and B negative
862  for (int k=0; k<nquad2; k++)
863  {
864  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+(nquad0*nquad1*(nquad2-1-k)),
865  nquad0,&(o_tmp[0])+(k*nquad0),1);
866  }
867  }
868  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
869  {
870  //Direction A negative and B negative
871  for (int k=0; k<nquad2; k++)
872  {
873  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1)+(nquad0*nquad1*(nquad2-1-k)),
874  -nquad0,&(o_tmp[0])+(k*nquad0),1);
875  }
876  }
877  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
878  {
879  //Transposed, Direction A and B positive
880  for (int j=0; j<nquad1; j++)
881  {
882  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0-1)+(j*nquad0),
883  nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
884  }
885  }
886  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
887  {
888  //Transposed, Direction A negative and B positive
889  for (int j=0; j<nquad0; j++)
890  {
891  Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*nquad1*(nquad2-1)+nquad0+j*nquad0,
892  -nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
893  }
894  }
895  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
896  {
897  //Transposed, Direction A positive and B negative
898  for (int j=0; j<nquad0; j++)
899  {
900  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*nquad1-1-j*nquad0),
901  nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
902  }
903  }
904  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
905  {
906  //Transposed, Direction A and B negative
907  for (int j=0; j<nquad0; j++)
908  {
909  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*nquad1*nquad2-1-j*nquad0),
910  -nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
911  }
912  }
913  //interpolate
914  LibUtilities::Interp2D(m_base[1]->GetPointsKey(),
915  m_base[2]->GetPointsKey(), o_tmp,
916  FaceExp->GetBasis(0)->GetPointsKey(),
917  FaceExp->GetBasis(1)->GetPointsKey(),
918  outarray);
919  break;
920  case 3:
921  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
922  {
923  //Directions A and B positive
924  for (int k=0; k<nquad2; k++)
925  {
926  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*(nquad1-1))+(k*nquad0*nquad1),
927  1,&(o_tmp[0])+(k*nquad0),1);
928  }
929  }
930  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
931  {
932  //Direction A negative 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,&(o_tmp[0])+(k*nquad0),1);
937  }
938  }
939  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
940  {
941  //Direction A positive and B negative
942  for (int k=0; k<nquad2; k++)
943  {
944  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*(nquad1-1))+(nquad0*nquad1*(nquad2-1-k)),
945  1,&(o_tmp[0])+(k*nquad0),1);
946  }
947  }
948  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
949  {
950  //Direction A negative 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,&(o_tmp[0])+(k*nquad0),1);
955  }
956  }
957  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
958  {
959  //Transposed, Direction A and B positive
960  for (int i=0; i<nquad0; i++)
961  {
962  Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*(nquad1-1)+i,nquad0*nquad1,
963  &(o_tmp[0])+(i*nquad2),1);
964  }
965  }
966  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
967  {
968  //Transposed, Direction A negative and B positive
969  for (int i=0; i<nquad0; i++)
970  {
971  Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*(nquad1*nquad2-1)+i,-nquad0*nquad1,
972  &(o_tmp[0])+(i*nquad2),1);
973  }
974  }
975  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
976  {
977  //Transposed, Direction A positive and B negative
978  for (int i=0; i<nquad0; i++)
979  {
980  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*nquad1-1-i),nquad0*nquad1,
981  &(o_tmp[0])+(i*nquad2),1);
982  }
983  }
984  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
985  {
986  //Transposed, Direction A and B negative
987  for (int i=0; i<nquad0; i++)
988  {
989  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*nquad1*nquad2-1-i),-nquad0*nquad1,
990  &(o_tmp[0])+(i*nquad2),1);
991  }
992  }
993  //interpolate
994  LibUtilities::Interp2D(m_base[0]->GetPointsKey(),
995  m_base[2]->GetPointsKey(), o_tmp,
996  FaceExp->GetBasis(0)->GetPointsKey(),
997  FaceExp->GetBasis(1)->GetPointsKey(),
998  outarray);
999  break;
1000  case 4:
1001  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
1002  {
1003  //Directions A and B positive
1004  Vmath::Vcopy(nquad0*nquad1,&(inarray[0]),nquad0,&(o_tmp[0]),1);
1005  }
1006  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
1007  {
1008  //Direction A negative and B positive
1009  for (int k=0; k<nquad2; k++)
1010  {
1011  Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1)+(k*nquad0*nquad1),
1012  -nquad0,&(o_tmp[0])+(k*nquad0),1);
1013  }
1014  }
1015  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
1016  {
1017  //Direction A positive and B negative
1018  for (int k=0; k<nquad2; k++)
1019  {
1020  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1*(nquad2-1-k)),
1021  nquad0,&(o_tmp[0])+(k*nquad0),1);
1022  }
1023  }
1024  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
1025  {
1026  //Direction A negative and B negative
1027  for (int k=0; k<nquad2; k++)
1028  {
1029  Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1)+(nquad0*nquad1*(nquad2-1-k)),
1030  -nquad0,&(o_tmp[0])+(k*nquad0),1);
1031  }
1032  }
1033  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
1034  {
1035  //Transposed, Direction A and B positive
1036  for (int j=0; j<nquad0; j++)
1037  {
1038  Vmath::Vcopy(nquad2,&(inarray[0])+j*nquad0,nquad0*nquad1,
1039  &(o_tmp[0])+(j*nquad2),1);
1040  }
1041  }
1042  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
1043  {
1044  //Transposed, Direction A negative and B positive
1045  for (int j=0; j<nquad0; j++)
1046  {
1047  Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*nquad1*(nquad2-1)+j*nquad0,
1048  -nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
1049  }
1050  }
1051  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
1052  {
1053  //Transposed, Direction A positive and B negative
1054  for (int j=0; j<nquad0; j++)
1055  {
1056  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*(nquad1-1)-j*nquad0),
1057  nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
1058  }
1059  }
1060  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
1061  {
1062  //Transposed, Direction A and B negative
1063  for (int j=0; j<nquad0; j++)
1064  {
1065  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*(nquad1*nquad2-1)-j*nquad0),
1066  -nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
1067  }
1068  }
1069  //interpolate
1070  LibUtilities::Interp2D(m_base[1]->GetPointsKey(),
1071  m_base[2]->GetPointsKey(), o_tmp,
1072  FaceExp->GetBasis(0)->GetPointsKey(),
1073  FaceExp->GetBasis(1)->GetPointsKey(),
1074  outarray);
1075  break;
1076  case 5:
1077  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
1078  {
1079  //Directions A and B positive
1080  Vmath::Vcopy(nquad0*nquad1,&(inarray[0])+nquad0*nquad1*(nquad2-1),1,&(o_tmp[0]),1);
1081  }
1082  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
1083  {
1084  //Direction A negative and B positive
1085  for (int j=0; j<nquad1; j++)
1086  {
1087  Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*nquad1*(nquad2-1)+(nquad0-1+j*nquad0),
1088  -1,&(o_tmp[0])+(j*nquad0),1);
1089  }
1090  }
1091  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
1092  {
1093  //Direction A positive and B negative
1094  for (int j=0; j<nquad1; j++)
1095  {
1096  Vmath::Vcopy(nquad0,&(inarray[0])+((nquad0*nquad1*nquad2-1)-(nquad0-1)-j*nquad0),
1097  1,&(o_tmp[0])+(j*nquad0),1);
1098  }
1099  }
1100  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
1101  {
1102  //Direction A negative and B negative
1103  for (int j=0; j<nquad1; j++)
1104  {
1105  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1*nquad2-1-j*nquad0),
1106  -1,&(o_tmp[0])+(j*nquad0),1);
1107  }
1108  }
1109  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
1110  {
1111  //Transposed, Direction A and B positive
1112  for (int i=0; i<nquad0; i++)
1113  {
1114  Vmath::Vcopy(nquad1,&(inarray[0])+nquad0*nquad1*(nquad2-1)+i,nquad0,
1115  &(o_tmp[0])+(i*nquad1),1);
1116  }
1117  }
1118  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
1119  {
1120  //Transposed, Direction A negative and B positive
1121  for (int i=0; i<nquad0; i++)
1122  {
1123  Vmath::Vcopy(nquad1,&(inarray[0])+nquad0*(nquad1*nquad2-1)+i,-nquad0,
1124  &(o_tmp[0])+(i*nquad1),1);
1125  }
1126  }
1127  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
1128  {
1129  //Transposed, Direction A positive and B negative
1130  for (int i=0; i<nquad0; i++)
1131  {
1132  Vmath::Vcopy(nquad1,&(inarray[0])+nquad0*nquad1*(nquad2-1)+(nquad0-1-i),
1133  nquad0,&(o_tmp[0])+(i*nquad1),1);
1134  }
1135  }
1136  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
1137  {
1138  //Transposed, Direction A and B negative
1139  for (int i=0; i<nquad0; i++)
1140  {
1141  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0*nquad1*nquad2-1-i),-nquad0,
1142  &(o_tmp[0])+(i*nquad1),1);
1143  }
1144  }
1145  //interpolate
1146  LibUtilities::Interp2D(m_base[0]->GetPointsKey(),
1147  m_base[1]->GetPointsKey(), o_tmp,
1148  FaceExp->GetBasis(0)->GetPointsKey(),
1149  FaceExp->GetBasis(1)->GetPointsKey(),
1150  outarray);
1151  break;
1152  default:
1153  ASSERTL0(false,"face value (> 5) is out of range");
1154  break;
1155  }
1156  }
1157 
1158 
1159  void HexExp::v_ComputeFaceNormal(const int face)
1160  {
1161  int i;
1162  const SpatialDomains::GeomFactorsSharedPtr & geomFactors =
1163  GetGeom()->GetMetricInfo();
1164  SpatialDomains::GeomType type = geomFactors->GetGtype();
1165  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1166  const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys);
1167  const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys);
1168 
1169  LibUtilities::BasisKey tobasis0 = DetFaceBasisKey(face,0);
1170  LibUtilities::BasisKey tobasis1 = DetFaceBasisKey(face,1);
1171 
1172  // Number of quadrature points in face expansion.
1173  int nq_face = tobasis0.GetNumPoints()*tobasis1.GetNumPoints();
1174 
1175  int vCoordDim = GetCoordim();
1176 
1177  m_faceNormals[face] = Array<OneD, Array<OneD, NekDouble> >(vCoordDim);
1178  Array<OneD, Array<OneD, NekDouble> > &normal = m_faceNormals[face];
1179  for (i = 0; i < vCoordDim; ++i)
1180  {
1181  normal[i] = Array<OneD, NekDouble>(nq_face);
1182  }
1183  // Regular geometry case
1184  if((type == SpatialDomains::eRegular)||(type == SpatialDomains::eMovingRegular))
1185  {
1186  NekDouble fac;
1187  // Set up normals
1188  switch(face)
1189  {
1190  case 0:
1191  for(i = 0; i < vCoordDim; ++i)
1192  {
1193  normal[i][0] = -df[3*i+2][0];
1194  }
1195  break;
1196  case 1:
1197  for(i = 0; i < vCoordDim; ++i)
1198  {
1199  normal[i][0] = -df[3*i+1][0];
1200  }
1201  break;
1202  case 2:
1203  for(i = 0; i < vCoordDim; ++i)
1204  {
1205  normal[i][0] = df[3*i][0];
1206  }
1207  break;
1208  case 3:
1209  for(i = 0; i < vCoordDim; ++i)
1210  {
1211  normal[i][0] = df[3*i+1][0];
1212  }
1213  break;
1214  case 4:
1215  for(i = 0; i < vCoordDim; ++i)
1216  {
1217  normal[i][0] = -df[3*i][0];
1218  }
1219  break;
1220  case 5:
1221  for(i = 0; i < vCoordDim; ++i)
1222  {
1223  normal[i][0] = df[3*i+2][0];
1224  }
1225  break;
1226  default:
1227  ASSERTL0(false,"face is out of range (edge < 5)");
1228  }
1229 
1230  // normalise
1231  fac = 0.0;
1232  for(i =0 ; i < vCoordDim; ++i)
1233  {
1234  fac += normal[i][0]*normal[i][0];
1235  }
1236  fac = 1.0/sqrt(fac);
1237  for (i = 0; i < vCoordDim; ++i)
1238  {
1239  Vmath::Fill(nq_face,fac*normal[i][0],normal[i],1);
1240  }
1241 
1242  }
1243  else // Set up deformed normals
1244  {
1245  int j, k;
1246 
1247  int nqe0 = m_base[0]->GetNumPoints();
1248  int nqe1 = m_base[1]->GetNumPoints();
1249  int nqe2 = m_base[2]->GetNumPoints();
1250  int nqe01 = nqe0*nqe1;
1251  int nqe02 = nqe0*nqe2;
1252  int nqe12 = nqe1*nqe2;
1253 
1254  int nqe;
1255  if (face == 0 || face == 5)
1256  {
1257  nqe = nqe01;
1258  }
1259  else if (face == 1 || face == 3)
1260  {
1261  nqe = nqe02;
1262  }
1263  else
1264  {
1265  nqe = nqe12;
1266  }
1267 
1268  LibUtilities::PointsKey points0;
1269  LibUtilities::PointsKey points1;
1270 
1271  Array<OneD, NekDouble> normals(vCoordDim*nqe,0.0);
1272 
1273  // Extract Jacobian along face and recover local
1274  // derivates (dx/dr) for polynomial interpolation by
1275  // multiplying m_gmat by jacobian
1276  switch(face)
1277  {
1278  case 0:
1279  for(j = 0; j < nqe; ++j)
1280  {
1281  normals[j] = -df[2][j]*jac[j];
1282  normals[nqe+j] = -df[5][j]*jac[j];
1283  normals[2*nqe+j] = -df[8][j]*jac[j];
1284  }
1285 
1286  points0 = ptsKeys[0];
1287  points1 = ptsKeys[1];
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  normals[j+k*nqe0] = -df[1][idx]*jac[idx];
1296  normals[nqe+j+k*nqe0] = -df[4][idx]*jac[idx];
1297  normals[2*nqe+j+k*nqe0] = -df[7][idx]*jac[idx];
1298  }
1299  }
1300  points0 = ptsKeys[0];
1301  points1 = ptsKeys[2];
1302  break;
1303  case 2:
1304  for (j = 0; j < nqe1; ++j)
1305  {
1306  for(k = 0; k < nqe2; ++k)
1307  {
1308  int idx = nqe0-1+nqe0*j+nqe01*k;
1309  normals[j+k*nqe0] = df[0][idx]*jac[idx];
1310  normals[nqe+j+k*nqe0] = df[3][idx]*jac[idx];
1311  normals[2*nqe+j+k*nqe0] = df[6][idx]*jac[idx];
1312  }
1313  }
1314  points0 = ptsKeys[1];
1315  points1 = ptsKeys[2];
1316  break;
1317  case 3:
1318  for (j = 0; j < nqe0; ++j)
1319  {
1320  for(k = 0; k < nqe2; ++k)
1321  {
1322  int idx = nqe0*(nqe1-1)+j+nqe01*k;
1323  normals[j+k*nqe0] = df[1][idx]*jac[idx];
1324  normals[nqe+j+k*nqe0] = df[4][idx]*jac[idx];
1325  normals[2*nqe+j+k*nqe0] = df[7][idx]*jac[idx];
1326  }
1327  }
1328  points0 = ptsKeys[0];
1329  points1 = ptsKeys[2];
1330  break;
1331  case 4:
1332  for (j = 0; j < nqe0; ++j)
1333  {
1334  for(k = 0; k < nqe2; ++k)
1335  {
1336  int idx = j*nqe0+nqe01*k;
1337  normals[j+k*nqe0] = -df[0][idx]*jac[idx];
1338  normals[nqe+j+k*nqe0] = -df[3][idx]*jac[idx];
1339  normals[2*nqe+j+k*nqe0] = -df[6][idx]*jac[idx];
1340  }
1341  }
1342  points0 = ptsKeys[1];
1343  points1 = ptsKeys[2];
1344  break;
1345  case 5:
1346  for (j = 0; j < nqe01; ++j)
1347  {
1348  int idx = j+nqe01*(nqe2-1);
1349  normals[j] = df[2][idx]*jac[idx];
1350  normals[nqe+j] = df[5][idx]*jac[idx];
1351  normals[2*nqe+j] = df[8][idx]*jac[idx];
1352  }
1353  points0 = ptsKeys[0];
1354  points1 = ptsKeys[1];
1355  break;
1356  default:
1357  ASSERTL0(false,"face is out of range (face < 5)");
1358  }
1359 
1360  Array<OneD, NekDouble> work (nq_face, 0.0);
1361  // Interpolate Jacobian and invert
1362  LibUtilities::Interp2D(points0, points1, jac,
1363  tobasis0.GetPointsKey(),
1364  tobasis1.GetPointsKey(),
1365  work);
1366 
1367  Vmath::Sdiv(nq_face,1.0,&work[0],1,&work[0],1);
1368 
1369  // interpolate
1370  for(i = 0; i < GetCoordim(); ++i)
1371  {
1372  LibUtilities::Interp2D(points0, points1,
1373  &normals[i*nqe],
1374  tobasis0.GetPointsKey(),
1375  tobasis1.GetPointsKey(),
1376  &normal[i][0]);
1377  Vmath::Vmul(nq_face,work,1,normal[i],1,normal[i],1);
1378  }
1379 
1380  //normalise normal vectors
1381  Vmath::Zero(nq_face,work,1);
1382  for(i = 0; i < GetCoordim(); ++i)
1383  {
1384  Vmath::Vvtvp(nq_face,normal[i],1, normal[i],1,work,1,work,1);
1385  }
1386 
1387  Vmath::Vsqrt(nq_face,work,1,work,1);
1388  Vmath::Sdiv(nq_face,1.0,work,1,work,1);
1389 
1390  for(i = 0; i < GetCoordim(); ++i)
1391  {
1392  Vmath::Vmul(nq_face,normal[i],1,work,1,normal[i],1);
1393  }
1394  }
1395  }
1396 
1397  //-----------------------------
1398  // Operator creation functions
1399  //-----------------------------
1400  void HexExp::v_MassMatrixOp(
1401  const Array<OneD, const NekDouble> &inarray,
1402  Array<OneD,NekDouble> &outarray,
1403  const StdRegions::StdMatrixKey &mkey)
1404  {
1405  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
1406  }
1407 
1408  void HexExp::v_LaplacianMatrixOp(
1409  const Array<OneD, const NekDouble> &inarray,
1410  Array<OneD,NekDouble> &outarray,
1411  const StdRegions::StdMatrixKey &mkey)
1412  {
1413  HexExp::v_LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1414  }
1415 
1416  void HexExp::v_LaplacianMatrixOp(
1417  const int k1,
1418  const int k2,
1419  const Array<OneD, const NekDouble> &inarray,
1420  Array<OneD,NekDouble> &outarray,
1421  const StdRegions::StdMatrixKey &mkey)
1422  {
1423  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,
1424  mkey);
1425  }
1426 
1427  void HexExp::v_WeakDerivMatrixOp(
1428  const int i,
1429  const Array<OneD, const NekDouble> &inarray,
1430  Array<OneD,NekDouble> &outarray,
1431  const StdRegions::StdMatrixKey &mkey)
1432  {
1433  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
1434  }
1435 
1436  void HexExp::v_WeakDirectionalDerivMatrixOp(
1437  const Array<OneD, const NekDouble> &inarray,
1438  Array<OneD,NekDouble> &outarray,
1439  const StdRegions::StdMatrixKey &mkey)
1440  {
1441  StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray,
1442  outarray,mkey);
1443  }
1444 
1445  void HexExp::v_MassLevelCurvatureMatrixOp(
1446  const Array<OneD, const NekDouble> &inarray,
1447  Array<OneD,NekDouble> &outarray,
1448  const StdRegions::StdMatrixKey &mkey)
1449  {
1450  StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray,
1451  outarray,mkey);
1452  }
1453 
1454  void HexExp::v_HelmholtzMatrixOp(
1455  const Array<OneD, const NekDouble> &inarray,
1456  Array<OneD,NekDouble> &outarray,
1457  const StdRegions::StdMatrixKey &mkey)
1458  {
1459  HexExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1460  }
1461 
1462 
1463  void HexExp::v_GeneralMatrixOp_MatOp(
1464  const Array<OneD, const NekDouble> &inarray,
1465  Array<OneD,NekDouble> &outarray,
1466  const StdRegions::StdMatrixKey &mkey)
1467  {
1468  //int nConsts = mkey.GetNconstants();
1469  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1470 
1471 // switch(nConsts)
1472 // {
1473 // case 0:
1474 // {
1475 // mat = GetLocMatrix(mkey.GetMatrixType());
1476 // }
1477 // break;
1478 // case 1:
1479 // {
1480 // mat = GetLocMatrix(mkey.GetMatrixType(),mkey.GetConstant(0));
1481 // }
1482 // break;
1483 // case 2:
1484 // {
1485 // mat = GetLocMatrix(mkey.GetMatrixType(),mkey.GetConstant(0),mkey.GetConstant(1));
1486 // }
1487 // break;
1488 //
1489 // default:
1490 // {
1491 // NEKERROR(ErrorUtil::efatal, "Unknown number of constants");
1492 // }
1493 // break;
1494 // }
1495 
1496  if(inarray.get() == outarray.get())
1497  {
1498  Array<OneD,NekDouble> tmp(m_ncoeffs);
1499  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1500 
1501  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1502  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1503  }
1504  else
1505  {
1506  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1507  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1508  }
1509  }
1510 
1511  /**
1512  * This function is used to compute exactly the advective numerical flux
1513  * on the interface of two elements with different expansions, hence an
1514  * appropriate number of Gauss points has to be used. The number of
1515  * Gauss points has to be equal to the number used by the highest
1516  * polynomial degree of the two adjacent elements
1517  *
1518  * @param numMin Is the reduced polynomial order
1519  * @param inarray Input array of coefficients
1520  * @param dumpVar Output array of reduced coefficients.
1521  */
1522  void HexExp::v_ReduceOrderCoeffs(
1523  int numMin,
1524  const Array<OneD, const NekDouble> &inarray,
1525  Array<OneD, NekDouble> &outarray)
1526  {
1527  int n_coeffs = inarray.num_elements();
1528  int nmodes0 = m_base[0]->GetNumModes();
1529  int nmodes1 = m_base[1]->GetNumModes();
1530  int nmodes2 = m_base[2]->GetNumModes();
1531  int numMax = nmodes0;
1532 
1533  Array<OneD, NekDouble> coeff (n_coeffs);
1534  Array<OneD, NekDouble> coeff_tmp1(nmodes0*nmodes1, 0.0);
1535  Array<OneD, NekDouble> coeff_tmp2(n_coeffs, 0.0);
1536  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
1537 
1538  Vmath::Vcopy(n_coeffs,inarray,1,coeff_tmp2,1);
1539 
1540  const LibUtilities::PointsKey Pkey0(
1541  nmodes0, LibUtilities::eGaussLobattoLegendre);
1542  const LibUtilities::PointsKey Pkey1(
1543  nmodes1, LibUtilities::eGaussLobattoLegendre);
1544  const LibUtilities::PointsKey Pkey2(
1545  nmodes2, LibUtilities::eGaussLobattoLegendre);
1546 
1547  LibUtilities::BasisKey b0(
1548  m_base[0]->GetBasisType(), nmodes0, Pkey0);
1549  LibUtilities::BasisKey b1(
1550  m_base[1]->GetBasisType(), nmodes1, Pkey1);
1551  LibUtilities::BasisKey b2(
1552  m_base[2]->GetBasisType(), nmodes2, Pkey2);
1553  LibUtilities::BasisKey bortho0(
1554  LibUtilities::eOrtho_A, nmodes0, Pkey0);
1555  LibUtilities::BasisKey bortho1(
1556  LibUtilities::eOrtho_A, nmodes1, Pkey1);
1557  LibUtilities::BasisKey bortho2(
1558  LibUtilities::eOrtho_A, nmodes2, Pkey2);
1559 
1560  LibUtilities::InterpCoeff3D(
1561  b0, b1, b2, coeff_tmp2,
1562  bortho0, bortho1, bortho2, coeff);
1563 
1564  Vmath::Zero(n_coeffs, coeff_tmp2, 1);
1565 
1566  int cnt = 0, cnt2 = 0;
1567 
1568  for (int u = 0; u < numMin+1; ++u)
1569  {
1570  for (int i = 0; i < numMin; ++i)
1571  {
1572  Vmath::Vcopy(numMin,
1573  tmp = coeff+cnt+cnt2,1,
1574  tmp2 = coeff_tmp1+cnt,1);
1575 
1576  cnt = i*numMax;
1577  }
1578 
1579  Vmath::Vcopy(nmodes0*nmodes1,
1580  tmp3 = coeff_tmp1,1,
1581  tmp4 = coeff_tmp2+cnt2,1);
1582 
1583  cnt2 = u*nmodes0*nmodes1;
1584  }
1585 
1586  LibUtilities::InterpCoeff3D(
1587  bortho0, bortho1, bortho2, coeff_tmp2,
1588  b0, b1, b2, outarray);
1589  }
1590 
1591  //-----------------------------
1592  // Matrix creation functions
1593  //-----------------------------
1594  DNekMatSharedPtr HexExp::v_GenMatrix(
1595  const StdRegions::StdMatrixKey &mkey)
1596  {
1597  DNekMatSharedPtr returnval;
1598 
1599  switch(mkey.GetMatrixType())
1600  {
1601  case StdRegions::eHybridDGHelmholtz:
1602  case StdRegions::eHybridDGLamToU:
1603  case StdRegions::eHybridDGLamToQ0:
1604  case StdRegions::eHybridDGLamToQ1:
1605  case StdRegions::eHybridDGLamToQ2:
1606  case StdRegions::eHybridDGHelmBndLam:
1607  case StdRegions::eInvLaplacianWithUnityMean:
1608  returnval = Expansion3D::v_GenMatrix(mkey);
1609  break;
1610  default:
1611  returnval = StdHexExp::v_GenMatrix(mkey);
1612  }
1613 
1614  return returnval;
1615  }
1616 
1617 
1618  DNekMatSharedPtr HexExp::v_CreateStdMatrix(
1619  const StdRegions::StdMatrixKey &mkey)
1620  {
1621  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1622  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1623  LibUtilities::BasisKey bkey2 = m_base[2]->GetBasisKey();
1624 
1625  StdRegions::StdHexExpSharedPtr tmp = MemoryManager<StdHexExp>
1626  ::AllocateSharedPtr(bkey0,bkey1,bkey2);
1627 
1628  return tmp->GetStdMatrix(mkey);
1629  }
1630 
1631 
1632  DNekScalMatSharedPtr HexExp::CreateMatrix(const MatrixKey &mkey)
1633  {
1634  DNekScalMatSharedPtr returnval;
1635  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1636 
1637  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,
1638  "Geometric information is not set up");
1639 
1640  switch(mkey.GetMatrixType())
1641  {
1642  case StdRegions::eMass:
1643  {
1644  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1645  mkey.GetNVarCoeff())
1646  {
1647  NekDouble one = 1.0;
1648  DNekMatSharedPtr mat = GenMatrix(mkey);
1649  returnval = MemoryManager<DNekScalMat>
1650  ::AllocateSharedPtr(one,mat);
1651  }
1652  else
1653  {
1654  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1655  DNekMatSharedPtr mat
1656  = GetStdMatrix(mkey);
1657  returnval = MemoryManager<DNekScalMat>
1658  ::AllocateSharedPtr(jac,mat);
1659  }
1660  }
1661  break;
1662  case StdRegions::eInvMass:
1663  {
1664  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1665  {
1666  NekDouble one = 1.0;
1667  StdRegions::StdMatrixKey masskey(StdRegions::eMass,
1668  DetShapeType(), *this);
1669  DNekMatSharedPtr mat = GenMatrix(masskey);
1670  mat->Invert();
1671 
1672  returnval = MemoryManager<DNekScalMat>
1673  ::AllocateSharedPtr(one,mat);
1674  }
1675  else
1676  {
1677  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1678  DNekMatSharedPtr mat
1679  = GetStdMatrix(mkey);
1680  returnval = MemoryManager<DNekScalMat>
1681  ::AllocateSharedPtr(fac,mat);
1682  }
1683  }
1684  break;
1685  case StdRegions::eWeakDeriv0:
1686  case StdRegions::eWeakDeriv1:
1687  case StdRegions::eWeakDeriv2:
1688  {
1689  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1690  mkey.GetNVarCoeff())
1691  {
1692  NekDouble one = 1.0;
1693  DNekMatSharedPtr mat = GenMatrix(mkey);
1694 
1695  returnval = MemoryManager<DNekScalMat>
1696  ::AllocateSharedPtr(one,mat);
1697  }
1698  else
1699  {
1700  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1701  Array<TwoD, const NekDouble> df
1702  = m_metricinfo->GetDerivFactors(ptsKeys);
1703  int dir = 0;
1704 
1705  switch(mkey.GetMatrixType())
1706  {
1707  case StdRegions::eWeakDeriv0:
1708  dir = 0;
1709  break;
1710  case StdRegions::eWeakDeriv1:
1711  dir = 1;
1712  break;
1713  case StdRegions::eWeakDeriv2:
1714  dir = 2;
1715  break;
1716  default:
1717  break;
1718  }
1719 
1720  MatrixKey deriv0key(StdRegions::eWeakDeriv0,
1721  mkey.GetShapeType(), *this);
1722  MatrixKey deriv1key(StdRegions::eWeakDeriv1,
1723  mkey.GetShapeType(), *this);
1724  MatrixKey deriv2key(StdRegions::eWeakDeriv2,
1725  mkey.GetShapeType(), *this);
1726 
1727  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
1728  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
1729  DNekMat &deriv2 = *GetStdMatrix(deriv2key);
1730 
1731  int rows = deriv0.GetRows();
1732  int cols = deriv1.GetColumns();
1733 
1734  DNekMatSharedPtr WeakDeriv = MemoryManager<DNekMat>
1735  ::AllocateSharedPtr(rows,cols);
1736 
1737  (*WeakDeriv) = df[3*dir ][0]*deriv0
1738  + df[3*dir+1][0]*deriv1
1739  + df[3*dir+2][0]*deriv2;
1740 
1741  returnval = MemoryManager<DNekScalMat>
1742  ::AllocateSharedPtr(jac,WeakDeriv);
1743  }
1744  }
1745  break;
1746  case StdRegions::eLaplacian:
1747  {
1748  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1749  mkey.GetNVarCoeff()||
1750  mkey.ConstFactorExists(
1751  StdRegions::eFactorSVVCutoffRatio))
1752  {
1753  NekDouble one = 1.0;
1754  DNekMatSharedPtr mat = GenMatrix(mkey);
1755 
1756  returnval = MemoryManager<DNekScalMat>
1757  ::AllocateSharedPtr(one,mat);
1758  }
1759  else
1760  {
1761  MatrixKey lap00key(StdRegions::eLaplacian00,
1762  mkey.GetShapeType(), *this);
1763  MatrixKey lap01key(StdRegions::eLaplacian01,
1764  mkey.GetShapeType(), *this);
1765  MatrixKey lap02key(StdRegions::eLaplacian02,
1766  mkey.GetShapeType(), *this);
1767  MatrixKey lap11key(StdRegions::eLaplacian11,
1768  mkey.GetShapeType(), *this);
1769  MatrixKey lap12key(StdRegions::eLaplacian12,
1770  mkey.GetShapeType(), *this);
1771  MatrixKey lap22key(StdRegions::eLaplacian22,
1772  mkey.GetShapeType(), *this);
1773 
1774  DNekMat &lap00 = *GetStdMatrix(lap00key);
1775  DNekMat &lap01 = *GetStdMatrix(lap01key);
1776  DNekMat &lap02 = *GetStdMatrix(lap02key);
1777  DNekMat &lap11 = *GetStdMatrix(lap11key);
1778  DNekMat &lap12 = *GetStdMatrix(lap12key);
1779  DNekMat &lap22 = *GetStdMatrix(lap22key);
1780 
1781  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1782  Array<TwoD, const NekDouble> gmat
1783  = m_metricinfo->GetGmat(ptsKeys);
1784 
1785  int rows = lap00.GetRows();
1786  int cols = lap00.GetColumns();
1787 
1788  DNekMatSharedPtr lap = MemoryManager<DNekMat>
1789  ::AllocateSharedPtr(rows,cols);
1790 
1791  (*lap) = gmat[0][0]*lap00
1792  + gmat[4][0]*lap11
1793  + gmat[8][0]*lap22
1794  + gmat[3][0]*(lap01 + Transpose(lap01))
1795  + gmat[6][0]*(lap02 + Transpose(lap02))
1796  + gmat[7][0]*(lap12 + Transpose(lap12));
1797 
1798  returnval = MemoryManager<DNekScalMat>
1799  ::AllocateSharedPtr(jac,lap);
1800  }
1801  }
1802  break;
1803  case StdRegions::eHelmholtz:
1804  {
1805  NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
1806  MatrixKey masskey(StdRegions::eMass,
1807  mkey.GetShapeType(), *this);
1808  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1809  MatrixKey lapkey(StdRegions::eLaplacian,
1810  mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1811  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1812 
1813  int rows = LapMat.GetRows();
1814  int cols = LapMat.GetColumns();
1815 
1816  DNekMatSharedPtr helm = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1817 
1818  NekDouble one = 1.0;
1819  (*helm) = LapMat + lambda*MassMat;
1820 
1821  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
1822  }
1823  break;
1824  case StdRegions::eIProductWRTBase:
1825  {
1826  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1827  {
1828  NekDouble one = 1.0;
1829  DNekMatSharedPtr mat = GenMatrix(mkey);
1830  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1831  }
1832  else
1833  {
1834  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1835  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1836  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1837  }
1838  }
1839  break;
1840  case StdRegions::eHybridDGHelmholtz:
1841  case StdRegions::eHybridDGLamToU:
1842  case StdRegions::eHybridDGLamToQ0:
1843  case StdRegions::eHybridDGLamToQ1:
1844  case StdRegions::eHybridDGLamToQ2:
1845  case StdRegions::eHybridDGHelmBndLam:
1846  case StdRegions::eInvLaplacianWithUnityMean:
1847  {
1848  NekDouble one = 1.0;
1849 
1850  DNekMatSharedPtr mat = GenMatrix(mkey);
1851  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1852  }
1853  break;
1854  case StdRegions::eInvHybridDGHelmholtz:
1855  {
1856  NekDouble one = 1.0;
1857 
1858 // StdRegions::StdMatrixKey hkey(StdRegions::eHybridDGHelmholtz,
1859 // DetShapeType(),*this,
1860 // mkey.GetConstant(0),
1861 // mkey.GetConstant(1));
1862  MatrixKey hkey(StdRegions::eHybridDGHelmholtz, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1863  DNekMatSharedPtr mat = GenMatrix(hkey);
1864 
1865  mat->Invert();
1866  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1867  }
1868  break;
1869  case StdRegions::ePreconLinearSpace:
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  DNekMatSharedPtr R=BuildVertexMatrix(A);
1876 
1877  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1878  }
1879  break;
1880  case StdRegions::ePreconLinearSpaceMass:
1881  {
1882  NekDouble one = 1.0;
1883  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1884  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1885  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1886  DNekMatSharedPtr R=BuildVertexMatrix(A);
1887 
1888  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1889  }
1890  break;
1891  case StdRegions::ePreconR:
1892  {
1893  NekDouble one = 1.0;
1894  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this,mkey.GetConstFactors(), mkey.GetVarCoeffs());
1895  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1896  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1897 
1898  DNekScalMatSharedPtr Atmp;
1899  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1900 
1901  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1902  }
1903  break;
1904  case StdRegions::ePreconRT:
1905  {
1906  NekDouble one = 1.0;
1907  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this,mkey.GetConstFactors(), mkey.GetVarCoeffs());
1908  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1909  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1910 
1911  DNekScalMatSharedPtr Atmp;
1912  DNekMatSharedPtr RT=BuildTransformationMatrix(A,mkey.GetMatrixType());
1913 
1914  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,RT);
1915  }
1916  break;
1917  case StdRegions::ePreconRMass:
1918  {
1919  NekDouble one = 1.0;
1920  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1921  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1922  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1923 
1924  DNekScalMatSharedPtr Atmp;
1925  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1926 
1927  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1928  }
1929  break;
1930  case StdRegions::ePreconRTMass:
1931  {
1932  NekDouble one = 1.0;
1933  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1934  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1935  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1936 
1937  DNekScalMatSharedPtr Atmp;
1938  DNekMatSharedPtr RT=BuildTransformationMatrix(A,mkey.GetMatrixType());
1939 
1940  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,RT);
1941  }
1942  break;
1943  default:
1944  {
1945  NekDouble one = 1.0;
1946  DNekMatSharedPtr mat = GenMatrix(mkey);
1947 
1948  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1949  }
1950  break;
1951  }
1952 
1953  return returnval;
1954  }
1955 
1956 
1957  DNekScalBlkMatSharedPtr HexExp::CreateStaticCondMatrix(const MatrixKey &mkey)
1958  {
1959  DNekScalBlkMatSharedPtr returnval;
1960 
1961  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
1962 
1963  // set up block matrix system
1964  unsigned int nbdry = NumBndryCoeffs();
1965  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
1966  unsigned int exp_size[] = {nbdry,nint};
1967  unsigned int nblks = 2;
1968  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks,nblks,exp_size,exp_size); //Really need a constructor which takes Arrays
1969  NekDouble factor = 1.0;
1970 
1971  switch(mkey.GetMatrixType())
1972  {
1973  case StdRegions::eLaplacian:
1974  case StdRegions::eHelmholtz: // special case since Helmholtz not defined in StdRegions
1975 
1976  // use Deformed case for both regular and deformed geometries
1977  factor = 1.0;
1978  goto UseLocRegionsMatrix;
1979  break;
1980  default:
1981  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1982  mkey.GetNVarCoeff())
1983  {
1984  factor = 1.0;
1985  goto UseLocRegionsMatrix;
1986  }
1987  else
1988  {
1989  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1990  factor = mat->Scale();
1991  goto UseStdRegionsMatrix;
1992  }
1993  break;
1994  UseStdRegionsMatrix:
1995  {
1996  NekDouble invfactor = 1.0/factor;
1997  NekDouble one = 1.0;
1998  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
1999  DNekScalMatSharedPtr Atmp;
2000  DNekMatSharedPtr Asubmat;
2001 
2002  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
2003  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
2004  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
2005  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
2006  }
2007  break;
2008  UseLocRegionsMatrix:
2009  {
2010  int i,j;
2011  NekDouble invfactor = 1.0/factor;
2012  NekDouble one = 1.0;
2013  DNekScalMat &mat = *GetLocMatrix(mkey);
2014  DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
2015  DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
2016  DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
2017  DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
2018 
2019  Array<OneD,unsigned int> bmap(nbdry);
2020  Array<OneD,unsigned int> imap(nint);
2021  GetBoundaryMap(bmap);
2022  GetInteriorMap(imap);
2023 
2024  for(i = 0; i < nbdry; ++i)
2025  {
2026  for(j = 0; j < nbdry; ++j)
2027  {
2028  (*A)(i,j) = mat(bmap[i],bmap[j]);
2029  }
2030 
2031  for(j = 0; j < nint; ++j)
2032  {
2033  (*B)(i,j) = mat(bmap[i],imap[j]);
2034  }
2035  }
2036 
2037  for(i = 0; i < nint; ++i)
2038  {
2039  for(j = 0; j < nbdry; ++j)
2040  {
2041  (*C)(i,j) = mat(imap[i],bmap[j]);
2042  }
2043 
2044  for(j = 0; j < nint; ++j)
2045  {
2046  (*D)(i,j) = mat(imap[i],imap[j]);
2047  }
2048  }
2049 
2050  // Calculate static condensed system
2051  if(nint)
2052  {
2053  D->Invert();
2054  (*B) = (*B)*(*D);
2055  (*A) = (*A) - (*B)*(*C);
2056  }
2057 
2058  DNekScalMatSharedPtr Atmp;
2059 
2060  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
2061  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
2062  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
2063  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
2064 
2065  }
2066  }
2067  return returnval;
2068  }
2069 
2070 
2071  DNekScalMatSharedPtr HexExp::v_GetLocMatrix(const MatrixKey &mkey)
2072  {
2073  return m_matrixManager[mkey];
2074  }
2075 
2076 
2077  DNekScalBlkMatSharedPtr HexExp::v_GetLocStaticCondMatrix(
2078  const MatrixKey &mkey)
2079  {
2080  return m_staticCondMatrixManager[mkey];
2081  }
2082 
2083  void HexExp::v_DropLocStaticCondMatrix(const MatrixKey &mkey)
2084  {
2085  m_staticCondMatrixManager.DeleteObject(mkey);
2086  }
2087 
2088  void HexExp::v_LaplacianMatrixOp_MatFree_Kernel(
2089  const Array<OneD, const NekDouble> &inarray,
2090  Array<OneD, NekDouble> &outarray,
2091  Array<OneD, NekDouble> &wsp)
2092  {
2093  // This implementation is only valid when there are no
2094  // coefficients associated to the Laplacian operator
2095  if (m_metrics.count(MetricLaplacian00) == 0)
2096  {
2097  ComputeLaplacianMetric();
2098  }
2099 
2100  int nquad0 = m_base[0]->GetNumPoints();
2101  int nquad1 = m_base[1]->GetNumPoints();
2102  int nquad2 = m_base[2]->GetNumPoints();
2103  int nqtot = nquad0*nquad1*nquad2;
2104 
2105  ASSERTL1(wsp.num_elements() >= 6*nqtot,
2106  "Insufficient workspace size.");
2107 
2108  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
2109  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
2110  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata();
2111  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
2112  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
2113  const Array<OneD, const NekDouble>& dbase2 = m_base[2]->GetDbdata();
2114  const Array<OneD, const NekDouble>& metric00 = m_metrics[MetricLaplacian00];
2115  const Array<OneD, const NekDouble>& metric01 = m_metrics[MetricLaplacian01];
2116  const Array<OneD, const NekDouble>& metric02 = m_metrics[MetricLaplacian02];
2117  const Array<OneD, const NekDouble>& metric11 = m_metrics[MetricLaplacian11];
2118  const Array<OneD, const NekDouble>& metric12 = m_metrics[MetricLaplacian12];
2119  const Array<OneD, const NekDouble>& metric22 = m_metrics[MetricLaplacian22];
2120 
2121  // Allocate temporary storage
2122  Array<OneD,NekDouble> wsp0(wsp);
2123  Array<OneD,NekDouble> wsp1(wsp+1*nqtot);
2124  Array<OneD,NekDouble> wsp2(wsp+2*nqtot);
2125  Array<OneD,NekDouble> wsp3(wsp+3*nqtot);
2126  Array<OneD,NekDouble> wsp4(wsp+4*nqtot);
2127  Array<OneD,NekDouble> wsp5(wsp+5*nqtot);
2128 
2129  StdExpansion3D::PhysTensorDeriv(inarray,wsp0,wsp1,wsp2);
2130 
2131  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
2132  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
2133  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
2134  // especially for this purpose
2135  Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp0[0],1,&metric01[0],1,&wsp1[0],1,&wsp3[0],1);
2136  Vmath::Vvtvp (nqtot,&metric02[0],1,&wsp2[0],1,&wsp3[0],1,&wsp3[0],1);
2137  Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp0[0],1,&metric11[0],1,&wsp1[0],1,&wsp4[0],1);
2138  Vmath::Vvtvp (nqtot,&metric12[0],1,&wsp2[0],1,&wsp4[0],1,&wsp4[0],1);
2139  Vmath::Vvtvvtp(nqtot,&metric02[0],1,&wsp0[0],1,&metric12[0],1,&wsp1[0],1,&wsp5[0],1);
2140  Vmath::Vvtvp (nqtot,&metric22[0],1,&wsp2[0],1,&wsp5[0],1,&wsp5[0],1);
2141 
2142  // outarray = m = (D_xi1 * B)^T * k
2143  // wsp1 = n = (D_xi2 * B)^T * l
2144  IProductWRTBase_SumFacKernel(dbase0,base1,base2,wsp3,outarray,wsp0,false,true,true);
2145  IProductWRTBase_SumFacKernel(base0,dbase1,base2,wsp4,wsp2, wsp0,true,false,true);
2146  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
2147  IProductWRTBase_SumFacKernel(base0,base1,dbase2,wsp5,wsp2, wsp0,true,true,false);
2148  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
2149  }
2150 
2151 
2152  void HexExp::v_ComputeLaplacianMetric()
2153  {
2154  if (m_metrics.count(MetricQuadrature) == 0)
2155  {
2156  ComputeQuadratureMetric();
2157  }
2158 
2159  const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
2160  const unsigned int nqtot = GetTotPoints();
2161  const unsigned int dim = 3;
2162  const MetricType m[3][3] = { {MetricLaplacian00, MetricLaplacian01, MetricLaplacian02},
2163  {MetricLaplacian01, MetricLaplacian11, MetricLaplacian12},
2164  {MetricLaplacian02, MetricLaplacian12, MetricLaplacian22}
2165  };
2166 
2167  for (unsigned int i = 0; i < dim; ++i)
2168  {
2169  for (unsigned int j = i; j < dim; ++j)
2170  {
2171  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
2172  const Array<TwoD, const NekDouble> &gmat =
2173  m_metricinfo->GetGmat(GetPointsKeys());
2174  if (type == SpatialDomains::eDeformed)
2175  {
2176  Vmath::Vcopy(nqtot, &gmat[i*dim+j][0], 1,
2177  &m_metrics[m[i][j]][0], 1);
2178  }
2179  else
2180  {
2181  Vmath::Fill(nqtot, gmat[i*dim+j][0],
2182  &m_metrics[m[i][j]][0], 1);
2183  }
2184  MultiplyByQuadratureMetric(m_metrics[m[i][j]],
2185  m_metrics[m[i][j]]);
2186 
2187  }
2188  }
2189  }
2190 
2191  }//end of namespace
2192 }//end of namespace
Generated on Wed Mar 4 2015 16:33:15 for Nektar++ by   doxygen 1.8.1.2