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.size())
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.size())
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.size())
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.size())
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.size())
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.size())
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 
379  Array<OneD, NekDouble> tmp1 (nqtot);
380  Array<OneD, NekDouble> tmp2 (nqtot);
381  Array<OneD, NekDouble> tmp3 (nqtot);
382  Array<OneD, NekDouble> tmp4 (nqtot);
383  Array<OneD, NekDouble> tmp6 (m_ncoeffs);
384  Array<OneD, NekDouble> wsp (nquad1*nquad2*order0 +
385  nquad2*order0*(order1+1)/2);
386 
387  MultiplyByQuadratureMetric(inarray,tmp1);
388 
389  Array<OneD, Array<OneD, NekDouble>> tmp2D{3};
390  tmp2D[0] = tmp2;
391  tmp2D[1] = tmp3;
392  tmp2D[2] = tmp4;
393 
394  TetExp::v_AlignVectorToCollapsedDir(dir, tmp1, tmp2D);
395 
396  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
397  m_base[1]->GetBdata (),
398  m_base[2]->GetBdata (),
399  tmp2,outarray,wsp,
400  false,true,true);
401 
402  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata (),
403  m_base[1]->GetDbdata(),
404  m_base[2]->GetBdata (),
405  tmp3,tmp6,wsp,
406  true,false,true);
407 
408  Vmath::Vadd(m_ncoeffs, tmp6, 1, outarray, 1, outarray, 1);
409 
410  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata (),
411  m_base[1]->GetBdata (),
412  m_base[2]->GetDbdata(),
413  tmp4,tmp6,wsp,
414  true,true,false);
415 
416  Vmath::Vadd(m_ncoeffs, tmp6, 1, outarray, 1, outarray, 1);
417  }
418 
419  void TetExp::v_AlignVectorToCollapsedDir(
420  const int dir,
421  const Array<OneD, const NekDouble> &inarray,
422  Array<OneD, Array<OneD, NekDouble> > &outarray)
423  {
424  int i, j;
425 
426  const int nquad0 = m_base[0]->GetNumPoints();
427  const int nquad1 = m_base[1]->GetNumPoints();
428  const int nquad2 = m_base[2]->GetNumPoints();
429  const int nqtot = nquad0*nquad1*nquad2;
430 
431  const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
432  const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
433  const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
434 
435  Array<OneD, NekDouble> tmp2 (nqtot);
436  Array<OneD, NekDouble> tmp3 (nqtot);
437 
438  const Array<TwoD, const NekDouble>& df =
439  m_metricinfo->GetDerivFactors(GetPointsKeys());
440 
441  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
442  {
443  Vmath::Vmul(nqtot,&df[3*dir][0], 1,inarray.get(),1,
444  tmp2.get(),1);
445  Vmath::Vmul(nqtot,&df[3*dir+1][0],1,inarray.get(),1,
446  tmp3.get(),1);
447  Vmath::Vmul(nqtot,&df[3*dir+2][0],1,inarray.get(),1,
448  outarray[2].get(),1);
449  }
450  else
451  {
452  Vmath::Smul(nqtot, df[3*dir ][0],inarray.get(),1,
453  tmp2.get(), 1);
454  Vmath::Smul(nqtot, df[3*dir+1][0],inarray.get(),1,
455  tmp3.get(), 1);
456  Vmath::Smul(nqtot, df[3*dir+2][0],inarray.get(),1,
457  outarray[2].get(), 1);
458  }
459 
460  NekDouble g0,g1,g1a,g2,g3;
461  int k,cnt;
462 
463  for(cnt = 0, k = 0; k < nquad2; ++k)
464  {
465  for(j = 0; j < nquad1; ++j)
466  {
467  g2 = 2.0/(1.0-z2[k]);
468  g1 = g2/(1.0-z1[j]);
469  g0 = 2.0*g1;
470  g3 = (1.0+z1[j])*g2*0.5;
471 
472  for(i = 0; i < nquad0; ++i,++cnt)
473  {
474  g1a = g1*(1+z0[i]);
475 
476  outarray[0][cnt] = g0*tmp2[cnt] + g1a*(tmp3[cnt] +
477  outarray[2][cnt]);
478 
479  outarray[1][cnt] = g2*tmp3[cnt] + g3*outarray[2][cnt];
480 
481  }
482  }
483  }
484  }
485 
486 
487  //-----------------------------
488  // Evaluation functions
489  //-----------------------------
490 
491  /**
492  * Given the local cartesian coordinate \a Lcoord evaluate the
493  * value of physvals at this point by calling through to the
494  * StdExpansion method
495  */
496  NekDouble TetExp::v_StdPhysEvaluate(
497  const Array<OneD, const NekDouble> &Lcoord,
498  const Array<OneD, const NekDouble> &physvals)
499  {
500  // Evaluate point in local (eta) coordinates.
501  return StdTetExp::v_PhysEvaluate(Lcoord,physvals);
502  }
503 
504  /**
505  * @param coord Physical space coordinate
506  * @returns Evaluation of expansion at given coordinate.
507  */
508  NekDouble TetExp::v_PhysEvaluate(
509  const Array<OneD, const NekDouble> &coord,
510  const Array<OneD, const NekDouble> & physvals)
511  {
512  ASSERTL0(m_geom,"m_geom not defined");
513 
514  Array<OneD,NekDouble> Lcoord = Array<OneD,NekDouble>(3);
515 
516  // Get the local (eta) coordinates of the point
517  m_geom->GetLocCoords(coord,Lcoord);
518 
519  // Evaluate point in local (eta) coordinates.
520  return StdTetExp::v_PhysEvaluate(Lcoord,physvals);
521  }
522 
523  /**
524  * \brief Get the coordinates "coords" at the local coordinates "Lcoords"
525  */
526  void TetExp::v_GetCoord(
527  const Array<OneD, const NekDouble> &Lcoords,
528  Array<OneD,NekDouble> &coords)
529  {
530  int i;
531 
532  ASSERTL1(Lcoords[0] <= -1.0 && Lcoords[0] >= 1.0 &&
533  Lcoords[1] <= -1.0 && Lcoords[1] >= 1.0 &&
534  Lcoords[2] <= -1.0 && Lcoords[2] >= 1.0,
535  "Local coordinates are not in region [-1,1]");
536 
537  // m_geom->FillGeom(); // TODO: implement FillGeom()
538 
539  for(i = 0; i < m_geom->GetCoordim(); ++i)
540  {
541  coords[i] = m_geom->GetCoord(i,Lcoords);
542  }
543  }
544 
545  void TetExp::v_GetCoords(
546  Array<OneD, NekDouble> &coords_0,
547  Array<OneD, NekDouble> &coords_1,
548  Array<OneD, NekDouble> &coords_2)
549  {
550  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
551  }
552 
553 
554  //-----------------------------
555  // Helper functions
556  //-----------------------------
557 
558  /**
559  * \brief Return Shape of region, using ShapeType enum list.
560  */
561  LibUtilities::ShapeType TetExp::v_DetShapeType() const
562  {
563  return LibUtilities::eTetrahedron;
564  }
565 
566  StdRegions::StdExpansionSharedPtr TetExp::v_GetStdExp(void) const
567  {
568  return MemoryManager<StdRegions::StdTetExp>
569  ::AllocateSharedPtr(m_base[0]->GetBasisKey(),
570  m_base[1]->GetBasisKey(),
571  m_base[2]->GetBasisKey());
572  }
573 
574 
575  StdRegions::StdExpansionSharedPtr TetExp::v_GetLinStdExp(void) const
576  {
577  LibUtilities::BasisKey bkey0(m_base[0]->GetBasisType(),
578  2, m_base[0]->GetPointsKey());
579  LibUtilities::BasisKey bkey1(m_base[1]->GetBasisType(),
580  2, m_base[1]->GetPointsKey());
581  LibUtilities::BasisKey bkey2(m_base[2]->GetBasisType(),
582  2, m_base[2]->GetPointsKey());
583 
584  return MemoryManager<StdRegions::StdTetExp>
585  ::AllocateSharedPtr( bkey0, bkey1, bkey2);
586  }
587 
588 
589  int TetExp::v_GetCoordim()
590  {
591  return m_geom->GetCoordim();
592  }
593 
594  void TetExp::v_ExtractDataToCoeffs(
595  const NekDouble *data,
596  const std::vector<unsigned int > &nummodes,
597  const int mode_offset,
598  NekDouble * coeffs,
599  std::vector<LibUtilities::BasisType> &fromType)
600  {
601  boost::ignore_unused(fromType);
602 
603  int data_order0 = nummodes[mode_offset];
604  int fillorder0 = min(m_base[0]->GetNumModes(),data_order0);
605  int data_order1 = nummodes[mode_offset+1];
606  int order1 = m_base[1]->GetNumModes();
607  int fillorder1 = min(order1,data_order1);
608  int data_order2 = nummodes[mode_offset+2];
609  int order2 = m_base[2]->GetNumModes();
610  int fillorder2 = min(order2,data_order2);
611 
612  switch(m_base[0]->GetBasisType())
613  {
614  case LibUtilities::eModified_A:
615  {
616  int i,j;
617  int cnt = 0;
618  int cnt1 = 0;
619 
620  ASSERTL1(m_base[1]->GetBasisType() ==
621  LibUtilities::eModified_B,
622  "Extraction routine not set up for this basis");
623  ASSERTL1(m_base[2]->GetBasisType() ==
624  LibUtilities::eModified_C,
625  "Extraction routine not set up for this basis");
626 
627  Vmath::Zero(m_ncoeffs,coeffs,1);
628  for(j = 0; j < fillorder0; ++j)
629  {
630  for(i = 0; i < fillorder1-j; ++i)
631  {
632  Vmath::Vcopy(fillorder2-j-i, &data[cnt], 1,
633  &coeffs[cnt1], 1);
634  cnt += data_order2-j-i;
635  cnt1 += order2-j-i;
636  }
637 
638  // count out data for j iteration
639  for(i = fillorder1-j; i < data_order1-j; ++i)
640  {
641  cnt += data_order2-j-i;
642  }
643 
644  for(i = fillorder1-j; i < order1-j; ++i)
645  {
646  cnt1 += order2-j-i;
647  }
648 
649  }
650  }
651  break;
652  default:
653  ASSERTL0(false, "basis is either not set up or not "
654  "hierarchicial");
655  }
656  }
657 
658  /**
659  * \brief Returns the physical values at the quadrature points of a face
660  */
661  void TetExp::v_GetTracePhysMap(const int face,
662  Array<OneD, int> &outarray)
663  {
664  int nquad0 = m_base[0]->GetNumPoints();
665  int nquad1 = m_base[1]->GetNumPoints();
666  int nquad2 = m_base[2]->GetNumPoints();
667 
668  int nq0 = 0;
669  int nq1 = 0;
670 
671  // get forward aligned faces.
672  switch(face)
673  {
674  case 0:
675  {
676  nq0 = nquad0;
677  nq1 = nquad1;
678  if(outarray.size()!=nq0*nq1)
679  {
680  outarray = Array<OneD, int>(nq0*nq1);
681  }
682 
683  for (int i = 0; i < nquad0*nquad1; ++i)
684  {
685  outarray[i] = i;
686  }
687 
688  break;
689  }
690  case 1:
691  {
692  nq0 = nquad0;
693  nq1 = nquad2;
694  if(outarray.size()!=nq0*nq1)
695  {
696  outarray = Array<OneD, int>(nq0*nq1);
697  }
698 
699  //Direction A and B positive
700  for (int k=0; k<nquad2; k++)
701  {
702  for(int i = 0; i < nquad0; ++i)
703  {
704  outarray[k*nquad0+i] = (nquad0*nquad1*k)+i;
705  }
706  }
707  break;
708  }
709  case 2:
710  {
711  nq0 = nquad1;
712  nq1 = nquad2;
713  if(outarray.size()!=nq0*nq1)
714  {
715  outarray = Array<OneD, int>(nq0*nq1);
716  }
717 
718  //Directions A and B positive
719  for(int j = 0; j < nquad1*nquad2; ++j)
720  {
721  outarray[j] = nquad0-1 + j*nquad0;
722  }
723  break;
724  }
725  case 3:
726  {
727  nq0 = nquad1;
728  nq1 = nquad2;
729  if(outarray.size() != nq0*nq1)
730  {
731  outarray = Array<OneD, int>(nq0*nq1);
732  }
733 
734  //Directions A and B positive
735  for(int j = 0; j < nquad1*nquad2; ++j)
736  {
737  outarray[j] = j*nquad0;
738  }
739  }
740  break;
741  default:
742  ASSERTL0(false,"face value (> 3) is out of range");
743  break;
744  }
745  }
746 
747 
748  /**
749  * \brief Compute the normal of a triangular face
750  */
751  void TetExp::v_ComputeTraceNormal(const int face)
752  {
753  int i;
754  const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
755  GetGeom()->GetMetricInfo();
756 
757  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
758  for(int i = 0; i < ptsKeys.size(); ++i)
759  {
760  // Need at least 2 points for computing normals
761  if (ptsKeys[i].GetNumPoints() == 1)
762  {
763  LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
764  ptsKeys[i] = pKey;
765  }
766  }
767 
768  SpatialDomains::GeomType type = geomFactors->GetGtype();
769  const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys);
770  const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
771 
772  LibUtilities::BasisKey tobasis0 = GetTraceBasisKey(face,0);
773  LibUtilities::BasisKey tobasis1 = GetTraceBasisKey(face,1);
774 
775  // number of face quadrature points
776  int nq_face = tobasis0.GetNumPoints()*tobasis1.GetNumPoints();
777 
778  int vCoordDim = GetCoordim();
779 
780  m_faceNormals[face] = Array<OneD, Array<OneD, NekDouble> >(vCoordDim);
781  Array<OneD, Array<OneD, NekDouble> > &normal = m_faceNormals[face];
782  for (i = 0; i < vCoordDim; ++i)
783  {
784  normal[i] = Array<OneD, NekDouble>(nq_face);
785  }
786 
787  size_t nqb = nq_face;
788  size_t nbnd= face;
789  m_elmtBndNormDirElmtLen[nbnd] = Array<OneD, NekDouble> {nqb, 0.0};
790  Array<OneD, NekDouble> &length = m_elmtBndNormDirElmtLen[nbnd];
791 
792  // Regular geometry case
793  if (type == SpatialDomains::eRegular ||
794  type == SpatialDomains::eMovingRegular)
795  {
796  NekDouble fac;
797 
798  // Set up normals
799  switch (face)
800  {
801  case 0:
802  {
803  for (i = 0; i < vCoordDim; ++i)
804  {
805  normal[i][0] = -df[3*i+2][0];
806  }
807 
808  break;
809  }
810  case 1:
811  {
812  for (i = 0; i < vCoordDim; ++i)
813  {
814  normal[i][0] = -df[3*i+1][0];
815  }
816 
817  break;
818  }
819  case 2:
820  {
821  for (i = 0; i < vCoordDim; ++i)
822  {
823  normal[i][0] = df[3*i][0]+df[3*i+1][0]+
824  df[3*i+2][0];
825  }
826 
827  break;
828  }
829  case 3:
830  {
831  for(i = 0; i < vCoordDim; ++i)
832  {
833  normal[i][0] = -df[3*i][0];
834  }
835  break;
836  }
837  default:
838  ASSERTL0(false,"face is out of range (edge < 3)");
839  }
840 
841  // normalise
842  fac = 0.0;
843  for (i = 0; i < vCoordDim; ++i)
844  {
845  fac += normal[i][0]*normal[i][0];
846  }
847  fac = 1.0/sqrt(fac);
848  Vmath::Fill(nqb, fac, length, 1);
849 
850  for (i = 0; i < vCoordDim; ++i)
851  {
852  Vmath::Fill(nq_face,fac*normal[i][0],normal[i],1);
853  }
854  }
855  else
856  {
857  // Set up deformed normals
858  int j, k;
859 
860  int nq0 = ptsKeys[0].GetNumPoints();
861  int nq1 = ptsKeys[1].GetNumPoints();
862  int nq2 = ptsKeys[2].GetNumPoints();
863  int nqtot;
864  int nq01 =nq0*nq1;
865 
866  // number of elemental quad points
867  if (face == 0)
868  {
869  nqtot = nq01;
870  }
871  else if (face == 1)
872  {
873  nqtot = nq0*nq2;
874  }
875  else
876  {
877  nqtot = nq1*nq2;
878  }
879 
880  LibUtilities::PointsKey points0;
881  LibUtilities::PointsKey points1;
882 
883  Array<OneD, NekDouble> faceJac(nqtot);
884  Array<OneD,NekDouble> normals(vCoordDim*nqtot, 0.0);
885 
886  // Extract Jacobian along face and recover local derivates
887  // (dx/dr) for polynomial interpolation by multiplying m_gmat by
888  // jacobian
889  switch (face)
890  {
891  case 0:
892  {
893  for(j = 0; j < nq01; ++j)
894  {
895  normals[j] = -df[2][j]*jac[j];
896  normals[nqtot+j] = -df[5][j]*jac[j];
897  normals[2*nqtot+j] = -df[8][j]*jac[j];
898  faceJac[j] = jac[j];
899  }
900 
901  points0 = ptsKeys[0];
902  points1 = ptsKeys[1];
903  break;
904  }
905 
906  case 1:
907  {
908  for (j = 0; j < nq0; ++j)
909  {
910  for(k = 0; k < nq2; ++k)
911  {
912  int tmp = j+nq01*k;
913  normals[j+k*nq0] =
914  -df[1][tmp]*jac[tmp];
915  normals[nqtot+j+k*nq0] =
916  -df[4][tmp]*jac[tmp];
917  normals[2*nqtot+j+k*nq0] =
918  -df[7][tmp]*jac[tmp];
919  faceJac[j+k*nq0] = jac[tmp];
920  }
921  }
922 
923  points0 = ptsKeys[0];
924  points1 = ptsKeys[2];
925  break;
926  }
927 
928  case 2:
929  {
930  for (j = 0; j < nq1; ++j)
931  {
932  for(k = 0; k < nq2; ++k)
933  {
934  int tmp = nq0-1+nq0*j+nq01*k;
935  normals[j+k*nq1] =
936  (df[0][tmp]+df[1][tmp]+df[2][tmp])*
937  jac[tmp];
938  normals[nqtot+j+k*nq1] =
939  (df[3][tmp]+df[4][tmp]+df[5][tmp])*
940  jac[tmp];
941  normals[2*nqtot+j+k*nq1] =
942  (df[6][tmp]+df[7][tmp]+df[8][tmp])*
943  jac[tmp];
944  faceJac[j+k*nq1] = jac[tmp];
945  }
946  }
947 
948  points0 = ptsKeys[1];
949  points1 = ptsKeys[2];
950  break;
951  }
952 
953  case 3:
954  {
955  for (j = 0; j < nq1; ++j)
956  {
957  for(k = 0; k < nq2; ++k)
958  {
959  int tmp = j*nq0+nq01*k;
960  normals[j+k*nq1] =
961  -df[0][tmp]*jac[tmp];
962  normals[nqtot+j+k*nq1] =
963  -df[3][tmp]*jac[tmp];
964  normals[2*nqtot+j+k*nq1] =
965  -df[6][tmp]*jac[tmp];
966  faceJac[j+k*nq1] = jac[tmp];
967  }
968  }
969 
970  points0 = ptsKeys[1];
971  points1 = ptsKeys[2];
972  break;
973  }
974 
975  default:
976  ASSERTL0(false,"face is out of range (face < 3)");
977  }
978 
979  Array<OneD,NekDouble> work (nq_face, 0.0);
980  // Interpolate Jacobian and invert
981  LibUtilities::Interp2D(points0, points1, faceJac,
982  tobasis0.GetPointsKey(),
983  tobasis1.GetPointsKey(),
984  work);
985  Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
986 
987  // Interpolate normal and multiply by inverse Jacobian.
988  for(i = 0; i < vCoordDim; ++i)
989  {
990  LibUtilities::Interp2D(points0, points1,
991  &normals[i*nqtot],
992  tobasis0.GetPointsKey(),
993  tobasis1.GetPointsKey(),
994  &normal[i][0]);
995  Vmath::Vmul(nq_face,work,1,normal[i],1,normal[i],1);
996  }
997 
998  // Normalise to obtain unit normals.
999  Vmath::Zero(nq_face,work,1);
1000  for(i = 0; i < GetCoordim(); ++i)
1001  {
1002  Vmath::Vvtvp(nq_face,normal[i],1,normal[i],1,work,1,work,1);
1003  }
1004 
1005  Vmath::Vsqrt(nq_face,work,1,work,1);
1006  Vmath::Sdiv (nq_face,1.0,work,1,work,1);
1007 
1008  Vmath::Vcopy(nqb, work, 1, length, 1);
1009 
1010  for(i = 0; i < GetCoordim(); ++i)
1011  {
1012  Vmath::Vmul(nq_face,normal[i],1,work,1,normal[i],1);
1013  }
1014  }
1015  }
1016 
1017  //-----------------------------
1018  // Operator creation functions
1019  //-----------------------------
1020  void TetExp::v_HelmholtzMatrixOp(
1021  const Array<OneD, const NekDouble> &inarray,
1022  Array<OneD,NekDouble> &outarray,
1023  const StdRegions::StdMatrixKey &mkey)
1024  {
1025  TetExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1026  }
1027 
1028 
1029  void TetExp::v_LaplacianMatrixOp(
1030  const Array<OneD, const NekDouble> &inarray,
1031  Array<OneD,NekDouble> &outarray,
1032  const StdRegions::StdMatrixKey &mkey)
1033  {
1034  TetExp::v_LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1035  }
1036 
1037  void TetExp::v_LaplacianMatrixOp(
1038  const int k1,
1039  const int k2,
1040  const Array<OneD, const NekDouble> &inarray,
1041  Array<OneD,NekDouble> &outarray,
1042  const StdRegions::StdMatrixKey &mkey)
1043  {
1044  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,
1045  mkey);
1046  }
1047 
1048  void TetExp::v_SVVLaplacianFilter(
1049  Array<OneD, NekDouble> &array,
1050  const StdRegions::StdMatrixKey &mkey)
1051  {
1052  int nq = GetTotPoints();
1053 
1054  // Calculate sqrt of the Jacobian
1055  Array<OneD, const NekDouble> jac =
1056  m_metricinfo->GetJac(GetPointsKeys());
1057  Array<OneD, NekDouble> sqrt_jac(nq);
1058  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1059  {
1060  Vmath::Vsqrt(nq,jac,1,sqrt_jac,1);
1061  }
1062  else
1063  {
1064  Vmath::Fill(nq,sqrt(jac[0]),sqrt_jac,1);
1065  }
1066 
1067  // Multiply array by sqrt(Jac)
1068  Vmath::Vmul(nq,sqrt_jac,1,array,1,array,1);
1069 
1070  // Apply std region filter
1071  StdTetExp::v_SVVLaplacianFilter( array, mkey);
1072 
1073  // Divide by sqrt(Jac)
1074  Vmath::Vdiv(nq,array,1,sqrt_jac,1,array,1);
1075  }
1076 
1077 
1078  //-----------------------------
1079  // Matrix creation functions
1080  //-----------------------------
1081  DNekMatSharedPtr TetExp::v_GenMatrix(
1082  const StdRegions::StdMatrixKey &mkey)
1083  {
1084  DNekMatSharedPtr returnval;
1085 
1086  switch(mkey.GetMatrixType())
1087  {
1088  case StdRegions::eHybridDGHelmholtz:
1089  case StdRegions::eHybridDGLamToU:
1090  case StdRegions::eHybridDGLamToQ0:
1091  case StdRegions::eHybridDGLamToQ1:
1092  case StdRegions::eHybridDGLamToQ2:
1093  case StdRegions::eHybridDGHelmBndLam:
1094  case StdRegions::eInvLaplacianWithUnityMean:
1095  returnval = Expansion3D::v_GenMatrix(mkey);
1096  break;
1097  default:
1098  returnval = StdTetExp::v_GenMatrix(mkey);
1099  }
1100 
1101  return returnval;
1102  }
1103 
1104 
1105  DNekScalMatSharedPtr TetExp::CreateMatrix(const MatrixKey &mkey)
1106  {
1107  DNekScalMatSharedPtr returnval;
1108  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1109 
1110  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
1111 
1112  switch(mkey.GetMatrixType())
1113  {
1114  case StdRegions::eMass:
1115  {
1116  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1117  mkey.GetNVarCoeff())
1118  {
1119  NekDouble one = 1.0;
1120  DNekMatSharedPtr mat = GenMatrix(mkey);
1121  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1122  }
1123  else
1124  {
1125  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1126  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1127  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1128  }
1129  }
1130  break;
1131  case StdRegions::eInvMass:
1132  {
1133  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1134  {
1135  NekDouble one = 1.0;
1136  StdRegions::StdMatrixKey masskey(StdRegions::eMass,DetShapeType(),
1137  *this);
1138  DNekMatSharedPtr mat = GenMatrix(masskey);
1139  mat->Invert();
1140  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1141  }
1142  else
1143  {
1144  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1145  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1146  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(fac,mat);
1147  }
1148  }
1149  break;
1150  case StdRegions::eWeakDeriv0:
1151  case StdRegions::eWeakDeriv1:
1152  case StdRegions::eWeakDeriv2:
1153  {
1154  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1155  mkey.GetNVarCoeff())
1156  {
1157  NekDouble one = 1.0;
1158  DNekMatSharedPtr mat = GenMatrix(mkey);
1159 
1160  returnval = MemoryManager<DNekScalMat>
1161  ::AllocateSharedPtr(one,mat);
1162  }
1163  else
1164  {
1165  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1166  Array<TwoD, const NekDouble> df
1167  = m_metricinfo->GetDerivFactors(ptsKeys);
1168  int dir = 0;
1169 
1170  switch(mkey.GetMatrixType())
1171  {
1172  case StdRegions::eWeakDeriv0:
1173  dir = 0;
1174  break;
1175  case StdRegions::eWeakDeriv1:
1176  dir = 1;
1177  break;
1178  case StdRegions::eWeakDeriv2:
1179  dir = 2;
1180  break;
1181  default:
1182  break;
1183  }
1184 
1185  MatrixKey deriv0key(StdRegions::eWeakDeriv0,
1186  mkey.GetShapeType(), *this);
1187  MatrixKey deriv1key(StdRegions::eWeakDeriv1,
1188  mkey.GetShapeType(), *this);
1189  MatrixKey deriv2key(StdRegions::eWeakDeriv2,
1190  mkey.GetShapeType(), *this);
1191 
1192  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
1193  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
1194  DNekMat &deriv2 = *GetStdMatrix(deriv2key);
1195 
1196  int rows = deriv0.GetRows();
1197  int cols = deriv1.GetColumns();
1198 
1199  DNekMatSharedPtr WeakDeriv = MemoryManager<DNekMat>
1200  ::AllocateSharedPtr(rows,cols);
1201  (*WeakDeriv) = df[3*dir][0]*deriv0
1202  + df[3*dir+1][0]*deriv1
1203  + df[3*dir+2][0]*deriv2;
1204 
1205  returnval = MemoryManager<DNekScalMat>
1206  ::AllocateSharedPtr(jac,WeakDeriv);
1207  }
1208  }
1209  break;
1210  case StdRegions::eLaplacian:
1211  {
1212  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1213  (mkey.GetNVarCoeff() > 0)||(mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio)))
1214  {
1215  NekDouble one = 1.0;
1216  DNekMatSharedPtr mat = GenMatrix(mkey);
1217 
1218  returnval = MemoryManager<DNekScalMat>
1219  ::AllocateSharedPtr(one,mat);
1220  }
1221  else
1222  {
1223  MatrixKey lap00key(StdRegions::eLaplacian00,
1224  mkey.GetShapeType(), *this);
1225  MatrixKey lap01key(StdRegions::eLaplacian01,
1226  mkey.GetShapeType(), *this);
1227  MatrixKey lap02key(StdRegions::eLaplacian02,
1228  mkey.GetShapeType(), *this);
1229  MatrixKey lap11key(StdRegions::eLaplacian11,
1230  mkey.GetShapeType(), *this);
1231  MatrixKey lap12key(StdRegions::eLaplacian12,
1232  mkey.GetShapeType(), *this);
1233  MatrixKey lap22key(StdRegions::eLaplacian22,
1234  mkey.GetShapeType(), *this);
1235 
1236  DNekMat &lap00 = *GetStdMatrix(lap00key);
1237  DNekMat &lap01 = *GetStdMatrix(lap01key);
1238  DNekMat &lap02 = *GetStdMatrix(lap02key);
1239  DNekMat &lap11 = *GetStdMatrix(lap11key);
1240  DNekMat &lap12 = *GetStdMatrix(lap12key);
1241  DNekMat &lap22 = *GetStdMatrix(lap22key);
1242 
1243  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1244  Array<TwoD, const NekDouble> gmat
1245  = m_metricinfo->GetGmat(ptsKeys);
1246 
1247  int rows = lap00.GetRows();
1248  int cols = lap00.GetColumns();
1249 
1250  DNekMatSharedPtr lap = MemoryManager<DNekMat>
1251  ::AllocateSharedPtr(rows,cols);
1252 
1253  (*lap) = gmat[0][0]*lap00
1254  + gmat[4][0]*lap11
1255  + gmat[8][0]*lap22
1256  + gmat[3][0]*(lap01 + Transpose(lap01))
1257  + gmat[6][0]*(lap02 + Transpose(lap02))
1258  + gmat[7][0]*(lap12 + Transpose(lap12));
1259 
1260  returnval = MemoryManager<DNekScalMat>
1261  ::AllocateSharedPtr(jac,lap);
1262  }
1263  }
1264  break;
1265  case StdRegions::eHelmholtz:
1266  {
1267  NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorLambda);
1268  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1269  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1270  MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1271  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1272 
1273  int rows = LapMat.GetRows();
1274  int cols = LapMat.GetColumns();
1275 
1276  DNekMatSharedPtr helm = MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
1277 
1278  NekDouble one = 1.0;
1279  (*helm) = LapMat + factor*MassMat;
1280 
1281  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, helm);
1282  }
1283  break;
1284  case StdRegions::eIProductWRTBase:
1285  {
1286  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1287  {
1288  NekDouble one = 1.0;
1289  DNekMatSharedPtr mat = GenMatrix(mkey);
1290  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1291  }
1292  else
1293  {
1294  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1295  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1296  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1297  }
1298  }
1299  break;
1300  case StdRegions::eHybridDGHelmholtz:
1301  case StdRegions::eHybridDGLamToU:
1302  case StdRegions::eHybridDGLamToQ0:
1303  case StdRegions::eHybridDGLamToQ1:
1304  case StdRegions::eHybridDGLamToQ2:
1305  case StdRegions::eHybridDGHelmBndLam:
1306  case StdRegions::eInvLaplacianWithUnityMean:
1307  {
1308  NekDouble one = 1.0;
1309 
1310  DNekMatSharedPtr mat = GenMatrix(mkey);
1311  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1312  }
1313  break;
1314  case StdRegions::eInvHybridDGHelmholtz:
1315  {
1316  NekDouble one = 1.0;
1317 
1318  MatrixKey hkey(StdRegions::eHybridDGHelmholtz, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1319  DNekMatSharedPtr mat = GenMatrix(hkey);
1320 
1321  mat->Invert();
1322  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1323  }
1324  break;
1325  case StdRegions::ePreconLinearSpace:
1326  {
1327  NekDouble one = 1.0;
1328  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1329  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1330  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1331  DNekMatSharedPtr R=BuildVertexMatrix(A);
1332 
1333  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1334  }
1335  break;
1336  case StdRegions::ePreconLinearSpaceMass:
1337  {
1338  NekDouble one = 1.0;
1339  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1340  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1341  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1342  DNekMatSharedPtr R=BuildVertexMatrix(A);
1343 
1344  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1345  }
1346  break;
1347  case StdRegions::ePreconR:
1348  {
1349  NekDouble one = 1.0;
1350  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this,mkey.GetConstFactors(), mkey.GetVarCoeffs());
1351  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1352  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1353 
1354  DNekScalMatSharedPtr Atmp;
1355  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1356 
1357  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1358  }
1359  break;
1360  case StdRegions::ePreconRMass:
1361  {
1362  NekDouble one = 1.0;
1363  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1364  DNekScalBlkMatSharedPtr StatCond = GetLocStaticCondMatrix(masskey);
1365  DNekScalMatSharedPtr A =StatCond->GetBlock(0,0);
1366 
1367  DNekScalMatSharedPtr Atmp;
1368  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1369 
1370  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1371  }
1372  break;
1373  default:
1374  {
1375  //ASSERTL0(false, "Missing definition for " + (*StdRegions::MatrixTypeMap[mkey.GetMatrixType()]));
1376  NekDouble one = 1.0;
1377  DNekMatSharedPtr mat = GenMatrix(mkey);
1378 
1379  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1380  }
1381  break;
1382  }
1383 
1384  return returnval;
1385  }
1386 
1387 
1388  DNekScalBlkMatSharedPtr TetExp::CreateStaticCondMatrix(
1389  const MatrixKey &mkey)
1390  {
1391  DNekScalBlkMatSharedPtr returnval;
1392 
1393  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
1394 
1395  // set up block matrix system
1396  unsigned int nbdry = NumBndryCoeffs();
1397  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
1398  unsigned int exp_size[] = {nbdry, nint};
1399  unsigned int nblks = 2;
1400  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks, nblks, exp_size, exp_size);
1401 
1402  NekDouble factor = 1.0;
1403  MatrixStorage AMatStorage = eFULL;
1404 
1405  switch(mkey.GetMatrixType())
1406  {
1407  case StdRegions::eLaplacian:
1408  case StdRegions::eHelmholtz: // special case since Helmholtz not defined in StdRegions
1409  // use Deformed case for both regular and deformed geometries
1410  factor = 1.0;
1411  goto UseLocRegionsMatrix;
1412  break;
1413  default:
1414  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1415  mkey.GetNVarCoeff())
1416  {
1417  factor = 1.0;
1418  goto UseLocRegionsMatrix;
1419  }
1420  else
1421  {
1422  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1423  factor = mat->Scale();
1424  goto UseStdRegionsMatrix;
1425  }
1426  break;
1427  UseStdRegionsMatrix:
1428  {
1429  NekDouble invfactor = 1.0/factor;
1430  NekDouble one = 1.0;
1431  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
1432  DNekScalMatSharedPtr Atmp;
1433  DNekMatSharedPtr Asubmat;
1434 
1435  //TODO: check below
1436  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
1437  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
1438  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
1439  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
1440  }
1441  break;
1442  UseLocRegionsMatrix:
1443  {
1444  int i,j;
1445  NekDouble invfactor = 1.0/factor;
1446  NekDouble one = 1.0;
1447  DNekScalMat &mat = *GetLocMatrix(mkey);
1448  DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry,AMatStorage);
1449  DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
1450  DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
1451  DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
1452 
1453  Array<OneD,unsigned int> bmap(nbdry);
1454  Array<OneD,unsigned int> imap(nint);
1455  GetBoundaryMap(bmap);
1456  GetInteriorMap(imap);
1457 
1458  for(i = 0; i < nbdry; ++i)
1459  {
1460  for(j = 0; j < nbdry; ++j)
1461  {
1462  (*A)(i,j) = mat(bmap[i],bmap[j]);
1463  }
1464 
1465  for(j = 0; j < nint; ++j)
1466  {
1467  (*B)(i,j) = mat(bmap[i],imap[j]);
1468  }
1469  }
1470 
1471  for(i = 0; i < nint; ++i)
1472  {
1473  for(j = 0; j < nbdry; ++j)
1474  {
1475  (*C)(i,j) = mat(imap[i],bmap[j]);
1476  }
1477 
1478  for(j = 0; j < nint; ++j)
1479  {
1480  (*D)(i,j) = mat(imap[i],imap[j]);
1481  }
1482  }
1483 
1484  // Calculate static condensed system
1485  if(nint)
1486  {
1487  D->Invert();
1488  (*B) = (*B)*(*D);
1489  (*A) = (*A) - (*B)*(*C);
1490  }
1491 
1492  DNekScalMatSharedPtr Atmp;
1493 
1494  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
1495  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
1496  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
1497  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
1498 
1499  }
1500  break;
1501  }
1502  return returnval;
1503  }
1504 
1505 
1506  DNekMatSharedPtr TetExp::v_CreateStdMatrix(
1507  const StdRegions::StdMatrixKey &mkey)
1508  {
1509  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1510  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1511  LibUtilities::BasisKey bkey2 = m_base[2]->GetBasisKey();
1512  StdRegions::StdTetExpSharedPtr tmp = MemoryManager<StdTetExp>::AllocateSharedPtr(bkey0, bkey1, bkey2);
1513 
1514  return tmp->GetStdMatrix(mkey);
1515  }
1516 
1517  DNekScalMatSharedPtr TetExp::v_GetLocMatrix(const MatrixKey &mkey)
1518  {
1519  return m_matrixManager[mkey];
1520  }
1521 
1522  DNekScalBlkMatSharedPtr TetExp::v_GetLocStaticCondMatrix(const MatrixKey &mkey)
1523  {
1524  return m_staticCondMatrixManager[mkey];
1525  }
1526 
1527  void TetExp::v_DropLocStaticCondMatrix(const MatrixKey &mkey)
1528  {
1529  m_staticCondMatrixManager.DeleteObject(mkey);
1530  }
1531 
1532  void TetExp::GeneralMatrixOp_MatOp(
1533  const Array<OneD, const NekDouble> &inarray,
1534  Array<OneD,NekDouble> &outarray,
1535  const StdRegions::StdMatrixKey &mkey)
1536  {
1537  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1538 
1539  if(inarray.get() == outarray.get())
1540  {
1541  Array<OneD,NekDouble> tmp(m_ncoeffs);
1542  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1543 
1544  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1545  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1546  }
1547  else
1548  {
1549  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1550  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1551  }
1552  }
1553 
1554 
1555  void TetExp::v_LaplacianMatrixOp_MatFree_Kernel(
1556  const Array<OneD, const NekDouble> &inarray,
1557  Array<OneD, NekDouble> &outarray,
1558  Array<OneD, NekDouble> &wsp)
1559  {
1560  // This implementation is only valid when there are no
1561  // coefficients associated to the Laplacian operator
1562  if (m_metrics.count(eMetricLaplacian00) == 0)
1563  {
1564  ComputeLaplacianMetric();
1565  }
1566 
1567  int nquad0 = m_base[0]->GetNumPoints();
1568  int nquad1 = m_base[1]->GetNumPoints();
1569  int nquad2 = m_base[2]->GetNumPoints();
1570  int nqtot = nquad0*nquad1*nquad2;
1571 
1572  ASSERTL1(wsp.size() >= 6*nqtot,
1573  "Insufficient workspace size.");
1574  ASSERTL1(m_ncoeffs <= nqtot,
1575  "Workspace not set up for ncoeffs > nqtot");
1576 
1577  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
1578  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
1579  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata();
1580  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
1581  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
1582  const Array<OneD, const NekDouble>& dbase2 = m_base[2]->GetDbdata();
1583  const Array<OneD, const NekDouble>& metric00 = m_metrics[eMetricLaplacian00];
1584  const Array<OneD, const NekDouble>& metric01 = m_metrics[eMetricLaplacian01];
1585  const Array<OneD, const NekDouble>& metric02 = m_metrics[eMetricLaplacian02];
1586  const Array<OneD, const NekDouble>& metric11 = m_metrics[eMetricLaplacian11];
1587  const Array<OneD, const NekDouble>& metric12 = m_metrics[eMetricLaplacian12];
1588  const Array<OneD, const NekDouble>& metric22 = m_metrics[eMetricLaplacian22];
1589 
1590  // Allocate temporary storage
1591  Array<OneD,NekDouble> wsp0 (2*nqtot, wsp);
1592  Array<OneD,NekDouble> wsp1 ( nqtot, wsp+1*nqtot);
1593  Array<OneD,NekDouble> wsp2 ( nqtot, wsp+2*nqtot);
1594  Array<OneD,NekDouble> wsp3 ( nqtot, wsp+3*nqtot);
1595  Array<OneD,NekDouble> wsp4 ( nqtot, wsp+4*nqtot);
1596  Array<OneD,NekDouble> wsp5 ( nqtot, wsp+5*nqtot);
1597 
1598  // LAPLACIAN MATRIX OPERATION
1599  // wsp1 = du_dxi1 = D_xi1 * inarray = D_xi1 * u
1600  // wsp2 = du_dxi2 = D_xi2 * inarray = D_xi2 * u
1601  // wsp2 = du_dxi3 = D_xi3 * inarray = D_xi3 * u
1602  StdExpansion3D::PhysTensorDeriv(inarray,wsp0,wsp1,wsp2);
1603 
1604  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1605  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1606  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1607  // especially for this purpose
1608  Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp0[0],1,&metric01[0],1,&wsp1[0],1,&wsp3[0],1);
1609  Vmath::Vvtvp (nqtot,&metric02[0],1,&wsp2[0],1,&wsp3[0],1,&wsp3[0],1);
1610  Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp0[0],1,&metric11[0],1,&wsp1[0],1,&wsp4[0],1);
1611  Vmath::Vvtvp (nqtot,&metric12[0],1,&wsp2[0],1,&wsp4[0],1,&wsp4[0],1);
1612  Vmath::Vvtvvtp(nqtot,&metric02[0],1,&wsp0[0],1,&metric12[0],1,&wsp1[0],1,&wsp5[0],1);
1613  Vmath::Vvtvp (nqtot,&metric22[0],1,&wsp2[0],1,&wsp5[0],1,&wsp5[0],1);
1614 
1615  // outarray = m = (D_xi1 * B)^T * k
1616  // wsp1 = n = (D_xi2 * B)^T * l
1617  IProductWRTBase_SumFacKernel(dbase0,base1,base2,wsp3,outarray,wsp0,false,true,true);
1618  IProductWRTBase_SumFacKernel(base0,dbase1,base2,wsp4,wsp2, wsp0,true,false,true);
1619  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
1620  IProductWRTBase_SumFacKernel(base0,base1,dbase2,wsp5,wsp2, wsp0,true,true,false);
1621  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
1622  }
1623 
1624 
1625  void TetExp::v_ComputeLaplacianMetric()
1626  {
1627  if (m_metrics.count(eMetricQuadrature) == 0)
1628  {
1629  ComputeQuadratureMetric();
1630  }
1631 
1632  int i, j;
1633  const unsigned int nqtot = GetTotPoints();
1634  const unsigned int dim = 3;
1635  const MetricType m[3][3] = { {eMetricLaplacian00, eMetricLaplacian01, eMetricLaplacian02},
1636  {eMetricLaplacian01, eMetricLaplacian11, eMetricLaplacian12},
1637  {eMetricLaplacian02, eMetricLaplacian12, eMetricLaplacian22}
1638  };
1639 
1640  for (unsigned int i = 0; i < dim; ++i)
1641  {
1642  for (unsigned int j = i; j < dim; ++j)
1643  {
1644  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1645  }
1646  }
1647 
1648  // Define shorthand synonyms for m_metrics storage
1649  Array<OneD,NekDouble> g0 (m_metrics[m[0][0]]);
1650  Array<OneD,NekDouble> g1 (m_metrics[m[1][1]]);
1651  Array<OneD,NekDouble> g2 (m_metrics[m[2][2]]);
1652  Array<OneD,NekDouble> g3 (m_metrics[m[0][1]]);
1653  Array<OneD,NekDouble> g4 (m_metrics[m[0][2]]);
1654  Array<OneD,NekDouble> g5 (m_metrics[m[1][2]]);
1655 
1656  // Allocate temporary storage
1657  Array<OneD,NekDouble> alloc(7*nqtot,0.0);
1658  Array<OneD,NekDouble> h0 (alloc); // h0
1659  Array<OneD,NekDouble> h1 (alloc+ 1*nqtot);// h1
1660  Array<OneD,NekDouble> h2 (alloc+ 2*nqtot);// h2
1661  Array<OneD,NekDouble> h3 (alloc+ 3*nqtot);// h3
1662  Array<OneD,NekDouble> wsp4 (alloc+ 4*nqtot);// wsp4
1663  Array<OneD,NekDouble> wsp5 (alloc+ 5*nqtot);// wsp5
1664  Array<OneD,NekDouble> wsp6 (alloc+ 6*nqtot);// wsp6
1665  // Reuse some of the storage as workspace
1666  Array<OneD,NekDouble> wsp7 (alloc); // wsp7
1667  Array<OneD,NekDouble> wsp8 (alloc+ 1*nqtot);// wsp8
1668  Array<OneD,NekDouble> wsp9 (alloc+ 2*nqtot);// wsp9
1669 
1670  const Array<TwoD, const NekDouble>& df =
1671  m_metricinfo->GetDerivFactors(GetPointsKeys());
1672  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
1673  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
1674  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1675  const unsigned int nquad0 = m_base[0]->GetNumPoints();
1676  const unsigned int nquad1 = m_base[1]->GetNumPoints();
1677  const unsigned int nquad2 = m_base[2]->GetNumPoints();
1678 
1679  for(j = 0; j < nquad2; ++j)
1680  {
1681  for(i = 0; i < nquad1; ++i)
1682  {
1683  Vmath::Fill(nquad0, 4.0/(1.0-z1[i])/(1.0-z2[j]), &h0[0]+i*nquad0 + j*nquad0*nquad1,1);
1684  Vmath::Fill(nquad0, 2.0/(1.0-z1[i])/(1.0-z2[j]), &h1[0]+i*nquad0 + j*nquad0*nquad1,1);
1685  Vmath::Fill(nquad0, 2.0/(1.0-z2[j]), &h2[0]+i*nquad0 + j*nquad0*nquad1,1);
1686  Vmath::Fill(nquad0, (1.0+z1[i])/(1.0-z2[j]), &h3[0]+i*nquad0 + j*nquad0*nquad1,1);
1687  }
1688  }
1689  for(i = 0; i < nquad0; i++)
1690  {
1691  Blas::Dscal(nquad1*nquad2, 1+z0[i], &h1[0]+i, nquad0);
1692  }
1693 
1694  // Step 3. Construct combined metric terms for physical space to
1695  // collapsed coordinate system.
1696  // Order of construction optimised to minimise temporary storage
1697  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1698  {
1699  // wsp4
1700  Vmath::Vadd(nqtot, &df[1][0], 1, &df[2][0], 1, &wsp4[0], 1);
1701  Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &wsp4[0], 1, &h1[0], 1, &wsp4[0], 1);
1702  // wsp5
1703  Vmath::Vadd(nqtot, &df[4][0], 1, &df[5][0], 1, &wsp5[0], 1);
1704  Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &wsp5[0], 1, &h1[0], 1, &wsp5[0], 1);
1705  // wsp6
1706  Vmath::Vadd(nqtot, &df[7][0], 1, &df[8][0], 1, &wsp6[0], 1);
1707  Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &wsp6[0], 1, &h1[0], 1, &wsp6[0], 1);
1708 
1709  // g0
1710  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0], 1, &g0[0], 1);
1711  Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1712 
1713  // g4
1714  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0], 1, &g4[0], 1);
1715  Vmath::Vvtvp (nqtot, &df[8][0], 1, &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1716 
1717  // overwrite h0, h1, h2
1718  // wsp7 (h2f1 + h3f2)
1719  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &h2[0], 1, &df[2][0], 1, &h3[0], 1, &wsp7[0], 1);
1720  // wsp8 (h2f4 + h3f5)
1721  Vmath::Vvtvvtp(nqtot, &df[4][0], 1, &h2[0], 1, &df[5][0], 1, &h3[0], 1, &wsp8[0], 1);
1722  // wsp9 (h2f7 + h3f8)
1723  Vmath::Vvtvvtp(nqtot, &df[7][0], 1, &h2[0], 1, &df[8][0], 1, &h3[0], 1, &wsp9[0], 1);
1724 
1725  // g3
1726  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp7[0], 1, &wsp5[0], 1, &wsp8[0], 1, &g3[0], 1);
1727  Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp9[0], 1, &g3[0], 1, &g3[0], 1);
1728 
1729  // overwrite wsp4, wsp5, wsp6
1730  // g1
1731  Vmath::Vvtvvtp(nqtot, &wsp7[0], 1, &wsp7[0], 1, &wsp8[0], 1, &wsp8[0], 1, &g1[0], 1);
1732  Vmath::Vvtvp (nqtot, &wsp9[0], 1, &wsp9[0], 1, &g1[0], 1, &g1[0], 1);
1733 
1734  // g5
1735  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp7[0], 1, &df[5][0], 1, &wsp8[0], 1, &g5[0], 1);
1736  Vmath::Vvtvp (nqtot, &df[8][0], 1, &wsp9[0], 1, &g5[0], 1, &g5[0], 1);
1737 
1738  // g2
1739  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &df[2][0], 1, &df[5][0], 1, &df[5][0], 1, &g2[0], 1);
1740  Vmath::Vvtvp (nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1741  }
1742  else
1743  {
1744  // wsp4
1745  Vmath::Svtsvtp(nqtot, df[0][0], &h0[0], 1, df[1][0] + df[2][0], &h1[0], 1, &wsp4[0], 1);
1746  // wsp5
1747  Vmath::Svtsvtp(nqtot, df[3][0], &h0[0], 1, df[4][0] + df[5][0], &h1[0], 1, &wsp5[0], 1);
1748  // wsp6
1749  Vmath::Svtsvtp(nqtot, df[6][0], &h0[0], 1, df[7][0] + df[8][0], &h1[0], 1, &wsp6[0], 1);
1750 
1751  // g0
1752  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0], 1, &g0[0], 1);
1753  Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1754 
1755  // g4
1756  Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1, &g4[0], 1);
1757  Vmath::Svtvp (nqtot, df[8][0], &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1758 
1759  // overwrite h0, h1, h2
1760  // wsp7 (h2f1 + h3f2)
1761  Vmath::Svtsvtp(nqtot, df[1][0], &h2[0], 1, df[2][0], &h3[0], 1, &wsp7[0], 1);
1762  // wsp8 (h2f4 + h3f5)
1763  Vmath::Svtsvtp(nqtot, df[4][0], &h2[0], 1, df[5][0], &h3[0], 1, &wsp8[0], 1);
1764  // wsp9 (h2f7 + h3f8)
1765  Vmath::Svtsvtp(nqtot, df[7][0], &h2[0], 1, df[8][0], &h3[0], 1, &wsp9[0], 1);
1766 
1767  // g3
1768  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp7[0], 1, &wsp5[0], 1, &wsp8[0], 1, &g3[0], 1);
1769  Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp9[0], 1, &g3[0], 1, &g3[0], 1);
1770 
1771  // overwrite wsp4, wsp5, wsp6
1772  // g1
1773  Vmath::Vvtvvtp(nqtot, &wsp7[0], 1, &wsp7[0], 1, &wsp8[0], 1, &wsp8[0], 1, &g1[0], 1);
1774  Vmath::Vvtvp (nqtot, &wsp9[0], 1, &wsp9[0], 1, &g1[0], 1, &g1[0], 1);
1775 
1776  // g5
1777  Vmath::Svtsvtp(nqtot, df[2][0], &wsp7[0], 1, df[5][0], &wsp8[0], 1, &g5[0], 1);
1778  Vmath::Svtvp (nqtot, df[8][0], &wsp9[0], 1, &g5[0], 1, &g5[0], 1);
1779 
1780  // g2
1781  Vmath::Fill(nqtot, df[2][0]*df[2][0] + df[5][0]*df[5][0] + df[8][0]*df[8][0], &g2[0], 1);
1782  }
1783 
1784  for (unsigned int i = 0; i < dim; ++i)
1785  {
1786  for (unsigned int j = i; j < dim; ++j)
1787  {
1788  MultiplyByQuadratureMetric(m_metrics[m[i][j]],
1789  m_metrics[m[i][j]]);
1790 
1791  }
1792  }
1793 
1794 
1795  }
1796  }//end of namespace
1797 }//end of namespace
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:274
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::LibUtilities::BasisKey::GetNumPoints
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:133
Nektar::LibUtilities::BasisKey::GetPointsKey
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:150
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:60
Nektar::LocalRegions::Expansion3D::m_faceNormals
std::map< int, NormalVector > m_faceNormals
Definition: Expansion3D.h:123
Nektar::LocalRegions::Expansion3D::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: Expansion3D.cpp:451
Nektar::LocalRegions::Expansion::m_elmtBndNormDirElmtLen
std::map< int, Array< OneD, NekDouble > > m_elmtBndNormDirElmtLen
the element length in each element boundary(Vertex, edge or face) normal direction calculated based o...
Definition: Expansion.h:284
Nektar::LocalRegions::Expansion::GetGeom
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:172
Nektar::LocalRegions::Expansion::m_geom
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:272
Nektar::LocalRegions::Expansion::BuildVertexMatrix
DNekMatSharedPtr BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
Definition: Expansion.cpp:103
Nektar::LocalRegions::Expansion::m_metricinfo
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:273
Nektar::LocalRegions::Expansion::v_GetCoords
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: Expansion.cpp:318
Nektar::LocalRegions::Expansion::GetLocMatrix
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:90
Nektar::LocalRegions::Expansion::BuildTransformationMatrix
DNekMatSharedPtr BuildTransformationMatrix(const DNekScalMatSharedPtr &r_bnd, const StdRegions::MatrixType matrixType)
Definition: Expansion.cpp:95
Nektar::LocalRegions::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:526
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::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
Nektar::LocalRegions::TetExp::v_ComputeLaplacianMetric
virtual void v_ComputeLaplacianMetric()
Definition: TetExp.cpp:1625
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:545
Nektar::LocalRegions::TetExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: TetExp.cpp:1081
Nektar::LocalRegions::TetExp::m_staticCondMatrixManager
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: TetExp.h:216
Nektar::LocalRegions::TetExp::v_GetStdExp
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
Definition: TetExp.cpp:566
Nektar::LocalRegions::TetExp::m_matrixManager
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: TetExp.h:215
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:594
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:1029
Nektar::LocalRegions::TetExp::CreateStaticCondMatrix
DNekScalBlkMatSharedPtr CreateStaticCondMatrix(const MatrixKey &mkey)
Definition: TetExp.cpp:1388
Nektar::LocalRegions::TetExp::v_GetCoordim
virtual int v_GetCoordim()
Definition: TetExp.cpp:589
Nektar::LocalRegions::TetExp::v_GetLocMatrix
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: TetExp.cpp:1517
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:1020
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::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: TetExp.cpp:1506
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::v_StdPhysEvaluate
virtual NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
Definition: TetExp.cpp:496
Nektar::LocalRegions::TetExp::v_GetLinStdExp
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
Definition: TetExp.cpp:575
Nektar::LocalRegions::TetExp::v_AlignVectorToCollapsedDir
virtual void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: TetExp.cpp:419
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::LocalRegions::TetExp::v_ComputeTraceNormal
void v_ComputeTraceNormal(const int face)
Compute the normal of a triangular face.
Definition: TetExp.cpp:751
Nektar::LocalRegions::TetExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
Definition: TetExp.cpp:1048
Nektar::LocalRegions::TetExp::v_DetShapeType
virtual LibUtilities::ShapeType v_DetShapeType() const
Return Shape of region, using ShapeType enum list.
Definition: TetExp.cpp:561
Nektar::LocalRegions::TetExp::v_GetLocStaticCondMatrix
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: TetExp.cpp:1522
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:1532
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::LocalRegions::TetExp::v_GetTracePhysMap
virtual void v_GetTracePhysMap(const int face, Array< OneD, int > &outarray)
Returns the physical values at the quadrature points of a face.
Definition: TetExp.cpp:661
Nektar::LocalRegions::TetExp::v_DropLocStaticCondMatrix
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: TetExp.cpp:1527
Nektar::LocalRegions::TetExp::v_PhysEvaluate
virtual NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals)
Definition: TetExp.cpp:508
Nektar::LocalRegions::TetExp::CreateMatrix
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition: TetExp.cpp:1105
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:1555
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:84
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:171
Nektar::StdRegions::StdExpansion3D::v_LaplacianMatrixOp_MatFree
virtual void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: StdExpansion3D.cpp:305
Nektar::StdRegions::StdExpansion3D::IProductWRTBase_SumFacKernel
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
Definition: StdExpansion3D.cpp:128
Nektar::StdRegions::StdExpansion3D::v_HelmholtzMatrixOp_MatFree
virtual void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: StdExpansion3D.cpp:347
Nektar::StdRegions::StdExpansion::GetBoundaryMap
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:687
Nektar::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:124
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:134
Nektar::StdRegions::StdExpansion::GetBasisType
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:158
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:617
Nektar::StdRegions::StdExpansion::GetPointsKeys
const LibUtilities::PointsKeyVector GetPointsKeys() const
Definition: StdExpansion.h:1078
Nektar::StdRegions::StdExpansion::GetStdStaticCondMatrix
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:622
Nektar::StdRegions::StdExpansion::GetLocStaticCondMatrix
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
Definition: StdExpansion.h:660
Nektar::StdRegions::StdExpansion::GetTraceBasisKey
const LibUtilities::BasisKey GetTraceBasisKey(const int i, int k=-1) const
This function returns the basis key belonging to the i-th trace.
Definition: StdExpansion.h:304
Nektar::StdRegions::StdExpansion::IProductWRTBase
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
this function calculates the inner product of a given function f with the different modes of the expa...
Definition: StdExpansion.h:537
Nektar::StdRegions::StdExpansion::GetInteriorMap
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:692
Nektar::StdRegions::StdExpansion::GenMatrix
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:850
Nektar::StdRegions::StdExpansion::GetPointsType
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:208
Nektar::StdRegions::StdExpansion::GetNumPoints
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:221
Nektar::StdRegions::StdExpansion::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:733
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1162
Nektar::StdRegions::StdMatrixKey::GetConstFactors
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:135
Nektar::StdRegions::StdMatrixKey::GetShapeType
LibUtilities::ShapeType GetShapeType() const
Definition: StdMatrixKey.h:86
Nektar::StdRegions::StdMatrixKey::GetVarCoeffs
const VarCoeffMap & GetVarCoeffs() const
Definition: StdMatrixKey.h:161
Nektar::StdRegions::StdMatrixKey::GetMatrixType
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:81
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:121
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:130
Nektar::StdRegions::StdTetExp::DetShapeType
LibUtilities::ShapeType DetShapeType() const
Definition: StdTetExp.h:69
Blas::Dgemv
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = A x where A[m x n].
Definition: Blas.hpp:265
Blas::Dscal
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:182
Blas::Daxpy
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
Definition: Blas.hpp:167
Nektar::LibUtilities::Interp2D
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
Definition: Interp.cpp:115
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
Nektar::LibUtilities::eModified_B
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::LibUtilities::eModified_C
@ eModified_C
Principle Modified Functions .
Definition: BasisType.h:50
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
Nektar::SpatialDomains::GeomFactorsSharedPtr
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
Nektar::SpatialDomains::GeomType
GeomType
Indicates the type of element geometry.
Definition: SpatialDomains.hpp:52
Nektar::SpatialDomains::eRegular
@ eRegular
Geometry is straight-sided with constant geometric factors.
Definition: SpatialDomains.hpp:54
Nektar::SpatialDomains::eNoGeomType
@ eNoGeomType
No type defined.
Definition: SpatialDomains.hpp:53
Nektar::SpatialDomains::eMovingRegular
@ eMovingRegular
Currently unused.
Definition: SpatialDomains.hpp:57
Nektar::SpatialDomains::eDeformed
@ eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:56
Nektar::SpatialDomains::TetGeomSharedPtr
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:82
Nektar::StdRegions::StdTetExpSharedPtr
std::shared_ptr< StdTetExp > StdTetExpSharedPtr
Definition: StdTetExp.h:279
Nektar::StdRegions::StdExpansionSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1714
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::DNekScalMatSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:69
Nektar::DNekBlkMatSharedPtr
std::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:71
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:217
Nektar::DNekScalBlkMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:73
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:475
Vmath::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:691
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:192
Vmath::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:565
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:513
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:322
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
Vmath::Sdiv
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:291
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:257
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:436
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
Vmath::Vvtvvtp
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:625
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267
Generated on Wed Nov 24 2021 10:33:29 for Nektar++ by doxygen 1.9.1