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