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