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