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