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.size())
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.size())
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.size())
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.size())
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.size())
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.size())
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.size());
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  Array<OneD, NekDouble> alloc(4*nq + m_ncoeffs + nm0*nq2*(nq1+nm1));
482  Array<OneD, NekDouble> tmp1 (alloc); // Quad metric
483  Array<OneD, NekDouble> tmp2 (alloc + nq); // Dir1 metric
484  Array<OneD, NekDouble> tmp3 (alloc + 2*nq); // Dir2 metric
485  Array<OneD, NekDouble> tmp4 (alloc + 3*nq); // Dir3 metric
486  Array<OneD, NekDouble> tmp5 (alloc + 4*nq); // iprod tmp
487  Array<OneD, NekDouble> wsp (tmp5 + m_ncoeffs); // Wsp
488 
489  MultiplyByQuadratureMetric(inarray, tmp1);
490 
491  Array<OneD, Array<OneD, NekDouble>> tmp2D{3};
492  tmp2D[0] = tmp2;
493  tmp2D[1] = tmp3;
494  tmp2D[2] = tmp4;
495 
496  HexExp::v_AlignVectorToCollapsedDir(dir, tmp1, tmp2D);
497 
498  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
499  m_base[1]->GetBdata(),
500  m_base[2]->GetBdata(),
501  tmp2,outarray,wsp,
502  false,true,true);
503 
504  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
505  m_base[1]->GetDbdata(),
506  m_base[2]->GetBdata(),
507  tmp3,tmp5,wsp,
508  true,false,true);
509  Vmath::Vadd(m_ncoeffs, tmp5, 1, outarray, 1, outarray, 1);
510 
511  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
512  m_base[1]->GetBdata(),
513  m_base[2]->GetDbdata(),
514  tmp4,tmp5,wsp,
515  true,true,false);
516  Vmath::Vadd(m_ncoeffs, tmp5, 1, outarray, 1, outarray, 1);
517  }
518 
519  void HexExp::v_AlignVectorToCollapsedDir(
520  const int dir,
521  const Array<OneD, const NekDouble> &inarray,
522  Array<OneD, Array<OneD, NekDouble> > &outarray)
523  {
524  ASSERTL1((dir==0)||(dir==1)||(dir==2),"Invalid direction.");
525 
526  const int nq0 = m_base[0]->GetNumPoints();
527  const int nq1 = m_base[1]->GetNumPoints();
528  const int nq2 = m_base[2]->GetNumPoints();
529  const int nq = nq0*nq1*nq2;
530 
531  const Array<TwoD, const NekDouble>& df =
532  m_metricinfo->GetDerivFactors(GetPointsKeys());
533 
534  Array<OneD, NekDouble> tmp1 (nq); // Quad metric
535 
536  Array<OneD, NekDouble> tmp2 = outarray[0]; // Dir1 metric
537  Array<OneD, NekDouble> tmp3 = outarray[1]; // Dir2 metric
538  Array<OneD, NekDouble> tmp4 = outarray[2];
539 
540  Vmath::Vcopy(nq,inarray,1,tmp1,1); // Dir3 metric
541 
542  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
543  {
544  Vmath::Vmul(nq,&df[3*dir][0], 1,tmp1.get(),1,tmp2.get(),1);
545  Vmath::Vmul(nq,&df[3*dir+1][0],1,tmp1.get(),1,tmp3.get(),1);
546  Vmath::Vmul(nq,&df[3*dir+2][0],1,tmp1.get(),1,tmp4.get(),1);
547  }
548  else
549  {
550  Vmath::Smul(nq, df[3*dir][0], tmp1.get(),1,tmp2.get(), 1);
551  Vmath::Smul(nq, df[3*dir+1][0],tmp1.get(),1,tmp3.get(), 1);
552  Vmath::Smul(nq, df[3*dir+2][0],tmp1.get(),1,tmp4.get(), 1);
553  }
554  }
555 
556 
557  void HexExp::IProductWRTDerivBase_MatOp(
558  const int dir,
559  const Array<OneD, const NekDouble>& inarray,
560  Array<OneD, NekDouble> &outarray)
561  {
562  int nq = GetTotPoints();
563  StdRegions::MatrixType mtype = StdRegions::eIProductWRTDerivBase0;
564 
565  switch(dir)
566  {
567  case 0:
568  {
569  mtype = StdRegions::eIProductWRTDerivBase0;
570  }
571  break;
572  case 1:
573  {
574  mtype = StdRegions::eIProductWRTDerivBase1;
575  }
576  break;
577  case 2:
578  {
579  mtype = StdRegions::eIProductWRTDerivBase2;
580  }
581  break;
582  default:
583  {
584  ASSERTL1(false,"input dir is out of range");
585  }
586  break;
587  }
588 
589  MatrixKey iprodmatkey(mtype,DetShapeType(),*this);
590  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
591 
592  Blas::Dgemv('N',m_ncoeffs,nq,iprodmat->Scale(),(iprodmat->GetOwnedMatrix())->GetPtr().get(),
593  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
594  }
595 
596 
597  /**
598  *
599  * @param dir Vector direction in which to take the derivative.
600  * @param inarray The function \f$ u \f$.
601  * @param outarray Value of the inner product.
602  */
603  void HexExp::IProductWRTDirectionalDerivBase_SumFac(
604  const Array<OneD, const NekDouble>& direction,
605  const Array<OneD, const NekDouble>& inarray,
606  Array<OneD, NekDouble> & outarray)
607  {
608  int shapedim = 3;
609  const int nq0 = m_base[0]->GetNumPoints();
610  const int nq1 = m_base[1]->GetNumPoints();
611  const int nq2 = m_base[2]->GetNumPoints();
612  const int nq = nq0*nq1*nq2;
613  const int nm0 = m_base[0]->GetNumModes();
614  const int nm1 = m_base[1]->GetNumModes();
615 
616  const Array<TwoD, const NekDouble>& df =
617  m_metricinfo->GetDerivFactors(GetPointsKeys());
618 
619  Array<OneD, NekDouble> alloc(4*nq + m_ncoeffs + nm0*nq2*(nq1+nm1));
620  Array<OneD, NekDouble> tmp1 (alloc); // Quad metric
621  Array<OneD, NekDouble> tmp2 (alloc + nq); // Dir1 metric
622  Array<OneD, NekDouble> tmp3 (alloc + 2*nq); // Dir2 metric
623  Array<OneD, NekDouble> tmp4 (alloc + 3*nq); // Dir3 metric
624  Array<OneD, NekDouble> tmp5 (alloc + 4*nq); // iprod tmp
625  Array<OneD, NekDouble> wsp (tmp5 + m_ncoeffs); // Wsp
626 
627  MultiplyByQuadratureMetric(inarray, tmp1);
628 
629  Array<OneD, Array<OneD, NekDouble> > dfdir(shapedim);
630  Expansion::ComputeGmatcdotMF(df,direction,dfdir);
631 
632  Vmath::Vmul(nq,&dfdir[0][0],1,tmp1.get(),1,tmp2.get(),1);
633  Vmath::Vmul(nq,&dfdir[1][0],1,tmp1.get(),1,tmp3.get(),1);
634  Vmath::Vmul(nq,&dfdir[2][0],1,tmp1.get(),1,tmp4.get(),1);
635 
636  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
637  m_base[1]->GetBdata(),
638  m_base[2]->GetBdata(),
639  tmp2,outarray,wsp,
640  false,true,true);
641 
642  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
643  m_base[1]->GetDbdata(),
644  m_base[2]->GetBdata(),
645  tmp3,tmp5,wsp,
646  true,false,true);
647 
648  Vmath::Vadd(m_ncoeffs, tmp5, 1, outarray, 1, outarray, 1);
649 
650  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
651  m_base[1]->GetBdata(),
652  m_base[2]->GetDbdata(),
653  tmp4,tmp5,wsp,
654  true,true,false);
655 
656  Vmath::Vadd(m_ncoeffs, tmp5, 1, outarray, 1, outarray, 1);
657  }
658 
659 
660  //-----------------------------
661  // Evaluation functions
662  //-----------------------------
663 
664 
665  /**
666  * Given the local cartesian coordinate \a Lcoord evaluate the
667  * value of physvals at this point by calling through to the
668  * StdExpansion method
669  */
670  NekDouble HexExp::v_StdPhysEvaluate(
671  const Array<OneD, const NekDouble> &Lcoord,
672  const Array<OneD, const NekDouble> &physvals)
673  {
674  // Evaluate point in local coordinates.
675  return StdHexExp::v_PhysEvaluate(Lcoord,physvals);
676  }
677 
678  NekDouble HexExp::v_PhysEvaluate(
679  const Array<OneD, const NekDouble> &coord,
680  const Array<OneD, const NekDouble> & physvals)
681  {
682  Array<OneD,NekDouble> Lcoord = Array<OneD,NekDouble>(3);
683 
684  ASSERTL0(m_geom,"m_geom not defined");
685  m_geom->GetLocCoords(coord,Lcoord);
686  return StdHexExp::v_PhysEvaluate(Lcoord, physvals);
687  }
688 
689  StdRegions::StdExpansionSharedPtr HexExp::v_GetStdExp(void) const
690  {
691  return MemoryManager<StdRegions::StdHexExp>
692  ::AllocateSharedPtr(m_base[0]->GetBasisKey(),
693  m_base[1]->GetBasisKey(),
694  m_base[2]->GetBasisKey());
695  }
696 
697 
698  StdRegions::StdExpansionSharedPtr HexExp::v_GetLinStdExp(void) const
699  {
700  LibUtilities::BasisKey bkey0(m_base[0]->GetBasisType(),
701  2, m_base[0]->GetPointsKey());
702  LibUtilities::BasisKey bkey1(m_base[1]->GetBasisType(),
703  2, m_base[1]->GetPointsKey());
704  LibUtilities::BasisKey bkey2(m_base[2]->GetBasisType(),
705  2, m_base[2]->GetPointsKey());
706 
707  return MemoryManager<StdRegions::StdHexExp>
708  ::AllocateSharedPtr( bkey0, bkey1, bkey2);
709  }
710 
711  /**
712  * \brief Retrieves the physical coordinates of a given set of
713  * reference coordinates.
714  *
715  * @param Lcoords Local coordinates in reference space.
716  * @param coords Corresponding coordinates in physical space.
717  */
718  void HexExp::v_GetCoord(
719  const Array<OneD, const NekDouble> &Lcoords,
720  Array<OneD,NekDouble> &coords)
721  {
722  int i;
723 
724  ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[0] <= 1.0 &&
725  Lcoords[1] >= -1.0 && Lcoords[1] <= 1.0 &&
726  Lcoords[2] >= -1.0 && Lcoords[2] <= 1.0,
727  "Local coordinates are not in region [-1,1]");
728 
729  m_geom->FillGeom();
730 
731  for(i = 0; i < m_geom->GetCoordim(); ++i)
732  {
733  coords[i] = m_geom->GetCoord(i,Lcoords);
734  }
735  }
736 
737  void HexExp::v_GetCoords(
738  Array<OneD, NekDouble> &coords_0,
739  Array<OneD, NekDouble> &coords_1,
740  Array<OneD, NekDouble> &coords_2)
741  {
742  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
743  }
744 
745  //-----------------------------
746  // Helper functions
747  //-----------------------------
748 
749  /// Return the region shape using the enum-list of ShapeType
750  LibUtilities::ShapeType HexExp::v_DetShapeType() const
751  {
752  return LibUtilities::eHexahedron;
753  }
754 
755 
756  void HexExp::v_ExtractDataToCoeffs(
757  const NekDouble *data,
758  const std::vector<unsigned int > &nummodes,
759  const int mode_offset,
760  NekDouble * coeffs,
761  std::vector<LibUtilities::BasisType> &fromType)
762  {
763  int data_order0 = nummodes[mode_offset];
764  int fillorder0 = min(m_base[0]->GetNumModes(),data_order0);
765  int data_order1 = nummodes[mode_offset+1];
766  int order1 = m_base[1]->GetNumModes();
767  int fillorder1 = min(order1,data_order1);
768  int data_order2 = nummodes[mode_offset+2];
769  int order2 = m_base[2]->GetNumModes();
770  int fillorder2 = min(order2,data_order2);
771 
772  // Check if same basis
773  if (fromType[0] != m_base[0]->GetBasisType() ||
774  fromType[1] != m_base[1]->GetBasisType() ||
775  fromType[2] != m_base[2]->GetBasisType())
776  {
777  // Construct a hex with the appropriate basis type at our
778  // quadrature points, and one more to do a forwards
779  // transform. We can then copy the output to coeffs.
780  StdRegions::StdHexExp tmpHex(
781  LibUtilities::BasisKey(
782  fromType[0], data_order0, m_base[0]->GetPointsKey()),
783  LibUtilities::BasisKey(
784  fromType[1], data_order1, m_base[1]->GetPointsKey()),
785  LibUtilities::BasisKey(
786  fromType[2], data_order2, m_base[2]->GetPointsKey()));
787  StdRegions::StdHexExp tmpHex2(m_base[0]->GetBasisKey(),
788  m_base[1]->GetBasisKey(),
789  m_base[2]->GetBasisKey());
790 
791  Array<OneD, const NekDouble> tmpData(tmpHex.GetNcoeffs(), data);
792  Array<OneD, NekDouble> tmpBwd(tmpHex2.GetTotPoints());
793  Array<OneD, NekDouble> tmpOut(tmpHex2.GetNcoeffs());
794 
795  tmpHex.BwdTrans(tmpData, tmpBwd);
796  tmpHex2.FwdTrans(tmpBwd, tmpOut);
797  Vmath::Vcopy(tmpOut.size(), &tmpOut[0], 1, coeffs, 1);
798 
799  return;
800  }
801 
802  switch(m_base[0]->GetBasisType())
803  {
804  case LibUtilities::eModified_A:
805  {
806  int i,j;
807  int cnt = 0;
808  int cnt1 = 0;
809 
810  ASSERTL1(m_base[1]->GetBasisType() ==
811  LibUtilities::eModified_A,
812  "Extraction routine not set up for this basis");
813  ASSERTL1(m_base[2]->GetBasisType() ==
814  LibUtilities::eModified_A,
815  "Extraction routine not set up for this basis");
816 
817  Vmath::Zero(m_ncoeffs,coeffs,1);
818  for(j = 0; j < fillorder0; ++j)
819  {
820  for(i = 0; i < fillorder1; ++i)
821  {
822  Vmath::Vcopy(fillorder2, &data[cnt], 1,
823  &coeffs[cnt1], 1);
824  cnt += data_order2;
825  cnt1 += order2;
826  }
827 
828  // count out data for j iteration
829  for(i = fillorder1; i < data_order1; ++i)
830  {
831  cnt += data_order2;
832  }
833 
834  for(i = fillorder1; i < order1; ++i)
835  {
836  cnt1 += order2;
837  }
838  }
839  break;
840  }
841  case LibUtilities::eGLL_Lagrange:
842  {
843  LibUtilities::PointsKey
844  p0(nummodes[0], LibUtilities::eGaussLobattoLegendre);
845  LibUtilities::PointsKey
846  p1(nummodes[1], LibUtilities::eGaussLobattoLegendre);
847  LibUtilities::PointsKey
848  p2(nummodes[2], LibUtilities::eGaussLobattoLegendre);
849  LibUtilities::PointsKey t0(
850  m_base[0]->GetNumModes(),
851  LibUtilities::eGaussLobattoLegendre);
852  LibUtilities::PointsKey t1(
853  m_base[1]->GetNumModes(),
854  LibUtilities::eGaussLobattoLegendre);
855  LibUtilities::PointsKey t2(
856  m_base[2]->GetNumModes(),
857  LibUtilities::eGaussLobattoLegendre);
858  LibUtilities::Interp3D(p0, p1, p2, data, t0, t1, t2, coeffs);
859  }
860  break;
861  default:
862  ASSERTL0(false, "basis is either not set up or not "
863  "hierarchicial");
864  }
865  }
866 
867  bool HexExp::v_GetFaceDGForwards(const int i) const
868  {
869  StdRegions::Orientation fo = GetGeom3D()->GetForient(i);
870 
871  return fo == StdRegions::eDir1FwdDir1_Dir2FwdDir2 ||
872  fo == StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
873  fo == StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
874  fo == StdRegions::eDir1FwdDir2_Dir2BwdDir1;
875  }
876 
877  void HexExp::v_GetTracePhysMap(const int face,
878  Array<OneD, int> &outarray)
879  {
880  int nquad0 = m_base[0]->GetNumPoints();
881  int nquad1 = m_base[1]->GetNumPoints();
882  int nquad2 = m_base[2]->GetNumPoints();
883 
884  int nq0 = 0;
885  int nq1 = 0;
886 
887  switch(face)
888  {
889  case 0:
890  nq0 = nquad0;
891  nq1 = nquad1;
892 
893  //Directions A and B positive
894  if(outarray.size()!=nq0*nq1)
895  {
896  outarray = Array<OneD, int>(nq0*nq1);
897  }
898 
899  for (int i = 0; i < nquad0*nquad1; ++i)
900  {
901  outarray[i] = i;
902  }
903 
904  break;
905  case 1:
906  nq0 = nquad0;
907  nq1 = nquad2;
908  //Direction A and B positive
909  if(outarray.size()!=nq0*nq1)
910  {
911  outarray = Array<OneD, int>(nq0*nq1);
912  }
913 
914  //Direction A and B positive
915  for (int k = 0; k < nquad2; k++)
916  {
917  for(int i = 0; i < nquad0; ++i)
918  {
919  outarray[k*nquad0 + i] = nquad0*nquad1*k + i;
920  }
921  }
922  break;
923  case 2:
924  nq0 = nquad1;
925  nq1 = nquad2;
926 
927  //Direction A and B positive
928  if(outarray.size()!=nq0*nq1)
929  {
930  outarray = Array<OneD, int>(nq0*nq1);
931  }
932 
933  for (int i = 0; i < nquad1*nquad2; i++)
934  {
935  outarray[i] = nquad0-1 + i*nquad0;
936  }
937  break;
938  case 3:
939  nq0 = nquad0;
940  nq1 = nquad2;
941 
942  //Direction A and B positive
943  if(outarray.size()!=nq0*nq1)
944  {
945  outarray = Array<OneD, int>(nq0*nq1);
946  }
947 
948  for (int k = 0; k < nquad2; k++)
949  {
950  for (int i = 0; i < nquad0; i++)
951  {
952  outarray[k*nquad0 + i] = (nquad0*(nquad1-1))+(k*nquad0*nquad1) + i;
953  }
954  }
955  break;
956  case 4:
957  nq0 = nquad1;
958  nq1 = nquad2;
959 
960  //Direction A and B positive
961  if(outarray.size()!=nq0*nq1)
962  {
963  outarray = Array<OneD, int>(nq0*nq1);
964  }
965 
966  for (int i = 0; i < nquad1*nquad2; i++)
967  {
968  outarray[i] = i*nquad0;
969  }
970  break;
971  case 5:
972  nq0 = nquad0;
973  nq1 = nquad1;
974  //Directions A and B positive
975  if(outarray.size()!=nq0*nq1)
976  {
977  outarray = Array<OneD, int>(nq0*nq1);
978  }
979 
980  for (int i = 0; i < nquad0*nquad1; i++)
981  {
982  outarray[i] = nquad0*nquad1*(nquad2-1) + i;
983  }
984 
985  break;
986  default:
987  ASSERTL0(false,"face value (> 5) is out of range");
988  break;
989  }
990 
991  }
992 
993  void HexExp::v_ComputeTraceNormal(const int face)
994  {
995  int i;
996  const SpatialDomains::GeomFactorsSharedPtr & geomFactors =
997  GetGeom()->GetMetricInfo();
998  SpatialDomains::GeomType type = geomFactors->GetGtype();
999 
1000  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1001  for(i = 0; i < ptsKeys.size(); ++i)
1002  {
1003  // Need at least 2 points for computing normals
1004  if (ptsKeys[i].GetNumPoints() == 1)
1005  {
1006  LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
1007  ptsKeys[i] = pKey;
1008  }
1009  }
1010 
1011  const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys);
1012  const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys);
1013 
1014  LibUtilities::BasisKey tobasis0 = GetTraceBasisKey(face,0);
1015  LibUtilities::BasisKey tobasis1 = GetTraceBasisKey(face,1);
1016 
1017  // Number of quadrature points in face expansion.
1018  int nq_face = tobasis0.GetNumPoints()*tobasis1.GetNumPoints();
1019 
1020  int vCoordDim = GetCoordim();
1021 
1022  m_faceNormals[face] = Array<OneD, Array<OneD, NekDouble> >(vCoordDim);
1023  Array<OneD, Array<OneD, NekDouble> > &normal = m_faceNormals[face];
1024  for (i = 0; i < vCoordDim; ++i)
1025  {
1026  normal[i] = Array<OneD, NekDouble>(nq_face);
1027  }
1028 
1029  size_t nqb = nq_face;
1030  size_t nbnd= face;
1031  m_elmtBndNormDirElmtLen[nbnd] = Array<OneD, NekDouble> {nqb, 0.0};
1032  Array<OneD, NekDouble> &length = m_elmtBndNormDirElmtLen[nbnd];
1033 
1034  // Regular geometry case
1035  if((type == SpatialDomains::eRegular)||(type == SpatialDomains::eMovingRegular))
1036  {
1037  NekDouble fac;
1038  // Set up normals
1039  switch(face)
1040  {
1041  case 0:
1042  for(i = 0; i < vCoordDim; ++i)
1043  {
1044  normal[i][0] = -df[3*i+2][0];
1045  }
1046  break;
1047  case 1:
1048  for(i = 0; i < vCoordDim; ++i)
1049  {
1050  normal[i][0] = -df[3*i+1][0];
1051  }
1052  break;
1053  case 2:
1054  for(i = 0; i < vCoordDim; ++i)
1055  {
1056  normal[i][0] = df[3*i][0];
1057  }
1058  break;
1059  case 3:
1060  for(i = 0; i < vCoordDim; ++i)
1061  {
1062  normal[i][0] = df[3*i+1][0];
1063  }
1064  break;
1065  case 4:
1066  for(i = 0; i < vCoordDim; ++i)
1067  {
1068  normal[i][0] = -df[3*i][0];
1069  }
1070  break;
1071  case 5:
1072  for(i = 0; i < vCoordDim; ++i)
1073  {
1074  normal[i][0] = df[3*i+2][0];
1075  }
1076  break;
1077  default:
1078  ASSERTL0(false,"face is out of range (edge < 5)");
1079  }
1080 
1081  // normalise
1082  fac = 0.0;
1083  for(i =0 ; i < vCoordDim; ++i)
1084  {
1085  fac += normal[i][0]*normal[i][0];
1086  }
1087  fac = 1.0/sqrt(fac);
1088 
1089  Vmath::Fill(nqb, fac, length, 1);
1090  for (i = 0; i < vCoordDim; ++i)
1091  {
1092  Vmath::Fill(nq_face, fac*normal[i][0], normal[i], 1);
1093  }
1094 
1095  }
1096  else // Set up deformed normals
1097  {
1098  int j, k;
1099 
1100  int nqe0 = ptsKeys[0].GetNumPoints();
1101  int nqe1 = ptsKeys[0].GetNumPoints();
1102  int nqe2 = ptsKeys[0].GetNumPoints();
1103  int nqe01 = nqe0*nqe1;
1104  int nqe02 = nqe0*nqe2;
1105  int nqe12 = nqe1*nqe2;
1106 
1107  int nqe;
1108  if (face == 0 || face == 5)
1109  {
1110  nqe = nqe01;
1111  }
1112  else if (face == 1 || face == 3)
1113  {
1114  nqe = nqe02;
1115  }
1116  else
1117  {
1118  nqe = nqe12;
1119  }
1120 
1121  LibUtilities::PointsKey points0;
1122  LibUtilities::PointsKey points1;
1123 
1124  Array<OneD, NekDouble> faceJac(nqe);
1125  Array<OneD, NekDouble> normals(vCoordDim*nqe,0.0);
1126 
1127  // Extract Jacobian along face and recover local
1128  // derivates (dx/dr) for polynomial interpolation by
1129  // multiplying m_gmat by jacobian
1130  switch(face)
1131  {
1132  case 0:
1133  for(j = 0; j < nqe; ++j)
1134  {
1135  normals[j] = -df[2][j]*jac[j];
1136  normals[nqe+j] = -df[5][j]*jac[j];
1137  normals[2*nqe+j] = -df[8][j]*jac[j];
1138  faceJac[j] = jac[j];
1139  }
1140 
1141  points0 = ptsKeys[0];
1142  points1 = ptsKeys[1];
1143  break;
1144  case 1:
1145  for (j = 0; j < nqe0; ++j)
1146  {
1147  for(k = 0; k < nqe2; ++k)
1148  {
1149  int idx = j + nqe01*k;
1150  normals[j+k*nqe0] = -df[1][idx]*jac[idx];
1151  normals[nqe+j+k*nqe0] = -df[4][idx]*jac[idx];
1152  normals[2*nqe+j+k*nqe0] = -df[7][idx]*jac[idx];
1153  faceJac[j+k*nqe0] = jac[idx];
1154  }
1155  }
1156  points0 = ptsKeys[0];
1157  points1 = ptsKeys[2];
1158  break;
1159  case 2:
1160  for (j = 0; j < nqe1; ++j)
1161  {
1162  for(k = 0; k < nqe2; ++k)
1163  {
1164  int idx = nqe0-1+nqe0*j+nqe01*k;
1165  normals[j+k*nqe1] = df[0][idx]*jac[idx];
1166  normals[nqe+j+k*nqe1] = df[3][idx]*jac[idx];
1167  normals[2*nqe+j+k*nqe1] = df[6][idx]*jac[idx];
1168  faceJac[j+k*nqe1] = jac[idx];
1169  }
1170  }
1171  points0 = ptsKeys[1];
1172  points1 = ptsKeys[2];
1173  break;
1174  case 3:
1175  for (j = 0; j < nqe0; ++j)
1176  {
1177  for(k = 0; k < nqe2; ++k)
1178  {
1179  int idx = nqe0*(nqe1-1)+j+nqe01*k;
1180  normals[j+k*nqe0] = df[1][idx]*jac[idx];
1181  normals[nqe+j+k*nqe0] = df[4][idx]*jac[idx];
1182  normals[2*nqe+j+k*nqe0] = df[7][idx]*jac[idx];
1183  faceJac[j+k*nqe0] = jac[idx];
1184  }
1185  }
1186  points0 = ptsKeys[0];
1187  points1 = ptsKeys[2];
1188  break;
1189  case 4:
1190  for (j = 0; j < nqe1; ++j)
1191  {
1192  for(k = 0; k < nqe2; ++k)
1193  {
1194  int idx = j*nqe0+nqe01*k;
1195  normals[j+k*nqe1] = -df[0][idx]*jac[idx];
1196  normals[nqe+j+k*nqe1] = -df[3][idx]*jac[idx];
1197  normals[2*nqe+j+k*nqe1] = -df[6][idx]*jac[idx];
1198  faceJac[j+k*nqe1] = jac[idx];
1199  }
1200  }
1201  points0 = ptsKeys[1];
1202  points1 = ptsKeys[2];
1203  break;
1204  case 5:
1205  for (j = 0; j < nqe01; ++j)
1206  {
1207  int idx = j+nqe01*(nqe2-1);
1208  normals[j] = df[2][idx]*jac[idx];
1209  normals[nqe+j] = df[5][idx]*jac[idx];
1210  normals[2*nqe+j] = df[8][idx]*jac[idx];
1211  faceJac[j] = jac[idx];
1212  }
1213  points0 = ptsKeys[0];
1214  points1 = ptsKeys[1];
1215  break;
1216  default:
1217  ASSERTL0(false,"face is out of range (face < 5)");
1218  }
1219 
1220  Array<OneD, NekDouble> work (nq_face, 0.0);
1221  // Interpolate Jacobian and invert
1222  LibUtilities::Interp2D(points0, points1, faceJac,
1223  tobasis0.GetPointsKey(),
1224  tobasis1.GetPointsKey(),
1225  work);
1226 
1227  Vmath::Sdiv(nq_face,1.0,&work[0],1,&work[0],1);
1228 
1229  // interpolate
1230  for(i = 0; i < GetCoordim(); ++i)
1231  {
1232  LibUtilities::Interp2D(points0, points1,
1233  &normals[i*nqe],
1234  tobasis0.GetPointsKey(),
1235  tobasis1.GetPointsKey(),
1236  &normal[i][0]);
1237  Vmath::Vmul(nq_face,work,1,normal[i],1,normal[i],1);
1238  }
1239 
1240  //normalise normal vectors
1241  Vmath::Zero(nq_face,work,1);
1242  for(i = 0; i < GetCoordim(); ++i)
1243  {
1244  Vmath::Vvtvp(nq_face,normal[i],1, normal[i],1,work,1,work,1);
1245  }
1246 
1247  Vmath::Vsqrt(nq_face,work,1,work,1);
1248  Vmath::Sdiv(nq_face,1.0,work,1,work,1);
1249 
1250  Vmath::Vcopy(nqb, work, 1, length, 1);
1251 
1252  for(i = 0; i < GetCoordim(); ++i)
1253  {
1254  Vmath::Vmul(nq_face,normal[i],1,work,1,normal[i],1);
1255  }
1256  }
1257  }
1258 
1259  //-----------------------------
1260  // Operator creation functions
1261  //-----------------------------
1262  void HexExp::v_MassMatrixOp(
1263  const Array<OneD, const NekDouble> &inarray,
1264  Array<OneD,NekDouble> &outarray,
1265  const StdRegions::StdMatrixKey &mkey)
1266  {
1267  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
1268  }
1269 
1270  void HexExp::v_LaplacianMatrixOp(
1271  const Array<OneD, const NekDouble> &inarray,
1272  Array<OneD,NekDouble> &outarray,
1273  const StdRegions::StdMatrixKey &mkey)
1274  {
1275  HexExp::v_LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1276  }
1277 
1278  void HexExp::v_LaplacianMatrixOp(
1279  const int k1,
1280  const int k2,
1281  const Array<OneD, const NekDouble> &inarray,
1282  Array<OneD,NekDouble> &outarray,
1283  const StdRegions::StdMatrixKey &mkey)
1284  {
1285  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,
1286  mkey);
1287  }
1288 
1289  void HexExp::v_WeakDerivMatrixOp(
1290  const int i,
1291  const Array<OneD, const NekDouble> &inarray,
1292  Array<OneD,NekDouble> &outarray,
1293  const StdRegions::StdMatrixKey &mkey)
1294  {
1295  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
1296  }
1297 
1298  void HexExp::v_WeakDirectionalDerivMatrixOp(
1299  const Array<OneD, const NekDouble> &inarray,
1300  Array<OneD,NekDouble> &outarray,
1301  const StdRegions::StdMatrixKey &mkey)
1302  {
1303  StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray,
1304  outarray,mkey);
1305  }
1306 
1307  void HexExp::v_MassLevelCurvatureMatrixOp(
1308  const Array<OneD, const NekDouble> &inarray,
1309  Array<OneD,NekDouble> &outarray,
1310  const StdRegions::StdMatrixKey &mkey)
1311  {
1312  StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray,
1313  outarray,mkey);
1314  }
1315 
1316  void HexExp::v_HelmholtzMatrixOp(
1317  const Array<OneD, const NekDouble> &inarray,
1318  Array<OneD,NekDouble> &outarray,
1319  const StdRegions::StdMatrixKey &mkey)
1320  {
1321  HexExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1322  }
1323 
1324 
1325  void HexExp::v_GeneralMatrixOp_MatOp(
1326  const Array<OneD, const NekDouble> &inarray,
1327  Array<OneD,NekDouble> &outarray,
1328  const StdRegions::StdMatrixKey &mkey)
1329  {
1330  //int nConsts = mkey.GetNconstants();
1331  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1332 
1333 // switch(nConsts)
1334 // {
1335 // case 0:
1336 // {
1337 // mat = GetLocMatrix(mkey.GetMatrixType());
1338 // }
1339 // break;
1340 // case 1:
1341 // {
1342 // mat = GetLocMatrix(mkey.GetMatrixType(),mkey.GetConstant(0));
1343 // }
1344 // break;
1345 // case 2:
1346 // {
1347 // mat = GetLocMatrix(mkey.GetMatrixType(),mkey.GetConstant(0),mkey.GetConstant(1));
1348 // }
1349 // break;
1350 //
1351 // default:
1352 // {
1353 // NEKERROR(ErrorUtil::efatal, "Unknown number of constants");
1354 // }
1355 // break;
1356 // }
1357 
1358  if(inarray.get() == outarray.get())
1359  {
1360  Array<OneD,NekDouble> tmp(m_ncoeffs);
1361  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1362 
1363  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1364  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1365  }
1366  else
1367  {
1368  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1369  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1370  }
1371  }
1372 
1373  /**
1374  * This function is used to compute exactly the advective numerical flux
1375  * on the interface of two elements with different expansions, hence an
1376  * appropriate number of Gauss points has to be used. The number of
1377  * Gauss points has to be equal to the number used by the highest
1378  * polynomial degree of the two adjacent elements
1379  *
1380  * @param numMin Is the reduced polynomial order
1381  * @param inarray Input array of coefficients
1382  * @param dumpVar Output array of reduced coefficients.
1383  */
1384  void HexExp::v_ReduceOrderCoeffs(
1385  int numMin,
1386  const Array<OneD, const NekDouble> &inarray,
1387  Array<OneD, NekDouble> &outarray)
1388  {
1389  int n_coeffs = inarray.size();
1390  int nmodes0 = m_base[0]->GetNumModes();
1391  int nmodes1 = m_base[1]->GetNumModes();
1392  int nmodes2 = m_base[2]->GetNumModes();
1393  int numMax = nmodes0;
1394 
1395  Array<OneD, NekDouble> coeff (n_coeffs);
1396  Array<OneD, NekDouble> coeff_tmp1(nmodes0*nmodes1, 0.0);
1397  Array<OneD, NekDouble> coeff_tmp2(n_coeffs, 0.0);
1398  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
1399 
1400  Vmath::Vcopy(n_coeffs,inarray,1,coeff_tmp2,1);
1401 
1402  const LibUtilities::PointsKey Pkey0(
1403  nmodes0, LibUtilities::eGaussLobattoLegendre);
1404  const LibUtilities::PointsKey Pkey1(
1405  nmodes1, LibUtilities::eGaussLobattoLegendre);
1406  const LibUtilities::PointsKey Pkey2(
1407  nmodes2, LibUtilities::eGaussLobattoLegendre);
1408 
1409  LibUtilities::BasisKey b0(
1410  m_base[0]->GetBasisType(), nmodes0, Pkey0);
1411  LibUtilities::BasisKey b1(
1412  m_base[1]->GetBasisType(), nmodes1, Pkey1);
1413  LibUtilities::BasisKey b2(
1414  m_base[2]->GetBasisType(), nmodes2, Pkey2);
1415  LibUtilities::BasisKey bortho0(
1416  LibUtilities::eOrtho_A, nmodes0, Pkey0);
1417  LibUtilities::BasisKey bortho1(
1418  LibUtilities::eOrtho_A, nmodes1, Pkey1);
1419  LibUtilities::BasisKey bortho2(
1420  LibUtilities::eOrtho_A, nmodes2, Pkey2);
1421 
1422  LibUtilities::InterpCoeff3D(
1423  b0, b1, b2, coeff_tmp2,
1424  bortho0, bortho1, bortho2, coeff);
1425 
1426  Vmath::Zero(n_coeffs, coeff_tmp2, 1);
1427 
1428  int cnt = 0, cnt2 = 0;
1429 
1430  for (int u = 0; u < numMin+1; ++u)
1431  {
1432  for (int i = 0; i < numMin; ++i)
1433  {
1434  Vmath::Vcopy(numMin,
1435  tmp = coeff+cnt+cnt2,1,
1436  tmp2 = coeff_tmp1+cnt,1);
1437 
1438  cnt = i*numMax;
1439  }
1440 
1441  Vmath::Vcopy(nmodes0*nmodes1,
1442  tmp3 = coeff_tmp1,1,
1443  tmp4 = coeff_tmp2+cnt2,1);
1444 
1445  cnt2 = u*nmodes0*nmodes1;
1446  }
1447 
1448  LibUtilities::InterpCoeff3D(
1449  bortho0, bortho1, bortho2, coeff_tmp2,
1450  b0, b1, b2, outarray);
1451  }
1452 
1453  void HexExp::v_SVVLaplacianFilter(
1454  Array<OneD, NekDouble> &array,
1455  const StdRegions::StdMatrixKey &mkey)
1456  {
1457  int nq = GetTotPoints();
1458 
1459  // Calculate sqrt of the Jacobian
1460  Array<OneD, const NekDouble> jac =
1461  m_metricinfo->GetJac(GetPointsKeys());
1462  Array<OneD, NekDouble> sqrt_jac(nq);
1463  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1464  {
1465  Vmath::Vsqrt(nq,jac,1,sqrt_jac,1);
1466  }
1467  else
1468  {
1469  Vmath::Fill(nq,sqrt(jac[0]),sqrt_jac,1);
1470  }
1471 
1472  // Multiply array by sqrt(Jac)
1473  Vmath::Vmul(nq,sqrt_jac,1,array,1,array,1);
1474 
1475  // Apply std region filter
1476  StdHexExp::v_SVVLaplacianFilter( array, mkey);
1477 
1478  // Divide by sqrt(Jac)
1479  Vmath::Vdiv(nq,array,1,sqrt_jac,1,array,1);
1480  }
1481 
1482  //-----------------------------
1483  // Matrix creation functions
1484  //-----------------------------
1485  DNekMatSharedPtr HexExp::v_GenMatrix(
1486  const StdRegions::StdMatrixKey &mkey)
1487  {
1488  DNekMatSharedPtr returnval;
1489 
1490  switch(mkey.GetMatrixType())
1491  {
1492  case StdRegions::eHybridDGHelmholtz:
1493  case StdRegions::eHybridDGLamToU:
1494  case StdRegions::eHybridDGLamToQ0:
1495  case StdRegions::eHybridDGLamToQ1:
1496  case StdRegions::eHybridDGLamToQ2:
1497  case StdRegions::eHybridDGHelmBndLam:
1498  case StdRegions::eInvLaplacianWithUnityMean:
1499  returnval = Expansion3D::v_GenMatrix(mkey);
1500  break;
1501  default:
1502  returnval = StdHexExp::v_GenMatrix(mkey);
1503  }
1504 
1505  return returnval;
1506  }
1507 
1508 
1509  DNekMatSharedPtr HexExp::v_CreateStdMatrix(
1510  const StdRegions::StdMatrixKey &mkey)
1511  {
1512  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1513  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1514  LibUtilities::BasisKey bkey2 = m_base[2]->GetBasisKey();
1515 
1516  StdRegions::StdHexExpSharedPtr tmp = MemoryManager<StdHexExp>
1517  ::AllocateSharedPtr(bkey0,bkey1,bkey2);
1518 
1519  return tmp->GetStdMatrix(mkey);
1520  }
1521 
1522 
1523  DNekScalMatSharedPtr HexExp::CreateMatrix(const MatrixKey &mkey)
1524  {
1525  DNekScalMatSharedPtr returnval;
1526  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1527 
1528  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,
1529  "Geometric information is not set up");
1530 
1531  switch(mkey.GetMatrixType())
1532  {
1533  case StdRegions::eMass:
1534  {
1535  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1536  mkey.GetNVarCoeff())
1537  {
1538  NekDouble one = 1.0;
1539  DNekMatSharedPtr mat = GenMatrix(mkey);
1540  returnval = MemoryManager<DNekScalMat>
1541  ::AllocateSharedPtr(one,mat);
1542  }
1543  else
1544  {
1545  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1546  DNekMatSharedPtr mat
1547  = GetStdMatrix(mkey);
1548  returnval = MemoryManager<DNekScalMat>
1549  ::AllocateSharedPtr(jac,mat);
1550  }
1551  }
1552  break;
1553  case StdRegions::eInvMass:
1554  {
1555  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1556  {
1557  NekDouble one = 1.0;
1558  StdRegions::StdMatrixKey masskey(StdRegions::eMass,
1559  DetShapeType(), *this);
1560  DNekMatSharedPtr mat = GenMatrix(masskey);
1561  mat->Invert();
1562 
1563  returnval = MemoryManager<DNekScalMat>
1564  ::AllocateSharedPtr(one,mat);
1565  }
1566  else
1567  {
1568  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1569  DNekMatSharedPtr mat
1570  = GetStdMatrix(mkey);
1571  returnval = MemoryManager<DNekScalMat>
1572  ::AllocateSharedPtr(fac,mat);
1573  }
1574  }
1575  break;
1576  case StdRegions::eWeakDeriv0:
1577  case StdRegions::eWeakDeriv1:
1578  case StdRegions::eWeakDeriv2:
1579  {
1580  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1581  mkey.GetNVarCoeff())
1582  {
1583  NekDouble one = 1.0;
1584  DNekMatSharedPtr mat = GenMatrix(mkey);
1585 
1586  returnval = MemoryManager<DNekScalMat>
1587  ::AllocateSharedPtr(one,mat);
1588  }
1589  else
1590  {
1591  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1592  Array<TwoD, const NekDouble> df
1593  = m_metricinfo->GetDerivFactors(ptsKeys);
1594  int dir = 0;
1595 
1596  switch(mkey.GetMatrixType())
1597  {
1598  case StdRegions::eWeakDeriv0:
1599  dir = 0;
1600  break;
1601  case StdRegions::eWeakDeriv1:
1602  dir = 1;
1603  break;
1604  case StdRegions::eWeakDeriv2:
1605  dir = 2;
1606  break;
1607  default:
1608  break;
1609  }
1610 
1611  MatrixKey deriv0key(StdRegions::eWeakDeriv0,
1612  mkey.GetShapeType(), *this);
1613  MatrixKey deriv1key(StdRegions::eWeakDeriv1,
1614  mkey.GetShapeType(), *this);
1615  MatrixKey deriv2key(StdRegions::eWeakDeriv2,
1616  mkey.GetShapeType(), *this);
1617 
1618  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
1619  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
1620  DNekMat &deriv2 = *GetStdMatrix(deriv2key);
1621 
1622  int rows = deriv0.GetRows();
1623  int cols = deriv1.GetColumns();
1624 
1625  DNekMatSharedPtr WeakDeriv = MemoryManager<DNekMat>
1626  ::AllocateSharedPtr(rows,cols);
1627 
1628  (*WeakDeriv) = df[3*dir ][0]*deriv0
1629  + df[3*dir+1][0]*deriv1
1630  + df[3*dir+2][0]*deriv2;
1631 
1632  returnval = MemoryManager<DNekScalMat>
1633  ::AllocateSharedPtr(jac,WeakDeriv);
1634  }
1635  }
1636  break;
1637  case StdRegions::eLaplacian:
1638  {
1639  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1640  mkey.GetNVarCoeff()||
1641  mkey.ConstFactorExists(
1642  StdRegions::eFactorSVVCutoffRatio))
1643  {
1644  NekDouble one = 1.0;
1645  DNekMatSharedPtr mat = GenMatrix(mkey);
1646 
1647  returnval = MemoryManager<DNekScalMat>
1648  ::AllocateSharedPtr(one,mat);
1649  }
1650  else
1651  {
1652  MatrixKey lap00key(StdRegions::eLaplacian00,
1653  mkey.GetShapeType(), *this);
1654  MatrixKey lap01key(StdRegions::eLaplacian01,
1655  mkey.GetShapeType(), *this);
1656  MatrixKey lap02key(StdRegions::eLaplacian02,
1657  mkey.GetShapeType(), *this);
1658  MatrixKey lap11key(StdRegions::eLaplacian11,
1659  mkey.GetShapeType(), *this);
1660  MatrixKey lap12key(StdRegions::eLaplacian12,
1661  mkey.GetShapeType(), *this);
1662  MatrixKey lap22key(StdRegions::eLaplacian22,
1663  mkey.GetShapeType(), *this);
1664 
1665  DNekMat &lap00 = *GetStdMatrix(lap00key);
1666  DNekMat &lap01 = *GetStdMatrix(lap01key);
1667  DNekMat &lap02 = *GetStdMatrix(lap02key);
1668  DNekMat &lap11 = *GetStdMatrix(lap11key);
1669  DNekMat &lap12 = *GetStdMatrix(lap12key);
1670  DNekMat &lap22 = *GetStdMatrix(lap22key);
1671 
1672  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1673  Array<TwoD, const NekDouble> gmat
1674  = m_metricinfo->GetGmat(ptsKeys);
1675 
1676  int rows = lap00.GetRows();
1677  int cols = lap00.GetColumns();
1678 
1679  DNekMatSharedPtr lap = MemoryManager<DNekMat>
1680  ::AllocateSharedPtr(rows,cols);
1681 
1682  (*lap) = gmat[0][0]*lap00
1683  + gmat[4][0]*lap11
1684  + gmat[8][0]*lap22
1685  + gmat[3][0]*(lap01 + Transpose(lap01))
1686  + gmat[6][0]*(lap02 + Transpose(lap02))
1687  + gmat[7][0]*(lap12 + Transpose(lap12));
1688 
1689  returnval = MemoryManager<DNekScalMat>
1690  ::AllocateSharedPtr(jac,lap);
1691  }
1692  }
1693  break;
1694  case StdRegions::eHelmholtz:
1695  {
1696  NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
1697  MatrixKey masskey(StdRegions::eMass,
1698  mkey.GetShapeType(), *this);
1699  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1700  MatrixKey lapkey(StdRegions::eLaplacian,
1701  mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1702  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1703 
1704  int rows = LapMat.GetRows();
1705  int cols = LapMat.GetColumns();
1706 
1707  DNekMatSharedPtr helm = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1708 
1709  NekDouble one = 1.0;
1710  (*helm) = LapMat + lambda*MassMat;
1711 
1712  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
1713  }
1714  break;
1715  case StdRegions::eIProductWRTBase:
1716  {
1717  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1718  {
1719  NekDouble one = 1.0;
1720  DNekMatSharedPtr mat = GenMatrix(mkey);
1721  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1722  }
1723  else
1724  {
1725  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1726  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1727  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1728  }
1729  }
1730  break;
1731  case StdRegions::eHybridDGHelmholtz:
1732  case StdRegions::eHybridDGLamToU:
1733  case StdRegions::eHybridDGLamToQ0:
1734  case StdRegions::eHybridDGLamToQ1:
1735  case StdRegions::eHybridDGLamToQ2:
1736  case StdRegions::eHybridDGHelmBndLam:
1737  case StdRegions::eInvLaplacianWithUnityMean:
1738  {
1739  NekDouble one = 1.0;
1740 
1741  DNekMatSharedPtr mat = GenMatrix(mkey);
1742  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1743  }
1744  break;
1745  case StdRegions::eInvHybridDGHelmholtz:
1746  {
1747  NekDouble one = 1.0;
1748 
1749 // StdRegions::StdMatrixKey hkey(StdRegions::eHybridDGHelmholtz,
1750 // DetShapeType(),*this,
1751 // mkey.GetConstant(0),
1752 // mkey.GetConstant(1));
1753  MatrixKey hkey(StdRegions::eHybridDGHelmholtz, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1754  DNekMatSharedPtr mat = GenMatrix(hkey);
1755 
1756  mat->Invert();
1757  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1758  }
1759  break;
1760  case StdRegions::ePreconLinearSpace:
1761  {
1762  NekDouble one = 1.0;
1763  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1764  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1765  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1766  DNekMatSharedPtr R=BuildVertexMatrix(A);
1767 
1768  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1769  }
1770  break;
1771  case StdRegions::ePreconLinearSpaceMass:
1772  {
1773  NekDouble one = 1.0;
1774  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1775  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1776  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1777  DNekMatSharedPtr R=BuildVertexMatrix(A);
1778 
1779  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1780  }
1781  break;
1782  case StdRegions::ePreconR:
1783  {
1784  NekDouble one = 1.0;
1785  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this,mkey.GetConstFactors(), mkey.GetVarCoeffs());
1786  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1787  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1788 
1789  DNekScalMatSharedPtr Atmp;
1790  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1791 
1792  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1793  }
1794  break;
1795  case StdRegions::ePreconRMass:
1796  {
1797  NekDouble one = 1.0;
1798  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1799  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1800  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1801 
1802  DNekScalMatSharedPtr Atmp;
1803  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1804 
1805  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1806  }
1807  break;
1808  default:
1809  {
1810  NekDouble one = 1.0;
1811  DNekMatSharedPtr mat = GenMatrix(mkey);
1812 
1813  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1814  }
1815  break;
1816  }
1817 
1818  return returnval;
1819  }
1820 
1821 
1822  DNekScalBlkMatSharedPtr HexExp::CreateStaticCondMatrix(const MatrixKey &mkey)
1823  {
1824  DNekScalBlkMatSharedPtr returnval;
1825 
1826  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
1827 
1828  // set up block matrix system
1829  unsigned int nbdry = NumBndryCoeffs();
1830  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
1831  unsigned int exp_size[] = {nbdry,nint};
1832  unsigned int nblks = 2;
1833  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks,nblks,exp_size,exp_size); //Really need a constructor which takes Arrays
1834  NekDouble factor = 1.0;
1835 
1836  switch(mkey.GetMatrixType())
1837  {
1838  case StdRegions::eLaplacian:
1839  case StdRegions::eHelmholtz: // special case since Helmholtz not defined in StdRegions
1840 
1841  // use Deformed case for both regular and deformed geometries
1842  factor = 1.0;
1843  goto UseLocRegionsMatrix;
1844  break;
1845  default:
1846  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1847  mkey.GetNVarCoeff())
1848  {
1849  factor = 1.0;
1850  goto UseLocRegionsMatrix;
1851  }
1852  else
1853  {
1854  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1855  factor = mat->Scale();
1856  goto UseStdRegionsMatrix;
1857  }
1858  break;
1859  UseStdRegionsMatrix:
1860  {
1861  NekDouble invfactor = 1.0/factor;
1862  NekDouble one = 1.0;
1863  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
1864  DNekScalMatSharedPtr Atmp;
1865  DNekMatSharedPtr Asubmat;
1866 
1867  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
1868  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
1869  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
1870  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
1871  }
1872  break;
1873  UseLocRegionsMatrix:
1874  {
1875  int i,j;
1876  NekDouble invfactor = 1.0/factor;
1877  NekDouble one = 1.0;
1878  DNekScalMat &mat = *GetLocMatrix(mkey);
1879  DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
1880  DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
1881  DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
1882  DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
1883 
1884  Array<OneD,unsigned int> bmap(nbdry);
1885  Array<OneD,unsigned int> imap(nint);
1886  GetBoundaryMap(bmap);
1887  GetInteriorMap(imap);
1888 
1889  for(i = 0; i < nbdry; ++i)
1890  {
1891  for(j = 0; j < nbdry; ++j)
1892  {
1893  (*A)(i,j) = mat(bmap[i],bmap[j]);
1894  }
1895 
1896  for(j = 0; j < nint; ++j)
1897  {
1898  (*B)(i,j) = mat(bmap[i],imap[j]);
1899  }
1900  }
1901 
1902  for(i = 0; i < nint; ++i)
1903  {
1904  for(j = 0; j < nbdry; ++j)
1905  {
1906  (*C)(i,j) = mat(imap[i],bmap[j]);
1907  }
1908 
1909  for(j = 0; j < nint; ++j)
1910  {
1911  (*D)(i,j) = mat(imap[i],imap[j]);
1912  }
1913  }
1914 
1915  // Calculate static condensed system
1916  if(nint)
1917  {
1918  D->Invert();
1919  (*B) = (*B)*(*D);
1920  (*A) = (*A) - (*B)*(*C);
1921  }
1922 
1923  DNekScalMatSharedPtr Atmp;
1924 
1925  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
1926  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
1927  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
1928  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
1929 
1930  }
1931  }
1932  return returnval;
1933  }
1934 
1935 
1936  DNekScalMatSharedPtr HexExp::v_GetLocMatrix(const MatrixKey &mkey)
1937  {
1938  return m_matrixManager[mkey];
1939  }
1940 
1941 
1942  DNekScalBlkMatSharedPtr HexExp::v_GetLocStaticCondMatrix(
1943  const MatrixKey &mkey)
1944  {
1945  return m_staticCondMatrixManager[mkey];
1946  }
1947 
1948  void HexExp::v_DropLocStaticCondMatrix(const MatrixKey &mkey)
1949  {
1950  m_staticCondMatrixManager.DeleteObject(mkey);
1951  }
1952 
1953  void HexExp::v_LaplacianMatrixOp_MatFree_Kernel(
1954  const Array<OneD, const NekDouble> &inarray,
1955  Array<OneD, NekDouble> &outarray,
1956  Array<OneD, NekDouble> &wsp)
1957  {
1958  // This implementation is only valid when there are no
1959  // coefficients associated to the Laplacian operator
1960  if (m_metrics.count(eMetricLaplacian00) == 0)
1961  {
1962  ComputeLaplacianMetric();
1963  }
1964 
1965  int nquad0 = m_base[0]->GetNumPoints();
1966  int nquad1 = m_base[1]->GetNumPoints();
1967  int nquad2 = m_base[2]->GetNumPoints();
1968  int nqtot = nquad0*nquad1*nquad2;
1969 
1970  ASSERTL1(wsp.size() >= 6*nqtot,
1971  "Insufficient workspace size.");
1972 
1973  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
1974  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
1975  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata();
1976  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
1977  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
1978  const Array<OneD, const NekDouble>& dbase2 = m_base[2]->GetDbdata();
1979  const Array<OneD, const NekDouble>& metric00 = m_metrics[eMetricLaplacian00];
1980  const Array<OneD, const NekDouble>& metric01 = m_metrics[eMetricLaplacian01];
1981  const Array<OneD, const NekDouble>& metric02 = m_metrics[eMetricLaplacian02];
1982  const Array<OneD, const NekDouble>& metric11 = m_metrics[eMetricLaplacian11];
1983  const Array<OneD, const NekDouble>& metric12 = m_metrics[eMetricLaplacian12];
1984  const Array<OneD, const NekDouble>& metric22 = m_metrics[eMetricLaplacian22];
1985 
1986  // Allocate temporary storage
1987  Array<OneD,NekDouble> wsp0(wsp);
1988  Array<OneD,NekDouble> wsp1(wsp+1*nqtot);
1989  Array<OneD,NekDouble> wsp2(wsp+2*nqtot);
1990  Array<OneD,NekDouble> wsp3(wsp+3*nqtot);
1991  Array<OneD,NekDouble> wsp4(wsp+4*nqtot);
1992  Array<OneD,NekDouble> wsp5(wsp+5*nqtot);
1993 
1994  StdExpansion3D::PhysTensorDeriv(inarray,wsp0,wsp1,wsp2);
1995 
1996  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1997  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1998  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1999  // especially for this purpose
2000  Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp0[0],1,&metric01[0],1,&wsp1[0],1,&wsp3[0],1);
2001  Vmath::Vvtvp (nqtot,&metric02[0],1,&wsp2[0],1,&wsp3[0],1,&wsp3[0],1);
2002  Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp0[0],1,&metric11[0],1,&wsp1[0],1,&wsp4[0],1);
2003  Vmath::Vvtvp (nqtot,&metric12[0],1,&wsp2[0],1,&wsp4[0],1,&wsp4[0],1);
2004  Vmath::Vvtvvtp(nqtot,&metric02[0],1,&wsp0[0],1,&metric12[0],1,&wsp1[0],1,&wsp5[0],1);
2005  Vmath::Vvtvp (nqtot,&metric22[0],1,&wsp2[0],1,&wsp5[0],1,&wsp5[0],1);
2006 
2007  // outarray = m = (D_xi1 * B)^T * k
2008  // wsp1 = n = (D_xi2 * B)^T * l
2009  IProductWRTBase_SumFacKernel(dbase0,base1,base2,wsp3,outarray,wsp0,false,true,true);
2010  IProductWRTBase_SumFacKernel(base0,dbase1,base2,wsp4,wsp2, wsp0,true,false,true);
2011  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
2012  IProductWRTBase_SumFacKernel(base0,base1,dbase2,wsp5,wsp2, wsp0,true,true,false);
2013  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
2014  }
2015 
2016 
2017  void HexExp::v_ComputeLaplacianMetric()
2018  {
2019  if (m_metrics.count(eMetricQuadrature) == 0)
2020  {
2021  ComputeQuadratureMetric();
2022  }
2023 
2024  const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
2025  const unsigned int nqtot = GetTotPoints();
2026  const unsigned int dim = 3;
2027  const MetricType m[3][3] = { {eMetricLaplacian00, eMetricLaplacian01, eMetricLaplacian02},
2028  {eMetricLaplacian01, eMetricLaplacian11, eMetricLaplacian12},
2029  {eMetricLaplacian02, eMetricLaplacian12, eMetricLaplacian22}
2030  };
2031 
2032  for (unsigned int i = 0; i < dim; ++i)
2033  {
2034  for (unsigned int j = i; j < dim; ++j)
2035  {
2036  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
2037  const Array<TwoD, const NekDouble> &gmat =
2038  m_metricinfo->GetGmat(GetPointsKeys());
2039  if (type == SpatialDomains::eDeformed)
2040  {
2041  Vmath::Vcopy(nqtot, &gmat[i*dim+j][0], 1,
2042  &m_metrics[m[i][j]][0], 1);
2043  }
2044  else
2045  {
2046  Vmath::Fill(nqtot, gmat[i*dim+j][0],
2047  &m_metrics[m[i][j]][0], 1);
2048  }
2049  MultiplyByQuadratureMetric(m_metrics[m[i][j]],
2050  m_metrics[m[i][j]]);
2051 
2052  }
2053  }
2054  }
2055 
2056  }//end of namespace
2057 }//end of namespace
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:274
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::LibUtilities::BasisKey::GetNumPoints
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:133
Nektar::LibUtilities::BasisKey::GetPointsKey
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:150
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:60
Nektar::LocalRegions::Expansion3D::m_faceNormals
std::map< int, NormalVector > m_faceNormals
Definition: Expansion3D.h:123
Nektar::LocalRegions::Expansion3D::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: Expansion3D.cpp:451
Nektar::LocalRegions::Expansion3D::GetGeom3D
SpatialDomains::Geometry3DSharedPtr GetGeom3D() const
Definition: Expansion3D.h:191
Nektar::LocalRegions::Expansion::m_elmtBndNormDirElmtLen
std::map< int, Array< OneD, NekDouble > > m_elmtBndNormDirElmtLen
the element length in each element boundary(Vertex, edge or face) normal direction calculated based o...
Definition: Expansion.h:284
Nektar::LocalRegions::Expansion::GetGeom
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:172
Nektar::LocalRegions::Expansion::m_geom
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:272
Nektar::LocalRegions::Expansion::BuildVertexMatrix
DNekMatSharedPtr BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
Definition: Expansion.cpp:103
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:399
Nektar::LocalRegions::Expansion::m_metricinfo
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:273
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:318
Nektar::LocalRegions::Expansion::GetLocMatrix
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:90
Nektar::LocalRegions::Expansion::BuildTransformationMatrix
DNekMatSharedPtr BuildTransformationMatrix(const DNekScalMatSharedPtr &r_bnd, const StdRegions::MatrixType matrixType)
Definition: Expansion.cpp:95
Nektar::LocalRegions::HexExp::IProductWRTDerivBase_MatOp
void IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: HexExp.cpp:557
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::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:678
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:1953
Nektar::LocalRegions::HexExp::v_StdPhysEvaluate
virtual NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
Definition: HexExp.cpp:670
Nektar::LocalRegions::HexExp::m_matrixManager
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: HexExp.h:288
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:756
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:1270
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:1262
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:1298
Nektar::LocalRegions::HexExp::m_staticCondMatrixManager
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: HexExp.h:289
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_GetTracePhysMap
virtual void v_GetTracePhysMap(const int face, Array< OneD, int > &outarray)
Definition: HexExp.cpp:877
Nektar::LocalRegions::HexExp::v_DetShapeType
virtual LibUtilities::ShapeType v_DetShapeType() const
Return the region shape using the enum-list of ShapeType.
Definition: HexExp.cpp:750
Nektar::LocalRegions::HexExp::v_DropLocStaticCondMatrix
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: HexExp.cpp:1948
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::HexExp::v_AlignVectorToCollapsedDir
virtual void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: HexExp.cpp:519
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:603
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::LocalRegions::HexExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1453
Nektar::LocalRegions::HexExp::CreateStaticCondMatrix
DNekScalBlkMatSharedPtr CreateStaticCondMatrix(const MatrixKey &mkey)
Definition: HexExp.cpp:1822
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
Nektar::LocalRegions::HexExp::v_GetLocMatrix
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: HexExp.cpp:1936
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:1289
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:1316
Nektar::LocalRegions::HexExp::v_ComputeTraceNormal
void v_ComputeTraceNormal(const int face)
Definition: HexExp.cpp:993
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:737
Nektar::LocalRegions::HexExp::v_GetLinStdExp
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
Definition: HexExp.cpp:698
Nektar::LocalRegions::HexExp::v_ComputeLaplacianMetric
virtual void v_ComputeLaplacianMetric()
Definition: HexExp.cpp:2017
Nektar::LocalRegions::HexExp::v_GetStdExp
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
Definition: HexExp.cpp:689
Nektar::LocalRegions::HexExp::v_GetLocStaticCondMatrix
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: HexExp.cpp:1942
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:1325
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:1307
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::HexExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1485
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_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1509
Nektar::LocalRegions::HexExp::v_ReduceOrderCoeffs
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: HexExp.cpp:1384
Nektar::LocalRegions::HexExp::CreateMatrix
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition: HexExp.cpp:1523
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_GetFaceDGForwards
virtual bool v_GetFaceDGForwards(const int i) const
Definition: HexExp.cpp:867
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:718
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:84
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:171
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:305
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::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:347
Nektar::StdRegions::StdExpansion::GetBoundaryMap
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:687
Nektar::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:124
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:134
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:158
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:617
Nektar::StdRegions::StdExpansion::GetPointsKeys
const LibUtilities::PointsKeyVector GetPointsKeys() const
Definition: StdExpansion.h:1078
Nektar::StdRegions::StdExpansion::GetStdStaticCondMatrix
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:622
Nektar::StdRegions::StdExpansion::GetLocStaticCondMatrix
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
Definition: StdExpansion.h:660
Nektar::StdRegions::StdExpansion::GetTraceBasisKey
const LibUtilities::BasisKey GetTraceBasisKey(const int i, int k=-1) const
This function returns the basis key belonging to the i-th trace.
Definition: StdExpansion.h:304
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:537
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:376
Nektar::StdRegions::StdExpansion::GetInteriorMap
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:692
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:433
Nektar::StdRegions::StdExpansion::GenMatrix
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:850
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:208
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:1745
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:221
Nektar::StdRegions::StdExpansion::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:733
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:855
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1162
Nektar::StdRegions::StdHexExp
Class representing a hexehedral element in reference space.
Definition: StdHexExp.h:48
Nektar::StdRegions::StdMatrixKey::GetConstFactors
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:135
Nektar::StdRegions::StdMatrixKey::GetShapeType
LibUtilities::ShapeType GetShapeType() const
Definition: StdMatrixKey.h:86
Nektar::StdRegions::StdMatrixKey::GetVarCoeffs
const VarCoeffMap & GetVarCoeffs() const
Definition: StdMatrixKey.h:161
Nektar::StdRegions::StdMatrixKey::GetMatrixType
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:81
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:121
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:130
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:265
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:167
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::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::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::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
Nektar::LibUtilities::eGaussLobattoLegendre
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
Nektar::LibUtilities::eOrtho_A
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:45
Nektar::LibUtilities::eGLL_Lagrange
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:54
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
Nektar::SpatialDomains::GeomFactorsSharedPtr
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
Nektar::SpatialDomains::GeomType
GeomType
Indicates the type of element geometry.
Definition: SpatialDomains.hpp:52
Nektar::SpatialDomains::eRegular
@ eRegular
Geometry is straight-sided with constant geometric factors.
Definition: SpatialDomains.hpp:54
Nektar::SpatialDomains::eNoGeomType
@ eNoGeomType
No type defined.
Definition: SpatialDomains.hpp:53
Nektar::SpatialDomains::eMovingRegular
@ eMovingRegular
Currently unused.
Definition: SpatialDomains.hpp:57
Nektar::SpatialDomains::eDeformed
@ eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:56
Nektar::SpatialDomains::HexGeomSharedPtr
std::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:84
Nektar::StdRegions::StdExpansionSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1714
Nektar::StdRegions::StdHexExpSharedPtr
std::shared_ptr< StdHexExp > StdHexExpSharedPtr
Definition: StdHexExp.h:293
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::DNekScalMatSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:69
Nektar::DNekBlkMatSharedPtr
std::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:71
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:217
Nektar::DNekScalBlkMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:73
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:820
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:475
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:192
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:513
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:322
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
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:291
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:257
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:436
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
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:625
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267
Generated on Wed Nov 24 2021 10:33:27 for Nektar++ by doxygen 1.9.1