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