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