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