Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HexExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File HexExp.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Methods for Hex expansion in local regoins
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 
37 #include <LocalRegions/HexExp.h>
38 #include <LibUtilities/Foundations/Interp.h>
39 #include <LibUtilities/Foundations/InterpCoeff.h>
40 #include <SpatialDomains/HexGeom.h>
41 
42 namespace Nektar
43 {
44  namespace LocalRegions
45  {
46  /**
47  * @class HexExp
48  * Defines a hexahedral local expansion.
49  */
50 
51  /**
52  * \brief Constructor using BasisKey class for quadrature points and
53  * order definition
54  *
55  * @param Ba Basis key for first coordinate.
56  * @param Bb Basis key for second coordinate.
57  * @param Bc Basis key for third coordinate.
58  */
59  HexExp::HexExp(const LibUtilities::BasisKey &Ba,
60  const LibUtilities::BasisKey &Bb,
61  const LibUtilities::BasisKey &Bc,
62  const SpatialDomains::HexGeomSharedPtr &geom):
63  StdExpansion (Ba.GetNumModes()*Bb.GetNumModes()*Bc.GetNumModes(),3,Ba,Bb,Bc),
64  StdExpansion3D(Ba.GetNumModes()*Bb.GetNumModes()*Bc.GetNumModes(),Ba,Bb,Bc),
65  StdRegions::StdHexExp(Ba,Bb,Bc),
66  Expansion (geom),
67  Expansion3D (geom),
68  m_matrixManager(
69  boost::bind(&HexExp::CreateMatrix, this, _1),
70  std::string("HexExpMatrix")),
71  m_staticCondMatrixManager(
72  boost::bind(&HexExp::CreateStaticCondMatrix, this, _1),
73  std::string("HexExpStaticCondMatrix"))
74  {
75  }
76 
77 
78  /**
79  * \brief Copy Constructor
80  *
81  * @param T HexExp to copy.
82  */
83  HexExp::HexExp(const HexExp &T):
84  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  void HexExp::v_GetFacePhysMap(const int face,
684  Array<OneD, int> &outarray)
685  {
686  int nquad0 = m_base[0]->GetNumPoints();
687  int nquad1 = m_base[1]->GetNumPoints();
688  int nquad2 = m_base[2]->GetNumPoints();
689 
690  int nq0 = 0;
691  int nq1 = 0;
692 
693  switch(face)
694  {
695  case 0:
696  nq0 = nquad0;
697  nq1 = nquad1;
698 
699  //Directions A and B positive
700  if(outarray.num_elements()!=nq0*nq1)
701  {
702  outarray = Array<OneD, int>(nq0*nq1);
703  }
704 
705  for (int i = 0; i < nquad0*nquad1; ++i)
706  {
707  outarray[i] = i;
708  }
709 
710  break;
711  case 1:
712  nq0 = nquad0;
713  nq1 = nquad2;
714  //Direction A and B positive
715  if(outarray.num_elements()!=nq0*nq1)
716  {
717  outarray = Array<OneD, int>(nq0*nq1);
718  }
719 
720  //Direction A and B positive
721  for (int k = 0; k < nquad2; k++)
722  {
723  for(int i = 0; i < nquad0; ++i)
724  {
725  outarray[k*nquad0 + i] = nquad0*nquad1*k + i;
726  }
727  }
728  break;
729  case 2:
730  nq0 = nquad1;
731  nq1 = nquad2;
732 
733  //Direction A and B positive
734  if(outarray.num_elements()!=nq0*nq1)
735  {
736  outarray = Array<OneD, int>(nq0*nq1);
737  }
738 
739  for (int i = 0; i < nquad1*nquad2; i++)
740  {
741  outarray[i] = nquad0-1 + i*nquad0;
742  }
743  break;
744  case 3:
745  nq0 = nquad0;
746  nq1 = nquad2;
747 
748  //Direction A and B positive
749  if(outarray.num_elements()!=nq0*nq1)
750  {
751  outarray = Array<OneD, int>(nq0*nq1);
752  }
753 
754  for (int k = 0; k < nquad2; k++)
755  {
756  for (int i = 0; i < nquad0; i++)
757  {
758  outarray[k*nquad0 + i] = (nquad0*(nquad1-1))+(k*nquad0*nquad1) + i;
759  }
760  }
761  break;
762  case 4:
763  nq0 = nquad1;
764  nq1 = nquad2;
765 
766  //Direction A and B positive
767  if(outarray.num_elements()!=nq0*nq1)
768  {
769  outarray = Array<OneD, int>(nq0*nq1);
770  }
771 
772  for (int i = 0; i < nquad1*nquad2; i++)
773  {
774  outarray[i] = i*nquad0;
775  }
776  break;
777  case 5:
778  nq0 = nquad0;
779  nq1 = nquad1;
780  //Directions A and B positive
781  if(outarray.num_elements()!=nq0*nq1)
782  {
783  outarray = Array<OneD, int>(nq0*nq1);
784  }
785 
786  for (int i = 0; i < nquad0*nquad1; i++)
787  {
788  outarray[i] = nquad0*nquad1*(nquad2-1) + i;
789  }
790 
791  break;
792  default:
793  ASSERTL0(false,"face value (> 5) is out of range");
794  break;
795  }
796 
797  }
798 
799  void HexExp::v_ComputeFaceNormal(const int face)
800  {
801  int i;
802  const SpatialDomains::GeomFactorsSharedPtr & geomFactors =
803  GetGeom()->GetMetricInfo();
804  SpatialDomains::GeomType type = geomFactors->GetGtype();
805  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
806  const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys);
807  const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys);
808 
809  LibUtilities::BasisKey tobasis0 = DetFaceBasisKey(face,0);
810  LibUtilities::BasisKey tobasis1 = DetFaceBasisKey(face,1);
811 
812  // Number of quadrature points in face expansion.
813  int nq_face = tobasis0.GetNumPoints()*tobasis1.GetNumPoints();
814 
815  int vCoordDim = GetCoordim();
816 
817  m_faceNormals[face] = Array<OneD, Array<OneD, NekDouble> >(vCoordDim);
818  Array<OneD, Array<OneD, NekDouble> > &normal = m_faceNormals[face];
819  for (i = 0; i < vCoordDim; ++i)
820  {
821  normal[i] = Array<OneD, NekDouble>(nq_face);
822  }
823  // Regular geometry case
824  if((type == SpatialDomains::eRegular)||(type == SpatialDomains::eMovingRegular))
825  {
826  NekDouble fac;
827  // Set up normals
828  switch(face)
829  {
830  case 0:
831  for(i = 0; i < vCoordDim; ++i)
832  {
833  normal[i][0] = -df[3*i+2][0];
834  }
835  break;
836  case 1:
837  for(i = 0; i < vCoordDim; ++i)
838  {
839  normal[i][0] = -df[3*i+1][0];
840  }
841  break;
842  case 2:
843  for(i = 0; i < vCoordDim; ++i)
844  {
845  normal[i][0] = df[3*i][0];
846  }
847  break;
848  case 3:
849  for(i = 0; i < vCoordDim; ++i)
850  {
851  normal[i][0] = df[3*i+1][0];
852  }
853  break;
854  case 4:
855  for(i = 0; i < vCoordDim; ++i)
856  {
857  normal[i][0] = -df[3*i][0];
858  }
859  break;
860  case 5:
861  for(i = 0; i < vCoordDim; ++i)
862  {
863  normal[i][0] = df[3*i+2][0];
864  }
865  break;
866  default:
867  ASSERTL0(false,"face is out of range (edge < 5)");
868  }
869 
870  // normalise
871  fac = 0.0;
872  for(i =0 ; i < vCoordDim; ++i)
873  {
874  fac += normal[i][0]*normal[i][0];
875  }
876  fac = 1.0/sqrt(fac);
877  for (i = 0; i < vCoordDim; ++i)
878  {
879  Vmath::Fill(nq_face,fac*normal[i][0],normal[i],1);
880  }
881 
882  }
883  else // Set up deformed normals
884  {
885  int j, k;
886 
887  int nqe0 = m_base[0]->GetNumPoints();
888  int nqe1 = m_base[1]->GetNumPoints();
889  int nqe2 = m_base[2]->GetNumPoints();
890  int nqe01 = nqe0*nqe1;
891  int nqe02 = nqe0*nqe2;
892  int nqe12 = nqe1*nqe2;
893 
894  int nqe;
895  if (face == 0 || face == 5)
896  {
897  nqe = nqe01;
898  }
899  else if (face == 1 || face == 3)
900  {
901  nqe = nqe02;
902  }
903  else
904  {
905  nqe = nqe12;
906  }
907 
908  LibUtilities::PointsKey points0;
909  LibUtilities::PointsKey points1;
910 
911  Array<OneD, NekDouble> faceJac(nqe);
912  Array<OneD, NekDouble> normals(vCoordDim*nqe,0.0);
913 
914  // Extract Jacobian along face and recover local
915  // derivates (dx/dr) for polynomial interpolation by
916  // multiplying m_gmat by jacobian
917  switch(face)
918  {
919  case 0:
920  for(j = 0; j < nqe; ++j)
921  {
922  normals[j] = -df[2][j]*jac[j];
923  normals[nqe+j] = -df[5][j]*jac[j];
924  normals[2*nqe+j] = -df[8][j]*jac[j];
925  faceJac[j] = jac[j];
926  }
927 
928  points0 = ptsKeys[0];
929  points1 = ptsKeys[1];
930  break;
931  case 1:
932  for (j = 0; j < nqe0; ++j)
933  {
934  for(k = 0; k < nqe2; ++k)
935  {
936  int idx = j + nqe01*k;
937  normals[j+k*nqe0] = -df[1][idx]*jac[idx];
938  normals[nqe+j+k*nqe0] = -df[4][idx]*jac[idx];
939  normals[2*nqe+j+k*nqe0] = -df[7][idx]*jac[idx];
940  faceJac[j+k*nqe0] = jac[idx];
941  }
942  }
943  points0 = ptsKeys[0];
944  points1 = ptsKeys[2];
945  break;
946  case 2:
947  for (j = 0; j < nqe1; ++j)
948  {
949  for(k = 0; k < nqe2; ++k)
950  {
951  int idx = nqe0-1+nqe0*j+nqe01*k;
952  normals[j+k*nqe1] = df[0][idx]*jac[idx];
953  normals[nqe+j+k*nqe1] = df[3][idx]*jac[idx];
954  normals[2*nqe+j+k*nqe1] = df[6][idx]*jac[idx];
955  faceJac[j+k*nqe1] = jac[idx];
956  }
957  }
958  points0 = ptsKeys[1];
959  points1 = ptsKeys[2];
960  break;
961  case 3:
962  for (j = 0; j < nqe0; ++j)
963  {
964  for(k = 0; k < nqe2; ++k)
965  {
966  int idx = nqe0*(nqe1-1)+j+nqe01*k;
967  normals[j+k*nqe0] = df[1][idx]*jac[idx];
968  normals[nqe+j+k*nqe0] = df[4][idx]*jac[idx];
969  normals[2*nqe+j+k*nqe0] = df[7][idx]*jac[idx];
970  faceJac[j+k*nqe0] = jac[idx];
971  }
972  }
973  points0 = ptsKeys[0];
974  points1 = ptsKeys[2];
975  break;
976  case 4:
977  for (j = 0; j < nqe1; ++j)
978  {
979  for(k = 0; k < nqe2; ++k)
980  {
981  int idx = j*nqe0+nqe01*k;
982  normals[j+k*nqe1] = -df[0][idx]*jac[idx];
983  normals[nqe+j+k*nqe1] = -df[3][idx]*jac[idx];
984  normals[2*nqe+j+k*nqe1] = -df[6][idx]*jac[idx];
985  faceJac[j+k*nqe1] = jac[idx];
986  }
987  }
988  points0 = ptsKeys[1];
989  points1 = ptsKeys[2];
990  break;
991  case 5:
992  for (j = 0; j < nqe01; ++j)
993  {
994  int idx = j+nqe01*(nqe2-1);
995  normals[j] = df[2][idx]*jac[idx];
996  normals[nqe+j] = df[5][idx]*jac[idx];
997  normals[2*nqe+j] = df[8][idx]*jac[idx];
998  faceJac[j] = jac[idx];
999  }
1000  points0 = ptsKeys[0];
1001  points1 = ptsKeys[1];
1002  break;
1003  default:
1004  ASSERTL0(false,"face is out of range (face < 5)");
1005  }
1006 
1007  Array<OneD, NekDouble> work (nq_face, 0.0);
1008  // Interpolate Jacobian and invert
1009  LibUtilities::Interp2D(points0, points1, faceJac,
1010  tobasis0.GetPointsKey(),
1011  tobasis1.GetPointsKey(),
1012  work);
1013 
1014  Vmath::Sdiv(nq_face,1.0,&work[0],1,&work[0],1);
1015 
1016  // interpolate
1017  for(i = 0; i < GetCoordim(); ++i)
1018  {
1019  LibUtilities::Interp2D(points0, points1,
1020  &normals[i*nqe],
1021  tobasis0.GetPointsKey(),
1022  tobasis1.GetPointsKey(),
1023  &normal[i][0]);
1024  Vmath::Vmul(nq_face,work,1,normal[i],1,normal[i],1);
1025  }
1026 
1027  //normalise normal vectors
1028  Vmath::Zero(nq_face,work,1);
1029  for(i = 0; i < GetCoordim(); ++i)
1030  {
1031  Vmath::Vvtvp(nq_face,normal[i],1, normal[i],1,work,1,work,1);
1032  }
1033 
1034  Vmath::Vsqrt(nq_face,work,1,work,1);
1035  Vmath::Sdiv(nq_face,1.0,work,1,work,1);
1036 
1037  for(i = 0; i < GetCoordim(); ++i)
1038  {
1039  Vmath::Vmul(nq_face,normal[i],1,work,1,normal[i],1);
1040  }
1041  }
1042  }
1043 
1044  //-----------------------------
1045  // Operator creation functions
1046  //-----------------------------
1047  void HexExp::v_MassMatrixOp(
1048  const Array<OneD, const NekDouble> &inarray,
1049  Array<OneD,NekDouble> &outarray,
1050  const StdRegions::StdMatrixKey &mkey)
1051  {
1052  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
1053  }
1054 
1055  void HexExp::v_LaplacianMatrixOp(
1056  const Array<OneD, const NekDouble> &inarray,
1057  Array<OneD,NekDouble> &outarray,
1058  const StdRegions::StdMatrixKey &mkey)
1059  {
1060  HexExp::v_LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1061  }
1062 
1063  void HexExp::v_LaplacianMatrixOp(
1064  const int k1,
1065  const int k2,
1066  const Array<OneD, const NekDouble> &inarray,
1067  Array<OneD,NekDouble> &outarray,
1068  const StdRegions::StdMatrixKey &mkey)
1069  {
1070  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,
1071  mkey);
1072  }
1073 
1074  void HexExp::v_WeakDerivMatrixOp(
1075  const int i,
1076  const Array<OneD, const NekDouble> &inarray,
1077  Array<OneD,NekDouble> &outarray,
1078  const StdRegions::StdMatrixKey &mkey)
1079  {
1080  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
1081  }
1082 
1083  void HexExp::v_WeakDirectionalDerivMatrixOp(
1084  const Array<OneD, const NekDouble> &inarray,
1085  Array<OneD,NekDouble> &outarray,
1086  const StdRegions::StdMatrixKey &mkey)
1087  {
1088  StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray,
1089  outarray,mkey);
1090  }
1091 
1092  void HexExp::v_MassLevelCurvatureMatrixOp(
1093  const Array<OneD, const NekDouble> &inarray,
1094  Array<OneD,NekDouble> &outarray,
1095  const StdRegions::StdMatrixKey &mkey)
1096  {
1097  StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray,
1098  outarray,mkey);
1099  }
1100 
1101  void HexExp::v_HelmholtzMatrixOp(
1102  const Array<OneD, const NekDouble> &inarray,
1103  Array<OneD,NekDouble> &outarray,
1104  const StdRegions::StdMatrixKey &mkey)
1105  {
1106  HexExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1107  }
1108 
1109 
1110  void HexExp::v_GeneralMatrixOp_MatOp(
1111  const Array<OneD, const NekDouble> &inarray,
1112  Array<OneD,NekDouble> &outarray,
1113  const StdRegions::StdMatrixKey &mkey)
1114  {
1115  //int nConsts = mkey.GetNconstants();
1116  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1117 
1118 // switch(nConsts)
1119 // {
1120 // case 0:
1121 // {
1122 // mat = GetLocMatrix(mkey.GetMatrixType());
1123 // }
1124 // break;
1125 // case 1:
1126 // {
1127 // mat = GetLocMatrix(mkey.GetMatrixType(),mkey.GetConstant(0));
1128 // }
1129 // break;
1130 // case 2:
1131 // {
1132 // mat = GetLocMatrix(mkey.GetMatrixType(),mkey.GetConstant(0),mkey.GetConstant(1));
1133 // }
1134 // break;
1135 //
1136 // default:
1137 // {
1138 // NEKERROR(ErrorUtil::efatal, "Unknown number of constants");
1139 // }
1140 // break;
1141 // }
1142 
1143  if(inarray.get() == outarray.get())
1144  {
1145  Array<OneD,NekDouble> tmp(m_ncoeffs);
1146  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1147 
1148  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1149  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1150  }
1151  else
1152  {
1153  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1154  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1155  }
1156  }
1157 
1158  /**
1159  * This function is used to compute exactly the advective numerical flux
1160  * on the interface of two elements with different expansions, hence an
1161  * appropriate number of Gauss points has to be used. The number of
1162  * Gauss points has to be equal to the number used by the highest
1163  * polynomial degree of the two adjacent elements
1164  *
1165  * @param numMin Is the reduced polynomial order
1166  * @param inarray Input array of coefficients
1167  * @param dumpVar Output array of reduced coefficients.
1168  */
1169  void HexExp::v_ReduceOrderCoeffs(
1170  int numMin,
1171  const Array<OneD, const NekDouble> &inarray,
1172  Array<OneD, NekDouble> &outarray)
1173  {
1174  int n_coeffs = inarray.num_elements();
1175  int nmodes0 = m_base[0]->GetNumModes();
1176  int nmodes1 = m_base[1]->GetNumModes();
1177  int nmodes2 = m_base[2]->GetNumModes();
1178  int numMax = nmodes0;
1179 
1180  Array<OneD, NekDouble> coeff (n_coeffs);
1181  Array<OneD, NekDouble> coeff_tmp1(nmodes0*nmodes1, 0.0);
1182  Array<OneD, NekDouble> coeff_tmp2(n_coeffs, 0.0);
1183  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
1184 
1185  Vmath::Vcopy(n_coeffs,inarray,1,coeff_tmp2,1);
1186 
1187  const LibUtilities::PointsKey Pkey0(
1188  nmodes0, LibUtilities::eGaussLobattoLegendre);
1189  const LibUtilities::PointsKey Pkey1(
1190  nmodes1, LibUtilities::eGaussLobattoLegendre);
1191  const LibUtilities::PointsKey Pkey2(
1192  nmodes2, LibUtilities::eGaussLobattoLegendre);
1193 
1194  LibUtilities::BasisKey b0(
1195  m_base[0]->GetBasisType(), nmodes0, Pkey0);
1196  LibUtilities::BasisKey b1(
1197  m_base[1]->GetBasisType(), nmodes1, Pkey1);
1198  LibUtilities::BasisKey b2(
1199  m_base[2]->GetBasisType(), nmodes2, Pkey2);
1200  LibUtilities::BasisKey bortho0(
1201  LibUtilities::eOrtho_A, nmodes0, Pkey0);
1202  LibUtilities::BasisKey bortho1(
1203  LibUtilities::eOrtho_A, nmodes1, Pkey1);
1204  LibUtilities::BasisKey bortho2(
1205  LibUtilities::eOrtho_A, nmodes2, Pkey2);
1206 
1207  LibUtilities::InterpCoeff3D(
1208  b0, b1, b2, coeff_tmp2,
1209  bortho0, bortho1, bortho2, coeff);
1210 
1211  Vmath::Zero(n_coeffs, coeff_tmp2, 1);
1212 
1213  int cnt = 0, cnt2 = 0;
1214 
1215  for (int u = 0; u < numMin+1; ++u)
1216  {
1217  for (int i = 0; i < numMin; ++i)
1218  {
1219  Vmath::Vcopy(numMin,
1220  tmp = coeff+cnt+cnt2,1,
1221  tmp2 = coeff_tmp1+cnt,1);
1222 
1223  cnt = i*numMax;
1224  }
1225 
1226  Vmath::Vcopy(nmodes0*nmodes1,
1227  tmp3 = coeff_tmp1,1,
1228  tmp4 = coeff_tmp2+cnt2,1);
1229 
1230  cnt2 = u*nmodes0*nmodes1;
1231  }
1232 
1233  LibUtilities::InterpCoeff3D(
1234  bortho0, bortho1, bortho2, coeff_tmp2,
1235  b0, b1, b2, outarray);
1236  }
1237 
1238  void HexExp::v_SVVLaplacianFilter(
1239  Array<OneD, NekDouble> &array,
1240  const StdRegions::StdMatrixKey &mkey)
1241  {
1242  int nq = GetTotPoints();
1243 
1244  // Calculate sqrt of the Jacobian
1245  Array<OneD, const NekDouble> jac =
1246  m_metricinfo->GetJac(GetPointsKeys());
1247  Array<OneD, NekDouble> sqrt_jac(nq);
1248  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1249  {
1250  Vmath::Vsqrt(nq,jac,1,sqrt_jac,1);
1251  }
1252  else
1253  {
1254  Vmath::Fill(nq,sqrt(jac[0]),sqrt_jac,1);
1255  }
1256 
1257  // Multiply array by sqrt(Jac)
1258  Vmath::Vmul(nq,sqrt_jac,1,array,1,array,1);
1259 
1260  // Apply std region filter
1261  StdHexExp::v_SVVLaplacianFilter( array, mkey);
1262 
1263  // Divide by sqrt(Jac)
1264  Vmath::Vdiv(nq,array,1,sqrt_jac,1,array,1);
1265  }
1266 
1267  //-----------------------------
1268  // Matrix creation functions
1269  //-----------------------------
1270  DNekMatSharedPtr HexExp::v_GenMatrix(
1271  const StdRegions::StdMatrixKey &mkey)
1272  {
1273  DNekMatSharedPtr returnval;
1274 
1275  switch(mkey.GetMatrixType())
1276  {
1277  case StdRegions::eHybridDGHelmholtz:
1278  case StdRegions::eHybridDGLamToU:
1279  case StdRegions::eHybridDGLamToQ0:
1280  case StdRegions::eHybridDGLamToQ1:
1281  case StdRegions::eHybridDGLamToQ2:
1282  case StdRegions::eHybridDGHelmBndLam:
1283  case StdRegions::eInvLaplacianWithUnityMean:
1284  returnval = Expansion3D::v_GenMatrix(mkey);
1285  break;
1286  default:
1287  returnval = StdHexExp::v_GenMatrix(mkey);
1288  }
1289 
1290  return returnval;
1291  }
1292 
1293 
1294  DNekMatSharedPtr HexExp::v_CreateStdMatrix(
1295  const StdRegions::StdMatrixKey &mkey)
1296  {
1297  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1298  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1299  LibUtilities::BasisKey bkey2 = m_base[2]->GetBasisKey();
1300 
1301  StdRegions::StdHexExpSharedPtr tmp = MemoryManager<StdHexExp>
1302  ::AllocateSharedPtr(bkey0,bkey1,bkey2);
1303 
1304  return tmp->GetStdMatrix(mkey);
1305  }
1306 
1307 
1308  DNekScalMatSharedPtr HexExp::CreateMatrix(const MatrixKey &mkey)
1309  {
1310  DNekScalMatSharedPtr returnval;
1311  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1312 
1313  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,
1314  "Geometric information is not set up");
1315 
1316  switch(mkey.GetMatrixType())
1317  {
1318  case StdRegions::eMass:
1319  {
1320  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1321  mkey.GetNVarCoeff())
1322  {
1323  NekDouble one = 1.0;
1324  DNekMatSharedPtr mat = GenMatrix(mkey);
1325  returnval = MemoryManager<DNekScalMat>
1326  ::AllocateSharedPtr(one,mat);
1327  }
1328  else
1329  {
1330  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1331  DNekMatSharedPtr mat
1332  = GetStdMatrix(mkey);
1333  returnval = MemoryManager<DNekScalMat>
1334  ::AllocateSharedPtr(jac,mat);
1335  }
1336  }
1337  break;
1338  case StdRegions::eInvMass:
1339  {
1340  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1341  {
1342  NekDouble one = 1.0;
1343  StdRegions::StdMatrixKey masskey(StdRegions::eMass,
1344  DetShapeType(), *this);
1345  DNekMatSharedPtr mat = GenMatrix(masskey);
1346  mat->Invert();
1347 
1348  returnval = MemoryManager<DNekScalMat>
1349  ::AllocateSharedPtr(one,mat);
1350  }
1351  else
1352  {
1353  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1354  DNekMatSharedPtr mat
1355  = GetStdMatrix(mkey);
1356  returnval = MemoryManager<DNekScalMat>
1357  ::AllocateSharedPtr(fac,mat);
1358  }
1359  }
1360  break;
1361  case StdRegions::eWeakDeriv0:
1362  case StdRegions::eWeakDeriv1:
1363  case StdRegions::eWeakDeriv2:
1364  {
1365  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1366  mkey.GetNVarCoeff())
1367  {
1368  NekDouble one = 1.0;
1369  DNekMatSharedPtr mat = GenMatrix(mkey);
1370 
1371  returnval = MemoryManager<DNekScalMat>
1372  ::AllocateSharedPtr(one,mat);
1373  }
1374  else
1375  {
1376  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1377  Array<TwoD, const NekDouble> df
1378  = m_metricinfo->GetDerivFactors(ptsKeys);
1379  int dir = 0;
1380 
1381  switch(mkey.GetMatrixType())
1382  {
1383  case StdRegions::eWeakDeriv0:
1384  dir = 0;
1385  break;
1386  case StdRegions::eWeakDeriv1:
1387  dir = 1;
1388  break;
1389  case StdRegions::eWeakDeriv2:
1390  dir = 2;
1391  break;
1392  default:
1393  break;
1394  }
1395 
1396  MatrixKey deriv0key(StdRegions::eWeakDeriv0,
1397  mkey.GetShapeType(), *this);
1398  MatrixKey deriv1key(StdRegions::eWeakDeriv1,
1399  mkey.GetShapeType(), *this);
1400  MatrixKey deriv2key(StdRegions::eWeakDeriv2,
1401  mkey.GetShapeType(), *this);
1402 
1403  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
1404  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
1405  DNekMat &deriv2 = *GetStdMatrix(deriv2key);
1406 
1407  int rows = deriv0.GetRows();
1408  int cols = deriv1.GetColumns();
1409 
1410  DNekMatSharedPtr WeakDeriv = MemoryManager<DNekMat>
1411  ::AllocateSharedPtr(rows,cols);
1412 
1413  (*WeakDeriv) = df[3*dir ][0]*deriv0
1414  + df[3*dir+1][0]*deriv1
1415  + df[3*dir+2][0]*deriv2;
1416 
1417  returnval = MemoryManager<DNekScalMat>
1418  ::AllocateSharedPtr(jac,WeakDeriv);
1419  }
1420  }
1421  break;
1422  case StdRegions::eLaplacian:
1423  {
1424  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1425  mkey.GetNVarCoeff()||
1426  mkey.ConstFactorExists(
1427  StdRegions::eFactorSVVCutoffRatio))
1428  {
1429  NekDouble one = 1.0;
1430  DNekMatSharedPtr mat = GenMatrix(mkey);
1431 
1432  returnval = MemoryManager<DNekScalMat>
1433  ::AllocateSharedPtr(one,mat);
1434  }
1435  else
1436  {
1437  MatrixKey lap00key(StdRegions::eLaplacian00,
1438  mkey.GetShapeType(), *this);
1439  MatrixKey lap01key(StdRegions::eLaplacian01,
1440  mkey.GetShapeType(), *this);
1441  MatrixKey lap02key(StdRegions::eLaplacian02,
1442  mkey.GetShapeType(), *this);
1443  MatrixKey lap11key(StdRegions::eLaplacian11,
1444  mkey.GetShapeType(), *this);
1445  MatrixKey lap12key(StdRegions::eLaplacian12,
1446  mkey.GetShapeType(), *this);
1447  MatrixKey lap22key(StdRegions::eLaplacian22,
1448  mkey.GetShapeType(), *this);
1449 
1450  DNekMat &lap00 = *GetStdMatrix(lap00key);
1451  DNekMat &lap01 = *GetStdMatrix(lap01key);
1452  DNekMat &lap02 = *GetStdMatrix(lap02key);
1453  DNekMat &lap11 = *GetStdMatrix(lap11key);
1454  DNekMat &lap12 = *GetStdMatrix(lap12key);
1455  DNekMat &lap22 = *GetStdMatrix(lap22key);
1456 
1457  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1458  Array<TwoD, const NekDouble> gmat
1459  = m_metricinfo->GetGmat(ptsKeys);
1460 
1461  int rows = lap00.GetRows();
1462  int cols = lap00.GetColumns();
1463 
1464  DNekMatSharedPtr lap = MemoryManager<DNekMat>
1465  ::AllocateSharedPtr(rows,cols);
1466 
1467  (*lap) = gmat[0][0]*lap00
1468  + gmat[4][0]*lap11
1469  + gmat[8][0]*lap22
1470  + gmat[3][0]*(lap01 + Transpose(lap01))
1471  + gmat[6][0]*(lap02 + Transpose(lap02))
1472  + gmat[7][0]*(lap12 + Transpose(lap12));
1473 
1474  returnval = MemoryManager<DNekScalMat>
1475  ::AllocateSharedPtr(jac,lap);
1476  }
1477  }
1478  break;
1479  case StdRegions::eHelmholtz:
1480  {
1481  NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
1482  MatrixKey masskey(StdRegions::eMass,
1483  mkey.GetShapeType(), *this);
1484  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1485  MatrixKey lapkey(StdRegions::eLaplacian,
1486  mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1487  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1488 
1489  int rows = LapMat.GetRows();
1490  int cols = LapMat.GetColumns();
1491 
1492  DNekMatSharedPtr helm = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1493 
1494  NekDouble one = 1.0;
1495  (*helm) = LapMat + lambda*MassMat;
1496 
1497  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
1498  }
1499  break;
1500  case StdRegions::eIProductWRTBase:
1501  {
1502  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1503  {
1504  NekDouble one = 1.0;
1505  DNekMatSharedPtr mat = GenMatrix(mkey);
1506  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1507  }
1508  else
1509  {
1510  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1511  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1512  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1513  }
1514  }
1515  break;
1516  case StdRegions::eHybridDGHelmholtz:
1517  case StdRegions::eHybridDGLamToU:
1518  case StdRegions::eHybridDGLamToQ0:
1519  case StdRegions::eHybridDGLamToQ1:
1520  case StdRegions::eHybridDGLamToQ2:
1521  case StdRegions::eHybridDGHelmBndLam:
1522  case StdRegions::eInvLaplacianWithUnityMean:
1523  {
1524  NekDouble one = 1.0;
1525 
1526  DNekMatSharedPtr mat = GenMatrix(mkey);
1527  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1528  }
1529  break;
1530  case StdRegions::eInvHybridDGHelmholtz:
1531  {
1532  NekDouble one = 1.0;
1533 
1534 // StdRegions::StdMatrixKey hkey(StdRegions::eHybridDGHelmholtz,
1535 // DetShapeType(),*this,
1536 // mkey.GetConstant(0),
1537 // mkey.GetConstant(1));
1538  MatrixKey hkey(StdRegions::eHybridDGHelmholtz, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1539  DNekMatSharedPtr mat = GenMatrix(hkey);
1540 
1541  mat->Invert();
1542  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1543  }
1544  break;
1545  case StdRegions::ePreconLinearSpace:
1546  {
1547  NekDouble one = 1.0;
1548  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1549  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1550  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1551  DNekMatSharedPtr R=BuildVertexMatrix(A);
1552 
1553  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1554  }
1555  break;
1556  case StdRegions::ePreconLinearSpaceMass:
1557  {
1558  NekDouble one = 1.0;
1559  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1560  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1561  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1562  DNekMatSharedPtr R=BuildVertexMatrix(A);
1563 
1564  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1565  }
1566  break;
1567  case StdRegions::ePreconR:
1568  {
1569  NekDouble one = 1.0;
1570  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this,mkey.GetConstFactors(), mkey.GetVarCoeffs());
1571  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1572  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1573 
1574  DNekScalMatSharedPtr Atmp;
1575  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1576 
1577  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1578  }
1579  break;
1580  case StdRegions::ePreconRT:
1581  {
1582  NekDouble one = 1.0;
1583  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this,mkey.GetConstFactors(), mkey.GetVarCoeffs());
1584  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1585  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1586 
1587  DNekScalMatSharedPtr Atmp;
1588  DNekMatSharedPtr RT=BuildTransformationMatrix(A,mkey.GetMatrixType());
1589 
1590  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,RT);
1591  }
1592  break;
1593  case StdRegions::ePreconRMass:
1594  {
1595  NekDouble one = 1.0;
1596  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1597  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1598  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1599 
1600  DNekScalMatSharedPtr Atmp;
1601  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1602 
1603  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1604  }
1605  break;
1606  case StdRegions::ePreconRTMass:
1607  {
1608  NekDouble one = 1.0;
1609  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1610  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1611  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1612 
1613  DNekScalMatSharedPtr Atmp;
1614  DNekMatSharedPtr RT=BuildTransformationMatrix(A,mkey.GetMatrixType());
1615 
1616  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,RT);
1617  }
1618  break;
1619  default:
1620  {
1621  NekDouble one = 1.0;
1622  DNekMatSharedPtr mat = GenMatrix(mkey);
1623 
1624  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1625  }
1626  break;
1627  }
1628 
1629  return returnval;
1630  }
1631 
1632 
1633  DNekScalBlkMatSharedPtr HexExp::CreateStaticCondMatrix(const MatrixKey &mkey)
1634  {
1635  DNekScalBlkMatSharedPtr returnval;
1636 
1637  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
1638 
1639  // set up block matrix system
1640  unsigned int nbdry = NumBndryCoeffs();
1641  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
1642  unsigned int exp_size[] = {nbdry,nint};
1643  unsigned int nblks = 2;
1644  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks,nblks,exp_size,exp_size); //Really need a constructor which takes Arrays
1645  NekDouble factor = 1.0;
1646 
1647  switch(mkey.GetMatrixType())
1648  {
1649  case StdRegions::eLaplacian:
1650  case StdRegions::eHelmholtz: // special case since Helmholtz not defined in StdRegions
1651 
1652  // use Deformed case for both regular and deformed geometries
1653  factor = 1.0;
1654  goto UseLocRegionsMatrix;
1655  break;
1656  default:
1657  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1658  mkey.GetNVarCoeff())
1659  {
1660  factor = 1.0;
1661  goto UseLocRegionsMatrix;
1662  }
1663  else
1664  {
1665  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1666  factor = mat->Scale();
1667  goto UseStdRegionsMatrix;
1668  }
1669  break;
1670  UseStdRegionsMatrix:
1671  {
1672  NekDouble invfactor = 1.0/factor;
1673  NekDouble one = 1.0;
1674  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
1675  DNekScalMatSharedPtr Atmp;
1676  DNekMatSharedPtr Asubmat;
1677 
1678  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
1679  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
1680  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
1681  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
1682  }
1683  break;
1684  UseLocRegionsMatrix:
1685  {
1686  int i,j;
1687  NekDouble invfactor = 1.0/factor;
1688  NekDouble one = 1.0;
1689  DNekScalMat &mat = *GetLocMatrix(mkey);
1690  DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
1691  DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
1692  DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
1693  DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
1694 
1695  Array<OneD,unsigned int> bmap(nbdry);
1696  Array<OneD,unsigned int> imap(nint);
1697  GetBoundaryMap(bmap);
1698  GetInteriorMap(imap);
1699 
1700  for(i = 0; i < nbdry; ++i)
1701  {
1702  for(j = 0; j < nbdry; ++j)
1703  {
1704  (*A)(i,j) = mat(bmap[i],bmap[j]);
1705  }
1706 
1707  for(j = 0; j < nint; ++j)
1708  {
1709  (*B)(i,j) = mat(bmap[i],imap[j]);
1710  }
1711  }
1712 
1713  for(i = 0; i < nint; ++i)
1714  {
1715  for(j = 0; j < nbdry; ++j)
1716  {
1717  (*C)(i,j) = mat(imap[i],bmap[j]);
1718  }
1719 
1720  for(j = 0; j < nint; ++j)
1721  {
1722  (*D)(i,j) = mat(imap[i],imap[j]);
1723  }
1724  }
1725 
1726  // Calculate static condensed system
1727  if(nint)
1728  {
1729  D->Invert();
1730  (*B) = (*B)*(*D);
1731  (*A) = (*A) - (*B)*(*C);
1732  }
1733 
1734  DNekScalMatSharedPtr Atmp;
1735 
1736  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
1737  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
1738  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
1739  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
1740 
1741  }
1742  }
1743  return returnval;
1744  }
1745 
1746 
1747  DNekScalMatSharedPtr HexExp::v_GetLocMatrix(const MatrixKey &mkey)
1748  {
1749  return m_matrixManager[mkey];
1750  }
1751 
1752 
1753  DNekScalBlkMatSharedPtr HexExp::v_GetLocStaticCondMatrix(
1754  const MatrixKey &mkey)
1755  {
1756  return m_staticCondMatrixManager[mkey];
1757  }
1758 
1759  void HexExp::v_DropLocStaticCondMatrix(const MatrixKey &mkey)
1760  {
1761  m_staticCondMatrixManager.DeleteObject(mkey);
1762  }
1763 
1764  void HexExp::v_LaplacianMatrixOp_MatFree_Kernel(
1765  const Array<OneD, const NekDouble> &inarray,
1766  Array<OneD, NekDouble> &outarray,
1767  Array<OneD, NekDouble> &wsp)
1768  {
1769  // This implementation is only valid when there are no
1770  // coefficients associated to the Laplacian operator
1771  if (m_metrics.count(eMetricLaplacian00) == 0)
1772  {
1773  ComputeLaplacianMetric();
1774  }
1775 
1776  int nquad0 = m_base[0]->GetNumPoints();
1777  int nquad1 = m_base[1]->GetNumPoints();
1778  int nquad2 = m_base[2]->GetNumPoints();
1779  int nqtot = nquad0*nquad1*nquad2;
1780 
1781  ASSERTL1(wsp.num_elements() >= 6*nqtot,
1782  "Insufficient workspace size.");
1783 
1784  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
1785  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
1786  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata();
1787  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
1788  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
1789  const Array<OneD, const NekDouble>& dbase2 = m_base[2]->GetDbdata();
1790  const Array<OneD, const NekDouble>& metric00 = m_metrics[eMetricLaplacian00];
1791  const Array<OneD, const NekDouble>& metric01 = m_metrics[eMetricLaplacian01];
1792  const Array<OneD, const NekDouble>& metric02 = m_metrics[eMetricLaplacian02];
1793  const Array<OneD, const NekDouble>& metric11 = m_metrics[eMetricLaplacian11];
1794  const Array<OneD, const NekDouble>& metric12 = m_metrics[eMetricLaplacian12];
1795  const Array<OneD, const NekDouble>& metric22 = m_metrics[eMetricLaplacian22];
1796 
1797  // Allocate temporary storage
1798  Array<OneD,NekDouble> wsp0(wsp);
1799  Array<OneD,NekDouble> wsp1(wsp+1*nqtot);
1800  Array<OneD,NekDouble> wsp2(wsp+2*nqtot);
1801  Array<OneD,NekDouble> wsp3(wsp+3*nqtot);
1802  Array<OneD,NekDouble> wsp4(wsp+4*nqtot);
1803  Array<OneD,NekDouble> wsp5(wsp+5*nqtot);
1804 
1805  StdExpansion3D::PhysTensorDeriv(inarray,wsp0,wsp1,wsp2);
1806 
1807  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1808  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1809  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1810  // especially for this purpose
1811  Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp0[0],1,&metric01[0],1,&wsp1[0],1,&wsp3[0],1);
1812  Vmath::Vvtvp (nqtot,&metric02[0],1,&wsp2[0],1,&wsp3[0],1,&wsp3[0],1);
1813  Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp0[0],1,&metric11[0],1,&wsp1[0],1,&wsp4[0],1);
1814  Vmath::Vvtvp (nqtot,&metric12[0],1,&wsp2[0],1,&wsp4[0],1,&wsp4[0],1);
1815  Vmath::Vvtvvtp(nqtot,&metric02[0],1,&wsp0[0],1,&metric12[0],1,&wsp1[0],1,&wsp5[0],1);
1816  Vmath::Vvtvp (nqtot,&metric22[0],1,&wsp2[0],1,&wsp5[0],1,&wsp5[0],1);
1817 
1818  // outarray = m = (D_xi1 * B)^T * k
1819  // wsp1 = n = (D_xi2 * B)^T * l
1820  IProductWRTBase_SumFacKernel(dbase0,base1,base2,wsp3,outarray,wsp0,false,true,true);
1821  IProductWRTBase_SumFacKernel(base0,dbase1,base2,wsp4,wsp2, wsp0,true,false,true);
1822  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
1823  IProductWRTBase_SumFacKernel(base0,base1,dbase2,wsp5,wsp2, wsp0,true,true,false);
1824  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
1825  }
1826 
1827 
1828  void HexExp::v_ComputeLaplacianMetric()
1829  {
1830  if (m_metrics.count(eMetricQuadrature) == 0)
1831  {
1832  ComputeQuadratureMetric();
1833  }
1834 
1835  const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
1836  const unsigned int nqtot = GetTotPoints();
1837  const unsigned int dim = 3;
1838  const MetricType m[3][3] = { {eMetricLaplacian00, eMetricLaplacian01, eMetricLaplacian02},
1839  {eMetricLaplacian01, eMetricLaplacian11, eMetricLaplacian12},
1840  {eMetricLaplacian02, eMetricLaplacian12, eMetricLaplacian22}
1841  };
1842 
1843  for (unsigned int i = 0; i < dim; ++i)
1844  {
1845  for (unsigned int j = i; j < dim; ++j)
1846  {
1847  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1848  const Array<TwoD, const NekDouble> &gmat =
1849  m_metricinfo->GetGmat(GetPointsKeys());
1850  if (type == SpatialDomains::eDeformed)
1851  {
1852  Vmath::Vcopy(nqtot, &gmat[i*dim+j][0], 1,
1853  &m_metrics[m[i][j]][0], 1);
1854  }
1855  else
1856  {
1857  Vmath::Fill(nqtot, gmat[i*dim+j][0],
1858  &m_metrics[m[i][j]][0], 1);
1859  }
1860  MultiplyByQuadratureMetric(m_metrics[m[i][j]],
1861  m_metrics[m[i][j]]);
1862 
1863  }
1864  }
1865  }
1866 
1867  }//end of namespace
1868 }//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:161
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:1101
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: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:607
Nektar::LocalRegions::HexExp::v_ComputeFaceNormal
void v_ComputeFaceNormal(const int face)
Definition: HexExp.cpp:799
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:1828
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:1747
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:114
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:1759
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::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:1074
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: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:1753
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:1047
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:673
Nektar::LocalRegions::HexExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1238
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:1110
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:683
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: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:1308
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:1270
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: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:1083
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:1169
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:1764
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: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:381
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:213
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:1633
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:559
Nektar::LocalRegions::HexExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1294
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:1055
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: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:191
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: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:1092
Generated on Wed Mar 23 2016 20:14:20 for Nektar++ by   doxygen 1.8.8