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