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