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.num_elements())
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.num_elements())
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.num_elements())
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.num_elements())
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.num_elements())
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.num_elements())
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  int i;
364 
365  const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
366  const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
367 
368  Array<OneD, NekDouble> gfac0(nquad0 );
369  Array<OneD, NekDouble> gfac2(nquad2 );
370  Array<OneD, NekDouble> tmp1 (nqtot );
371  Array<OneD, NekDouble> tmp2 (nqtot );
372  Array<OneD, NekDouble> tmp3 (nqtot );
373  Array<OneD, NekDouble> tmp4 (nqtot );
374  Array<OneD, NekDouble> tmp5 (nqtot );
375  Array<OneD, NekDouble> tmp6 (m_ncoeffs);
376  Array<OneD, NekDouble> wsp (order0*nquad2*(nquad1+order1));
377 
378  const Array<TwoD, const NekDouble>& df =
379  m_metricinfo->GetDerivFactors(GetPointsKeys());
380 
381  MultiplyByQuadratureMetric(inarray, tmp1);
382 
383  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
384  {
385  Vmath::Vmul(nqtot,&df[3*dir][0], 1,tmp1.get(),1,tmp2.get(),1);
386  Vmath::Vmul(nqtot,&df[3*dir+1][0],1,tmp1.get(),1,tmp3.get(),1);
387  Vmath::Vmul(nqtot,&df[3*dir+2][0],1,tmp1.get(),1,tmp4.get(),1);
388  }
389  else
390  {
391  Vmath::Smul(nqtot, df[3*dir][0], tmp1.get(),1,tmp2.get(), 1);
392  Vmath::Smul(nqtot, df[3*dir+1][0],tmp1.get(),1,tmp3.get(), 1);
393  Vmath::Smul(nqtot, df[3*dir+2][0],tmp1.get(),1,tmp4.get(), 1);
394  }
395 
396  // set up geometric factor: (1+z0)/2
397  for (i = 0; i < nquad0; ++i)
398  {
399  gfac0[i] = 0.5*(1+z0[i]);
400  }
401 
402  // Set up geometric factor: 2/(1-z2)
403  for (i = 0; i < nquad2; ++i)
404  {
405  gfac2[i] = 2.0/(1-z2[i]);
406  }
407 
408  const int nq01 = nquad0*nquad1;
409 
410  for (i = 0; i < nquad2; ++i)
411  {
412  Vmath::Smul(nq01,gfac2[i],&tmp2[0]+i*nq01,1,&tmp2[0]+i*nq01,1);
413  Vmath::Smul(nq01,gfac2[i],&tmp4[0]+i*nq01,1,&tmp5[0]+i*nq01,1);
414  }
415 
416  for(i = 0; i < nquad1*nquad2; ++i)
417  {
418  Vmath::Vmul(nquad0,&gfac0[0],1,&tmp5[0]+i*nquad0,1,
419  &tmp5[0]+i*nquad0,1);
420  }
421 
422  Vmath::Vadd(nqtot, &tmp2[0], 1, &tmp5[0], 1, &tmp2[0], 1);
423 
424  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
425  m_base[1]->GetBdata (),
426  m_base[2]->GetBdata (),
427  tmp2,outarray,wsp,
428  true,true,true);
429 
430  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata (),
431  m_base[1]->GetDbdata(),
432  m_base[2]->GetBdata (),
433  tmp3,tmp6,wsp,
434  true,true,true);
435 
436  Vmath::Vadd(m_ncoeffs, tmp6, 1, outarray, 1, outarray, 1);
437 
438  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata (),
439  m_base[1]->GetBdata (),
440  m_base[2]->GetDbdata(),
441  tmp4,tmp6,wsp,
442  true,true,true);
443 
444  Vmath::Vadd(m_ncoeffs, tmp6, 1, outarray, 1, outarray, 1);
445  }
446 
447  //---------------------------------------
448  // Evaluation functions
449  //---------------------------------------
450 
451  StdRegions::StdExpansionSharedPtr PrismExp::v_GetStdExp(void) const
452  {
453  return MemoryManager<StdRegions::StdPrismExp>
454  ::AllocateSharedPtr(m_base[0]->GetBasisKey(),
455  m_base[1]->GetBasisKey(),
456  m_base[2]->GetBasisKey());
457  }
458 
459  StdRegions::StdExpansionSharedPtr PrismExp::v_GetLinStdExp(void) const
460  {
461  LibUtilities::BasisKey bkey0(m_base[0]->GetBasisType(),
462  2, m_base[0]->GetPointsKey());
463  LibUtilities::BasisKey bkey1(m_base[1]->GetBasisType(),
464  2, m_base[1]->GetPointsKey());
465  LibUtilities::BasisKey bkey2(m_base[2]->GetBasisType(),
466  2, m_base[2]->GetPointsKey());
467 
468  return MemoryManager<StdRegions::StdPrismExp>
469  ::AllocateSharedPtr( bkey0, bkey1, bkey2);
470  }
471 
472  /**
473  * @brief Get the coordinates #coords at the local coordinates
474  * #Lcoords.
475  */
476  void PrismExp::v_GetCoord(const Array<OneD, const NekDouble>& Lcoords,
477  Array<OneD, NekDouble>& coords)
478  {
479  int i;
480 
481  ASSERTL1(Lcoords[0] <= -1.0 && Lcoords[0] >= 1.0 &&
482  Lcoords[1] <= -1.0 && Lcoords[1] >= 1.0 &&
483  Lcoords[2] <= -1.0 && Lcoords[2] >= 1.0,
484  "Local coordinates are not in region [-1,1]");
485 
486  m_geom->FillGeom();
487 
488  for(i = 0; i < m_geom->GetCoordim(); ++i)
489  {
490  coords[i] = m_geom->GetCoord(i,Lcoords);
491  }
492  }
493 
494  void PrismExp::v_GetCoords(
495  Array<OneD, NekDouble> &coords_0,
496  Array<OneD, NekDouble> &coords_1,
497  Array<OneD, NekDouble> &coords_2)
498  {
499  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
500  }
501 
502  /**
503  * Given the local cartesian coordinate \a Lcoord evaluate the
504  * value of physvals at this point by calling through to the
505  * StdExpansion method
506  */
507  NekDouble PrismExp::v_StdPhysEvaluate(
508  const Array<OneD, const NekDouble> &Lcoord,
509  const Array<OneD, const NekDouble> &physvals)
510  {
511  // Evaluate point in local (eta) coordinates.
512  return StdPrismExp::v_PhysEvaluate(Lcoord,physvals);
513  }
514 
515  NekDouble PrismExp::v_PhysEvaluate(const Array<OneD, const NekDouble>& coord,
516  const Array<OneD, const NekDouble>& physvals)
517  {
518  Array<OneD, NekDouble> Lcoord(3);
519 
520  ASSERTL0(m_geom,"m_geom not defined");
521 
522  m_geom->GetLocCoords(coord, Lcoord);
523 
524  return StdPrismExp::v_PhysEvaluate(Lcoord, physvals);
525  }
526 
527 
528  //---------------------------------------
529  // Helper functions
530  //---------------------------------------
531 
532  int PrismExp::v_GetCoordim()
533  {
534  return m_geom->GetCoordim();
535  }
536 
537  void PrismExp::v_ExtractDataToCoeffs(
538  const NekDouble* data,
539  const std::vector<unsigned int >& nummodes,
540  const int mode_offset,
541  NekDouble* coeffs,
542  std::vector<LibUtilities::BasisType> &fromType)
543  {
544  boost::ignore_unused(fromType);
545 
546  int data_order0 = nummodes[mode_offset];
547  int fillorder0 = min(m_base[0]->GetNumModes(),data_order0);
548  int data_order1 = nummodes[mode_offset+1];
549  int order1 = m_base[1]->GetNumModes();
550  int fillorder1 = min(order1,data_order1);
551  int data_order2 = nummodes[mode_offset+2];
552  int order2 = m_base[2]->GetNumModes();
553  int fillorder2 = min(order2,data_order2);
554 
555  switch(m_base[0]->GetBasisType())
556  {
557  case LibUtilities::eModified_A:
558  {
559  int i,j;
560  int cnt = 0;
561  int cnt1 = 0;
562 
563  ASSERTL1(m_base[1]->GetBasisType() ==
564  LibUtilities::eModified_A,
565  "Extraction routine not set up for this basis");
566  ASSERTL1(m_base[2]->GetBasisType() ==
567  LibUtilities::eModified_B,
568  "Extraction routine not set up for this basis");
569 
570  Vmath::Zero(m_ncoeffs,coeffs,1);
571  for(j = 0; j < fillorder0; ++j)
572  {
573  for(i = 0; i < fillorder1; ++i)
574  {
575  Vmath::Vcopy(fillorder2-j, &data[cnt], 1,
576  &coeffs[cnt1], 1);
577  cnt += data_order2-j;
578  cnt1 += order2-j;
579  }
580 
581  // count out data for j iteration
582  for(i = fillorder1; i < data_order1; ++i)
583  {
584  cnt += data_order2-j;
585  }
586 
587  for(i = fillorder1; i < order1; ++i)
588  {
589  cnt1 += order2-j;
590  }
591  }
592  }
593  break;
594  default:
595  ASSERTL0(false, "basis is either not set up or not "
596  "hierarchicial");
597  }
598  }
599 
600  void PrismExp::v_GetFacePhysMap(const int face,
601  Array<OneD, int> &outarray)
602  {
603  int nquad0 = m_base[0]->GetNumPoints();
604  int nquad1 = m_base[1]->GetNumPoints();
605  int nquad2 = m_base[2]->GetNumPoints();
606  int nq0 = 0;
607  int nq1 = 0;
608 
609  switch(face)
610  {
611  case 0:
612  nq0 = nquad0;
613  nq1 = nquad1;
614  if(outarray.num_elements()!=nq0*nq1)
615  {
616  outarray = Array<OneD, int>(nq0*nq1);
617  }
618 
619  //Directions A and B positive
620  for(int i = 0; i < nquad0*nquad1; ++i)
621  {
622  outarray[i] = i;
623  }
624  break;
625  case 1:
626 
627  nq0 = nquad0;
628  nq1 = nquad2;
629  if(outarray.num_elements()!=nq0*nq1)
630  {
631  outarray = Array<OneD, int>(nq0*nq1);
632  }
633 
634  //Direction A and B positive
635  for (int k=0; k<nquad2; k++)
636  {
637  for(int i = 0; i < nquad0; ++i)
638  {
639  outarray[k*nquad0+i] = (nquad0*nquad1*k)+i;
640  }
641  }
642 
643  break;
644  case 2:
645 
646  nq0 = nquad1;
647  nq1 = nquad2;
648  if(outarray.num_elements()!=nq0*nq1)
649  {
650  outarray = Array<OneD, int>(nq0*nq1);
651  }
652 
653  //Directions A and B positive
654  for(int j = 0; j < nquad1*nquad2; ++j)
655  {
656  outarray[j] = nquad0-1 + j*nquad0;
657 
658  }
659  break;
660  case 3:
661  nq0 = nquad0;
662  nq1 = nquad2;
663  if(outarray.num_elements()!=nq0*nq1)
664  {
665  outarray = Array<OneD, int>(nq0*nq1);
666  }
667 
668  //Direction A and B positive
669  for (int k=0; k<nquad2; k++)
670  {
671  for(int i = 0; i < nquad0; ++i)
672  {
673  outarray[k*nquad0+i] = nquad0*(nquad1-1) + (nquad0*nquad1*k)+i;
674  }
675  }
676  break;
677  case 4:
678 
679  nq0 = nquad1;
680  nq1 = nquad2;
681  if(outarray.num_elements()!=nq0*nq1)
682  {
683  outarray = Array<OneD, int>(nq0*nq1);
684  }
685 
686  //Directions A and B positive
687  for(int j = 0; j < nquad1*nquad2; ++j)
688  {
689  outarray[j] = j*nquad0;
690 
691  }
692  break;
693  default:
694  ASSERTL0(false,"face value (> 4) is out of range");
695  break;
696  }
697 
698  }
699 
700  /** \brief Get the normals along specficied face
701  * Get the face normals interplated to a points0 x points 0
702  * type distribution
703  **/
704  void PrismExp::v_ComputeFaceNormal(const int face)
705  {
706  const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
707  GetGeom()->GetMetricInfo();
708 
709  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
710  for(int i = 0; i < ptsKeys.size(); ++i)
711  {
712  // Need at least 2 points for computing normals
713  if (ptsKeys[i].GetNumPoints() == 1)
714  {
715  LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
716  ptsKeys[i] = pKey;
717  }
718  }
719 
720  SpatialDomains::GeomType type = geomFactors->GetGtype();
721  const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys);
722  const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
723 
724  int nq0 = ptsKeys[0].GetNumPoints();
725  int nq1 = ptsKeys[1].GetNumPoints();
726  int nq2 = ptsKeys[2].GetNumPoints();
727  int nq01 = nq0*nq1;
728  int nqtot;
729 
730 
731  LibUtilities::BasisKey tobasis0 = DetFaceBasisKey(face,0);
732  LibUtilities::BasisKey tobasis1 = DetFaceBasisKey(face,1);
733 
734  // Number of quadrature points in face expansion.
735  int nq_face = tobasis0.GetNumPoints()*tobasis1.GetNumPoints();
736 
737  int vCoordDim = GetCoordim();
738  int i;
739 
740  m_faceNormals[face] = Array<OneD, Array<OneD, NekDouble> >(vCoordDim);
741  Array<OneD, Array<OneD, NekDouble> > &normal = m_faceNormals[face];
742  for (i = 0; i < vCoordDim; ++i)
743  {
744  normal[i] = Array<OneD, NekDouble>(nq_face);
745  }
746 
747  // Regular geometry case
748  if (type == SpatialDomains::eRegular ||
749  type == SpatialDomains::eMovingRegular)
750  {
751  NekDouble fac;
752  // Set up normals
753  switch(face)
754  {
755  case 0:
756  {
757  for(i = 0; i < vCoordDim; ++i)
758  {
759  normal[i][0] = -df[3*i+2][0];;
760  }
761  break;
762  }
763  case 1:
764  {
765  for(i = 0; i < vCoordDim; ++i)
766  {
767  normal[i][0] = -df[3*i+1][0];
768  }
769  break;
770  }
771  case 2:
772  {
773  for(i = 0; i < vCoordDim; ++i)
774  {
775  normal[i][0] = df[3*i][0]+df[3*i+2][0];
776  }
777  break;
778  }
779  case 3:
780  {
781  for(i = 0; i < vCoordDim; ++i)
782  {
783  normal[i][0] = df[3*i+1][0];
784  }
785  break;
786  }
787  case 4:
788  {
789  for(i = 0; i < vCoordDim; ++i)
790  {
791  normal[i][0] = -df[3*i][0];
792  }
793  break;
794  }
795  default:
796  ASSERTL0(false,"face is out of range (face < 4)");
797  }
798 
799  // Normalise resulting vector.
800  fac = 0.0;
801  for(i = 0; i < vCoordDim; ++i)
802  {
803  fac += normal[i][0]*normal[i][0];
804  }
805  fac = 1.0/sqrt(fac);
806  for (i = 0; i < vCoordDim; ++i)
807  {
808  Vmath::Fill(nq_face,fac*normal[i][0],normal[i],1);
809  }
810  }
811  else
812  {
813  // Set up deformed normals.
814  int j, k;
815 
816  // Determine number of quadrature points on the face of 3D elmt
817  if (face == 0)
818  {
819  nqtot = nq0*nq1;
820  }
821  else if (face == 1 || face == 3)
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 derivatives
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[2][tmp])*jac[tmp];
887  normals[nqtot+j+k*nq1] =
888  (df[3][tmp]+df[5][tmp])*jac[tmp];
889  normals[2*nqtot+j+k*nq1] =
890  (df[6][tmp]+df[8][tmp])*jac[tmp];
891  faceJac[j+k*nq1] = jac[tmp];
892  }
893  }
894 
895  points0 = ptsKeys[1];
896  points1 = ptsKeys[2];
897  break;
898  }
899 
900  case 3:
901  {
902  for (j = 0; j < nq0; ++j)
903  {
904  for(k = 0; k < nq2; ++k)
905  {
906  int tmp = nq0*(nq1-1) + j + nq01*k;
907  normals[j+k*nq0] =
908  df[1][tmp]*jac[tmp];
909  normals[nqtot+j+k*nq0] =
910  df[4][tmp]*jac[tmp];
911  normals[2*nqtot+j+k*nq0] =
912  df[7][tmp]*jac[tmp];
913  faceJac[j+k*nq0] = jac[tmp];
914  }
915  }
916 
917  points0 = ptsKeys[0];
918  points1 = ptsKeys[2];
919  break;
920  }
921 
922  case 4:
923  {
924  for (j = 0; j < nq1; ++j)
925  {
926  for(k = 0; k < nq2; ++k)
927  {
928  int tmp = j*nq0+nq01*k;
929  normals[j+k*nq1] =
930  -df[0][tmp]*jac[tmp];
931  normals[nqtot+j+k*nq1] =
932  -df[3][tmp]*jac[tmp];
933  normals[2*nqtot+j+k*nq1] =
934  -df[6][tmp]*jac[tmp];
935  faceJac[j+k*nq1] = jac[tmp];
936  }
937  }
938 
939  points0 = ptsKeys[1];
940  points1 = ptsKeys[2];
941  break;
942  }
943 
944  default:
945  ASSERTL0(false,"face is out of range (face < 4)");
946  }
947 
948 
949  Array<OneD, NekDouble> work (nq_face, 0.0);
950  // Interpolate Jacobian and invert
951  LibUtilities::Interp2D(points0, points1, faceJac,
952  tobasis0.GetPointsKey(),
953  tobasis1.GetPointsKey(),
954  work);
955  Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
956 
957  // Interpolate normal and multiply by inverse Jacobian.
958  for(i = 0; i < vCoordDim; ++i)
959  {
960  LibUtilities::Interp2D(points0, points1,
961  &normals[i*nqtot],
962  tobasis0.GetPointsKey(),
963  tobasis1.GetPointsKey(),
964  &normal[i][0]);
965  Vmath::Vmul(nq_face,work,1,normal[i],1,normal[i],1);
966  }
967 
968  // Normalise to obtain unit normals.
969  Vmath::Zero(nq_face,work,1);
970  for(i = 0; i < GetCoordim(); ++i)
971  {
972  Vmath::Vvtvp(nq_face,normal[i],1,normal[i],1,work,1,work,1);
973  }
974 
975  Vmath::Vsqrt(nq_face,work,1,work,1);
976  Vmath::Sdiv (nq_face,1.0,work,1,work,1);
977 
978  for(i = 0; i < GetCoordim(); ++i)
979  {
980  Vmath::Vmul(nq_face,normal[i],1,work,1,normal[i],1);
981  }
982  }
983  }
984 
985  void PrismExp::v_MassMatrixOp(
986  const Array<OneD, const NekDouble> &inarray,
987  Array<OneD, NekDouble> &outarray,
988  const StdRegions::StdMatrixKey &mkey)
989  {
990  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
991  }
992 
993  void PrismExp::v_LaplacianMatrixOp(
994  const Array<OneD, const NekDouble> &inarray,
995  Array<OneD, NekDouble> &outarray,
996  const StdRegions::StdMatrixKey &mkey)
997  {
998  PrismExp::LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
999  }
1000 
1001  void PrismExp::v_LaplacianMatrixOp(
1002  const int k1,
1003  const int k2,
1004  const Array<OneD, const NekDouble> &inarray,
1005  Array<OneD, NekDouble> &outarray,
1006  const StdRegions::StdMatrixKey &mkey)
1007  {
1008  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,mkey);
1009  }
1010 
1011  void PrismExp::v_HelmholtzMatrixOp(
1012  const Array<OneD, const NekDouble> &inarray,
1013  Array<OneD, NekDouble> &outarray,
1014  const StdRegions::StdMatrixKey &mkey)
1015  {
1016  PrismExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1017  }
1018 
1019  void PrismExp::v_GeneralMatrixOp_MatOp(
1020  const Array<OneD, const NekDouble> &inarray,
1021  Array<OneD, NekDouble> &outarray,
1022  const StdRegions::StdMatrixKey &mkey)
1023  {
1024  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1025 
1026  if(inarray.get() == outarray.get())
1027  {
1028  Array<OneD,NekDouble> tmp(m_ncoeffs);
1029  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1030 
1031  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1032  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1033  }
1034  else
1035  {
1036  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1037  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1038  }
1039  }
1040 
1041  void PrismExp::v_SVVLaplacianFilter(
1042  Array<OneD, NekDouble> &array,
1043  const StdRegions::StdMatrixKey &mkey)
1044  {
1045  int nq = GetTotPoints();
1046 
1047  // Calculate sqrt of the Jacobian
1048  Array<OneD, const NekDouble> jac =
1049  m_metricinfo->GetJac(GetPointsKeys());
1050  Array<OneD, NekDouble> sqrt_jac(nq);
1051  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1052  {
1053  Vmath::Vsqrt(nq,jac,1,sqrt_jac,1);
1054  }
1055  else
1056  {
1057  Vmath::Fill(nq,sqrt(jac[0]),sqrt_jac,1);
1058  }
1059 
1060  // Multiply array by sqrt(Jac)
1061  Vmath::Vmul(nq,sqrt_jac,1,array,1,array,1);
1062 
1063  // Apply std region filter
1064  StdPrismExp::v_SVVLaplacianFilter( array, mkey);
1065 
1066  // Divide by sqrt(Jac)
1067  Vmath::Vdiv(nq,array,1,sqrt_jac,1,array,1);
1068  }
1069 
1070 
1071  //---------------------------------------
1072  // Matrix creation functions
1073  //---------------------------------------
1074 
1075  DNekMatSharedPtr PrismExp::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
1076  {
1077  DNekMatSharedPtr returnval;
1078 
1079  switch(mkey.GetMatrixType())
1080  {
1081  case StdRegions::eHybridDGHelmholtz:
1082  case StdRegions::eHybridDGLamToU:
1083  case StdRegions::eHybridDGLamToQ0:
1084  case StdRegions::eHybridDGLamToQ1:
1085  case StdRegions::eHybridDGLamToQ2:
1086  case StdRegions::eHybridDGHelmBndLam:
1087  case StdRegions::eInvLaplacianWithUnityMean:
1088  returnval = Expansion3D::v_GenMatrix(mkey);
1089  break;
1090  default:
1091  returnval = StdPrismExp::v_GenMatrix(mkey);
1092  break;
1093  }
1094 
1095  return returnval;
1096  }
1097 
1098  DNekMatSharedPtr PrismExp::v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
1099  {
1100  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1101  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1102  LibUtilities::BasisKey bkey2 = m_base[2]->GetBasisKey();
1103  StdRegions::StdPrismExpSharedPtr tmp =
1104  MemoryManager<StdPrismExp>::AllocateSharedPtr(bkey0, bkey1, bkey2);
1105 
1106  return tmp->GetStdMatrix(mkey);
1107  }
1108 
1109  DNekScalMatSharedPtr PrismExp::v_GetLocMatrix(const MatrixKey &mkey)
1110  {
1111  return m_matrixManager[mkey];
1112  }
1113 
1114  DNekScalBlkMatSharedPtr PrismExp::v_GetLocStaticCondMatrix(const MatrixKey &mkey)
1115  {
1116  return m_staticCondMatrixManager[mkey];
1117  }
1118 
1119  void PrismExp::v_DropLocStaticCondMatrix(const MatrixKey &mkey)
1120  {
1121  m_staticCondMatrixManager.DeleteObject(mkey);
1122  }
1123 
1124  DNekScalMatSharedPtr PrismExp::CreateMatrix(const MatrixKey &mkey)
1125  {
1126  DNekScalMatSharedPtr returnval;
1127  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1128 
1129  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,
1130  "Geometric information is not set up");
1131 
1132  switch(mkey.GetMatrixType())
1133  {
1134  case StdRegions::eMass:
1135  {
1136  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1137  {
1138  NekDouble one = 1.0;
1139  DNekMatSharedPtr mat = GenMatrix(mkey);
1140  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1141  }
1142  else
1143  {
1144  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1145  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1146  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1147  }
1148  break;
1149  }
1150 
1151  case StdRegions::eInvMass:
1152  {
1153  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1154  {
1155  NekDouble one = 1.0;
1156  StdRegions::StdMatrixKey masskey(StdRegions::eMass,DetShapeType(),*this);
1157  DNekMatSharedPtr mat = GenMatrix(masskey);
1158  mat->Invert();
1159 
1160  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1161  }
1162  else
1163  {
1164  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1165  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1166  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(fac,mat);
1167  }
1168  break;
1169  }
1170 
1171  case StdRegions::eWeakDeriv0:
1172  case StdRegions::eWeakDeriv1:
1173  case StdRegions::eWeakDeriv2:
1174  {
1175  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1176  {
1177  NekDouble one = 1.0;
1178  DNekMatSharedPtr mat = GenMatrix(mkey);
1179 
1180  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1181  }
1182  else
1183  {
1184  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1185  Array<TwoD, const NekDouble> df =
1186  m_metricinfo->GetDerivFactors(ptsKeys);
1187  int dir = 0;
1188 
1189  switch(mkey.GetMatrixType())
1190  {
1191  case StdRegions::eWeakDeriv0:
1192  dir = 0;
1193  break;
1194  case StdRegions::eWeakDeriv1:
1195  dir = 1;
1196  break;
1197  case StdRegions::eWeakDeriv2:
1198  dir = 2;
1199  break;
1200  default:
1201  break;
1202  }
1203 
1204  MatrixKey deriv0key(StdRegions::eWeakDeriv0,
1205  mkey.GetShapeType(), *this);
1206  MatrixKey deriv1key(StdRegions::eWeakDeriv1,
1207  mkey.GetShapeType(), *this);
1208  MatrixKey deriv2key(StdRegions::eWeakDeriv2,
1209  mkey.GetShapeType(), *this);
1210 
1211  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
1212  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
1213  DNekMat &deriv2 = *GetStdMatrix(deriv2key);
1214 
1215  int rows = deriv0.GetRows();
1216  int cols = deriv1.GetColumns();
1217 
1218  DNekMatSharedPtr WeakDeriv = MemoryManager<DNekMat>
1219  ::AllocateSharedPtr(rows,cols);
1220 
1221  (*WeakDeriv) = df[3*dir ][0]*deriv0
1222  + df[3*dir+1][0]*deriv1
1223  + df[3*dir+2][0]*deriv2;
1224 
1225  returnval = MemoryManager<DNekScalMat>
1226  ::AllocateSharedPtr(jac,WeakDeriv);
1227  }
1228  break;
1229  }
1230 
1231  case StdRegions::eLaplacian:
1232  {
1233  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1234  mkey.GetNVarCoeff() > 0 ||
1235  mkey.ConstFactorExists(
1236  StdRegions::eFactorSVVCutoffRatio))
1237  {
1238  NekDouble one = 1.0;
1239  DNekMatSharedPtr mat = GenMatrix(mkey);
1240  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1241  }
1242  else
1243  {
1244  MatrixKey lap00key(StdRegions::eLaplacian00,
1245  mkey.GetShapeType(), *this);
1246  MatrixKey lap01key(StdRegions::eLaplacian01,
1247  mkey.GetShapeType(), *this);
1248  MatrixKey lap02key(StdRegions::eLaplacian02,
1249  mkey.GetShapeType(), *this);
1250  MatrixKey lap11key(StdRegions::eLaplacian11,
1251  mkey.GetShapeType(), *this);
1252  MatrixKey lap12key(StdRegions::eLaplacian12,
1253  mkey.GetShapeType(), *this);
1254  MatrixKey lap22key(StdRegions::eLaplacian22,
1255  mkey.GetShapeType(), *this);
1256 
1257  DNekMat &lap00 = *GetStdMatrix(lap00key);
1258  DNekMat &lap01 = *GetStdMatrix(lap01key);
1259  DNekMat &lap02 = *GetStdMatrix(lap02key);
1260  DNekMat &lap11 = *GetStdMatrix(lap11key);
1261  DNekMat &lap12 = *GetStdMatrix(lap12key);
1262  DNekMat &lap22 = *GetStdMatrix(lap22key);
1263 
1264  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1265  Array<TwoD, const NekDouble> gmat
1266  = m_metricinfo->GetGmat(ptsKeys);
1267 
1268  int rows = lap00.GetRows();
1269  int cols = lap00.GetColumns();
1270 
1271  DNekMatSharedPtr lap = MemoryManager<DNekMat>
1272  ::AllocateSharedPtr(rows,cols);
1273 
1274  (*lap) = gmat[0][0]*lap00
1275  + gmat[4][0]*lap11
1276  + gmat[8][0]*lap22
1277  + gmat[3][0]*(lap01 + Transpose(lap01))
1278  + gmat[6][0]*(lap02 + Transpose(lap02))
1279  + gmat[7][0]*(lap12 + Transpose(lap12));
1280 
1281  returnval = MemoryManager<DNekScalMat>
1282  ::AllocateSharedPtr(jac,lap);
1283  }
1284  break;
1285  }
1286 
1287  case StdRegions::eHelmholtz:
1288  {
1289  NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorLambda);
1290  MatrixKey masskey(StdRegions::eMass,
1291  mkey.GetShapeType(), *this);
1292  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1293  MatrixKey lapkey(StdRegions::eLaplacian,
1294  mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1295  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1296 
1297  int rows = LapMat.GetRows();
1298  int cols = LapMat.GetColumns();
1299 
1300  DNekMatSharedPtr helm = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1301 
1302  NekDouble one = 1.0;
1303  (*helm) = LapMat + factor*MassMat;
1304 
1305  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
1306  break;
1307  }
1308  case StdRegions::eHybridDGHelmholtz:
1309  case StdRegions::eHybridDGLamToU:
1310  case StdRegions::eHybridDGLamToQ0:
1311  case StdRegions::eHybridDGLamToQ1:
1312  case StdRegions::eHybridDGHelmBndLam:
1313  case StdRegions::eInvLaplacianWithUnityMean:
1314  {
1315  NekDouble one = 1.0;
1316 
1317  DNekMatSharedPtr mat = GenMatrix(mkey);
1318  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1319 
1320  break;
1321  }
1322 
1323  case StdRegions::eInvHybridDGHelmholtz:
1324  {
1325  NekDouble one = 1.0;
1326 
1327  MatrixKey hkey(StdRegions::eHybridDGHelmholtz, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1328 // StdRegions::StdMatrixKey hkey(StdRegions::eHybridDGHelmholtz,
1329 // DetExpansionType(),*this,
1330 // mkey.GetConstant(0),
1331 // mkey.GetConstant(1));
1332  DNekMatSharedPtr mat = GenMatrix(hkey);
1333 
1334  mat->Invert();
1335  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1336  break;
1337  }
1338 
1339  case StdRegions::eIProductWRTBase:
1340  {
1341  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1342  {
1343  NekDouble one = 1.0;
1344  DNekMatSharedPtr mat = GenMatrix(mkey);
1345  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1346  }
1347  else
1348  {
1349  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1350  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1351  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1352  }
1353  break;
1354  }
1355  case StdRegions::ePreconLinearSpace:
1356  {
1357  NekDouble one = 1.0;
1358  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1359  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1360  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1361  DNekMatSharedPtr R=BuildVertexMatrix(A);
1362 
1363  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1364  }
1365  break;
1366  case StdRegions::ePreconLinearSpaceMass:
1367  {
1368  NekDouble one = 1.0;
1369  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1370  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1371  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1372  DNekMatSharedPtr R=BuildVertexMatrix(A);
1373 
1374  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1375  }
1376  break;
1377  case StdRegions::ePreconR:
1378  {
1379  NekDouble one = 1.0;
1380  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this,mkey.GetConstFactors(), mkey.GetVarCoeffs());
1381  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1382  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1383 
1384  DNekScalMatSharedPtr Atmp;
1385  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1386 
1387  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1388  }
1389  break;
1390  case StdRegions::ePreconRMass:
1391  {
1392  NekDouble one = 1.0;
1393  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1394  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1395  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1396 
1397  DNekScalMatSharedPtr Atmp;
1398  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1399 
1400  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1401  }
1402  break;
1403  default:
1404  {
1405  NekDouble one = 1.0;
1406  DNekMatSharedPtr mat = GenMatrix(mkey);
1407 
1408  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1409  }
1410  }
1411 
1412  return returnval;
1413  }
1414 
1415  DNekScalBlkMatSharedPtr PrismExp::CreateStaticCondMatrix(const MatrixKey &mkey)
1416  {
1417  DNekScalBlkMatSharedPtr returnval;
1418 
1419  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
1420 
1421  // set up block matrix system
1422  unsigned int nbdry = NumBndryCoeffs();
1423  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
1424  unsigned int exp_size[] = {nbdry, nint};
1425  unsigned int nblks=2;
1426  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks, nblks, exp_size, exp_size);
1427  NekDouble factor = 1.0;
1428 
1429  switch(mkey.GetMatrixType())
1430  {
1431  case StdRegions::eLaplacian:
1432  case StdRegions::eHelmholtz: // special case since Helmholtz not defined in StdRegions
1433  // use Deformed case for both regular and deformed geometries
1434  factor = 1.0;
1435  goto UseLocRegionsMatrix;
1436  break;
1437  default:
1438  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1439  {
1440  factor = 1.0;
1441  goto UseLocRegionsMatrix;
1442  }
1443  else
1444  {
1445  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1446  factor = mat->Scale();
1447  goto UseStdRegionsMatrix;
1448  }
1449  break;
1450  UseStdRegionsMatrix:
1451  {
1452  NekDouble invfactor = 1.0/factor;
1453  NekDouble one = 1.0;
1454  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
1455  DNekScalMatSharedPtr Atmp;
1456  DNekMatSharedPtr Asubmat;
1457 
1458  //TODO: check below
1459  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
1460  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
1461  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
1462  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
1463  }
1464  break;
1465  UseLocRegionsMatrix:
1466  {
1467  int i,j;
1468  NekDouble invfactor = 1.0/factor;
1469  NekDouble one = 1.0;
1470  DNekScalMat &mat = *GetLocMatrix(mkey);
1471  DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
1472  DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
1473  DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
1474  DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
1475 
1476  Array<OneD,unsigned int> bmap(nbdry);
1477  Array<OneD,unsigned int> imap(nint);
1478  GetBoundaryMap(bmap);
1479  GetInteriorMap(imap);
1480 
1481  for(i = 0; i < nbdry; ++i)
1482  {
1483  for(j = 0; j < nbdry; ++j)
1484  {
1485  (*A)(i,j) = mat(bmap[i],bmap[j]);
1486  }
1487 
1488  for(j = 0; j < nint; ++j)
1489  {
1490  (*B)(i,j) = mat(bmap[i],imap[j]);
1491  }
1492  }
1493 
1494  for(i = 0; i < nint; ++i)
1495  {
1496  for(j = 0; j < nbdry; ++j)
1497  {
1498  (*C)(i,j) = mat(imap[i],bmap[j]);
1499  }
1500 
1501  for(j = 0; j < nint; ++j)
1502  {
1503  (*D)(i,j) = mat(imap[i],imap[j]);
1504  }
1505  }
1506 
1507  // Calculate static condensed system
1508  if(nint)
1509  {
1510  D->Invert();
1511  (*B) = (*B)*(*D);
1512  (*A) = (*A) - (*B)*(*C);
1513  }
1514 
1515  DNekScalMatSharedPtr Atmp;
1516 
1517  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
1518  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
1519  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
1520  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
1521 
1522  }
1523  break;
1524  }
1525  return returnval;
1526  }
1527 
1528 
1529  /**
1530  * @brief Calculate the Laplacian multiplication in a matrix-free
1531  * manner.
1532  *
1533  * This function is the kernel of the Laplacian matrix-free operator,
1534  * and is used in #v_HelmholtzMatrixOp_MatFree to determine the effect
1535  * of the Helmholtz operator in a similar fashion.
1536  *
1537  * The majority of the calculation is precisely the same as in the
1538  * hexahedral expansion; however the collapsed co-ordinate system must
1539  * be taken into account when constructing the geometric factors. How
1540  * this is done is detailed more exactly in the tetrahedral expansion.
1541  * On entry to this function, the input #inarray must be in its
1542  * backwards-transformed state (i.e. \f$\mathbf{u} =
1543  * \mathbf{B}\hat{\mathbf{u}}\f$). The output is in coefficient space.
1544  *
1545  * @see %TetExp::v_HelmholtzMatrixOp_MatFree
1546  */
1547  void PrismExp::v_LaplacianMatrixOp_MatFree_Kernel(
1548  const Array<OneD, const NekDouble> &inarray,
1549  Array<OneD, NekDouble> &outarray,
1550  Array<OneD, NekDouble> &wsp)
1551  {
1552  int nquad0 = m_base[0]->GetNumPoints();
1553  int nquad1 = m_base[1]->GetNumPoints();
1554  int nquad2 = m_base[2]->GetNumPoints();
1555  int nqtot = nquad0*nquad1*nquad2;
1556  int i;
1557 
1558  // Set up temporary storage.
1559  Array<OneD,NekDouble> alloc(11*nqtot,0.0);
1560  Array<OneD,NekDouble> wsp1 (alloc ); // TensorDeriv 1
1561  Array<OneD,NekDouble> wsp2 (alloc+ 1*nqtot); // TensorDeriv 2
1562  Array<OneD,NekDouble> wsp3 (alloc+ 2*nqtot); // TensorDeriv 3
1563  Array<OneD,NekDouble> g0 (alloc+ 3*nqtot); // g0
1564  Array<OneD,NekDouble> g1 (alloc+ 4*nqtot); // g1
1565  Array<OneD,NekDouble> g2 (alloc+ 5*nqtot); // g2
1566  Array<OneD,NekDouble> g3 (alloc+ 6*nqtot); // g3
1567  Array<OneD,NekDouble> g4 (alloc+ 7*nqtot); // g4
1568  Array<OneD,NekDouble> g5 (alloc+ 8*nqtot); // g5
1569  Array<OneD,NekDouble> h0 (alloc+ 3*nqtot); // h0 == g0
1570  Array<OneD,NekDouble> h1 (alloc+ 6*nqtot); // h1 == g3
1571  Array<OneD,NekDouble> wsp4 (alloc+ 4*nqtot); // wsp4 == g1
1572  Array<OneD,NekDouble> wsp5 (alloc+ 5*nqtot); // wsp5 == g2
1573  Array<OneD,NekDouble> wsp6 (alloc+ 8*nqtot); // wsp6 == g5
1574  Array<OneD,NekDouble> wsp7 (alloc+ 3*nqtot); // wsp7 == g0
1575  Array<OneD,NekDouble> wsp8 (alloc+ 9*nqtot); // wsp8
1576  Array<OneD,NekDouble> wsp9 (alloc+10*nqtot); // wsp9
1577 
1578  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
1579  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
1580  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata();
1581  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
1582  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
1583  const Array<OneD, const NekDouble>& dbase2 = m_base[2]->GetDbdata();
1584 
1585  // Step 1. LAPLACIAN MATRIX OPERATION
1586  // wsp1 = du_dxi1 = D_xi1 * wsp0 = D_xi1 * u
1587  // wsp2 = du_dxi2 = D_xi2 * wsp0 = D_xi2 * u
1588  // wsp3 = du_dxi3 = D_xi3 * wsp0 = D_xi3 * u
1589  StdExpansion3D::PhysTensorDeriv(inarray,wsp1,wsp2,wsp3);
1590 
1591  const Array<TwoD, const NekDouble>& df =
1592  m_metricinfo->GetDerivFactors(GetPointsKeys());
1593  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
1594  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1595 
1596  // Step 2. Calculate the metric terms of the collapsed
1597  // coordinate transformation (Spencer's book P152)
1598  for (i = 0; i < nquad2; ++i)
1599  {
1600  Vmath::Fill(nquad0*nquad1, 2.0/(1.0-z2[i]), &h0[0]+i*nquad0*nquad1,1);
1601  Vmath::Fill(nquad0*nquad1, 2.0/(1.0-z2[i]), &h1[0]+i*nquad0*nquad1,1);
1602  }
1603  for (i = 0; i < nquad0; i++)
1604  {
1605  Blas::Dscal(nquad1*nquad2, 0.5*(1+z0[i]), &h1[0]+i, nquad0);
1606  }
1607 
1608  // Step 3. Construct combined metric terms for physical space to
1609  // collapsed coordinate system. Order of construction optimised
1610  // to minimise temporary storage
1611  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1612  {
1613  // wsp4 = d eta_1/d x_1
1614  Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &df[2][0], 1, &h1[0], 1, &wsp4[0], 1);
1615  // wsp5 = d eta_2/d x_1
1616  Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &df[5][0], 1, &h1[0], 1, &wsp5[0], 1);
1617  // wsp6 = d eta_3/d x_1d
1618  Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &df[8][0], 1, &h1[0], 1, &wsp6[0], 1);
1619 
1620  // g0 (overwrites h0)
1621  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0], 1, &g0[0], 1);
1622  Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1623 
1624  // g3 (overwrites h1)
1625  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &wsp4[0], 1, &df[4][0], 1, &wsp5[0], 1, &g3[0], 1);
1626  Vmath::Vvtvp (nqtot, &df[7][0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1627 
1628  // g4
1629  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0], 1, &g4[0], 1);
1630  Vmath::Vvtvp (nqtot, &df[8][0], 1, &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1631 
1632  // Overwrite wsp4/5/6 with g1/2/5
1633  // g1
1634  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &df[1][0], 1, &df[4][0], 1, &df[4][0], 1, &g1[0], 1);
1635  Vmath::Vvtvp (nqtot, &df[7][0], 1, &df[7][0], 1, &g1[0], 1, &g1[0], 1);
1636 
1637  // g2
1638  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &df[2][0], 1, &df[5][0], 1, &df[5][0], 1, &g2[0], 1);
1639  Vmath::Vvtvp (nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1640 
1641  // g5
1642  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &df[2][0], 1, &df[4][0], 1, &df[5][0], 1, &g5[0], 1);
1643  Vmath::Vvtvp (nqtot, &df[7][0], 1, &df[8][0], 1, &g5[0], 1, &g5[0], 1);
1644  }
1645  else
1646  {
1647  // wsp4 = d eta_1/d x_1
1648  Vmath::Svtsvtp(nqtot, df[0][0], &h0[0], 1, df[2][0], &h1[0], 1, &wsp4[0], 1);
1649  // wsp5 = d eta_2/d x_1
1650  Vmath::Svtsvtp(nqtot, df[3][0], &h0[0], 1, df[5][0], &h1[0], 1, &wsp5[0], 1);
1651  // wsp6 = d eta_3/d x_1
1652  Vmath::Svtsvtp(nqtot, df[6][0], &h0[0], 1, df[8][0], &h1[0], 1, &wsp6[0], 1);
1653 
1654  // g0 (overwrites h0)
1655  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0], 1, &g0[0], 1);
1656  Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1657 
1658  // g3 (overwrites h1)
1659  Vmath::Svtsvtp(nqtot, df[1][0], &wsp4[0], 1, df[4][0], &wsp5[0], 1, &g3[0], 1);
1660  Vmath::Svtvp (nqtot, df[7][0], &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1661 
1662  // g4
1663  Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1, &g4[0], 1);
1664  Vmath::Svtvp (nqtot, df[8][0], &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1665 
1666  // Overwrite wsp4/5/6 with g1/2/5
1667  // g1
1668  Vmath::Fill(nqtot, df[1][0]*df[1][0] + df[4][0]*df[4][0] + df[7][0]*df[7][0], &g1[0], 1);
1669 
1670  // g2
1671  Vmath::Fill(nqtot, df[2][0]*df[2][0] + df[5][0]*df[5][0] + df[8][0]*df[8][0], &g2[0], 1);
1672 
1673  // g5
1674  Vmath::Fill(nqtot, df[1][0]*df[2][0] + df[4][0]*df[5][0] + df[7][0]*df[8][0], &g5[0], 1);
1675  }
1676  // Compute component derivatives into wsp7, 8, 9 (wsp7 overwrites
1677  // g0).
1678  Vmath::Vvtvvtp(nqtot,&g0[0],1,&wsp1[0],1,&g3[0],1,&wsp2[0],1,&wsp7[0],1);
1679  Vmath::Vvtvp (nqtot,&g4[0],1,&wsp3[0],1,&wsp7[0],1,&wsp7[0],1);
1680  Vmath::Vvtvvtp(nqtot,&g1[0],1,&wsp2[0],1,&g3[0],1,&wsp1[0],1,&wsp8[0],1);
1681  Vmath::Vvtvp (nqtot,&g5[0],1,&wsp3[0],1,&wsp8[0],1,&wsp8[0],1);
1682  Vmath::Vvtvvtp(nqtot,&g2[0],1,&wsp3[0],1,&g4[0],1,&wsp1[0],1,&wsp9[0],1);
1683  Vmath::Vvtvp (nqtot,&g5[0],1,&wsp2[0],1,&wsp9[0],1,&wsp9[0],1);
1684 
1685  // Step 4.
1686  // Multiply by quadrature metric
1687  MultiplyByQuadratureMetric(wsp7,wsp7);
1688  MultiplyByQuadratureMetric(wsp8,wsp8);
1689  MultiplyByQuadratureMetric(wsp9,wsp9);
1690 
1691  // Perform inner product w.r.t derivative bases.
1692  IProductWRTBase_SumFacKernel(dbase0,base1,base2,wsp7,wsp1, wsp,false,true,true);
1693  IProductWRTBase_SumFacKernel(base0,dbase1,base2,wsp8,wsp2, wsp,true,false,true);
1694  IProductWRTBase_SumFacKernel(base0,base1,dbase2,wsp9,outarray,wsp,true,true,false);
1695 
1696  // Step 5.
1697  // Sum contributions from wsp1, wsp2 and outarray.
1698  Vmath::Vadd(m_ncoeffs,wsp1.get(),1,outarray.get(),1,outarray.get(),1);
1699  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
1700  }
1701 
1702 
1703  void PrismExp::v_GetSimplexEquiSpacedConnectivity(
1704  Array<OneD, int> &conn,
1705  bool oldstandard)
1706  {
1707  boost::ignore_unused(oldstandard);
1708 
1709  int np0 = m_base[0]->GetNumPoints();
1710  int np1 = m_base[1]->GetNumPoints();
1711  int np2 = m_base[2]->GetNumPoints();
1712  int np = max(np0,max(np1,np2));
1713  Array<OneD, int> prismpt(6);
1714  bool standard = true;
1715 
1716  int vid0 = m_geom->GetVid(0);
1717  int vid1 = m_geom->GetVid(1);
1718  int vid2 = m_geom->GetVid(4);
1719  int rotate = 0;
1720 
1721  // sort out prism rotation according to
1722  if((vid2 < vid1)&&(vid2 < vid0)) // top triangle vertex is lowest id
1723  {
1724  rotate = 0;
1725  if(vid0 > vid1)
1726  {
1727  standard = false;// reverse base direction
1728  }
1729  }
1730  else if((vid1 < vid2)&&(vid1 < vid0))
1731  {
1732  rotate = 1;
1733  if(vid2 > vid0)
1734  {
1735  standard = false;// reverse base direction
1736  }
1737  }
1738  else if ((vid0 < vid2)&&(vid0 < vid1))
1739  {
1740  rotate = 2;
1741  if(vid1 > vid2)
1742  {
1743  standard = false; // reverse base direction
1744  }
1745  }
1746 
1747  conn = Array<OneD, int>(12*(np-1)*(np-1)*(np-1));
1748 
1749  int row = 0;
1750  int rowp1 = 0;
1751  int plane = 0;
1752  int row1 = 0;
1753  int row1p1 = 0;
1754  int planep1 = 0;
1755  int cnt = 0;
1756 
1757 
1758  Array<OneD, int> rot(3);
1759 
1760  rot[0] = (0+rotate)%3;
1761  rot[1] = (1+rotate)%3;
1762  rot[2] = (2+rotate)%3;
1763 
1764  // lower diagonal along 1-3 on base
1765  for(int i = 0; i < np-1; ++i)
1766  {
1767  planep1 += (np-i)*np;
1768  row = 0; // current plane row offset
1769  rowp1 = 0; // current plane row plus one offset
1770  row1 = 0; // next plane row offset
1771  row1p1 = 0; // nex plane row plus one offset
1772  if(standard == false)
1773  {
1774  for(int j = 0; j < np-1; ++j)
1775  {
1776  rowp1 += np-i;
1777  row1p1 += np-i-1;
1778  for(int k = 0; k < np-i-2; ++k)
1779  {
1780  // bottom prism block
1781  prismpt[rot[0]] = plane + row + k;
1782  prismpt[rot[1]] = plane + row + k+1;
1783  prismpt[rot[2]] = planep1 + row1 + k;
1784 
1785  prismpt[3+rot[0]] = plane + rowp1 + k;
1786  prismpt[3+rot[1]] = plane + rowp1 + k+1;
1787  prismpt[3+rot[2]] = planep1 + row1p1 + k;
1788 
1789  conn[cnt++] = prismpt[0];
1790  conn[cnt++] = prismpt[1];
1791  conn[cnt++] = prismpt[3];
1792  conn[cnt++] = prismpt[2];
1793 
1794  conn[cnt++] = prismpt[5];
1795  conn[cnt++] = prismpt[2];
1796  conn[cnt++] = prismpt[3];
1797  conn[cnt++] = prismpt[4];
1798 
1799  conn[cnt++] = prismpt[3];
1800  conn[cnt++] = prismpt[1];
1801  conn[cnt++] = prismpt[4];
1802  conn[cnt++] = prismpt[2];
1803 
1804  // upper prism block.
1805  prismpt[rot[0]] = planep1 + row1 + k+1;
1806  prismpt[rot[1]] = planep1 + row1 + k;
1807  prismpt[rot[2]] = plane + row + k+1;
1808 
1809  prismpt[3+rot[0]] = planep1 + row1p1 + k+1;
1810  prismpt[3+rot[1]] = planep1 + row1p1 + k;
1811  prismpt[3+rot[2]] = plane + rowp1 + k+1;
1812 
1813 
1814  conn[cnt++] = prismpt[0];
1815  conn[cnt++] = prismpt[1];
1816  conn[cnt++] = prismpt[2];
1817  conn[cnt++] = prismpt[5];
1818 
1819  conn[cnt++] = prismpt[5];
1820  conn[cnt++] = prismpt[0];
1821  conn[cnt++] = prismpt[4];
1822  conn[cnt++] = prismpt[1];
1823 
1824  conn[cnt++] = prismpt[3];
1825  conn[cnt++] = prismpt[4];
1826  conn[cnt++] = prismpt[0];
1827  conn[cnt++] = prismpt[5];
1828 
1829  }
1830 
1831  // bottom prism block
1832  prismpt[rot[0]] = plane + row + np-i-2;
1833  prismpt[rot[1]] = plane + row + np-i-1;
1834  prismpt[rot[2]] = planep1 + row1 + np-i-2;
1835 
1836  prismpt[3+rot[0]] = plane + rowp1 + np-i-2;
1837  prismpt[3+rot[1]] = plane + rowp1 + np-i-1;
1838  prismpt[3+rot[2]] = planep1 + row1p1 + np-i-2;
1839 
1840  conn[cnt++] = prismpt[0];
1841  conn[cnt++] = prismpt[1];
1842  conn[cnt++] = prismpt[3];
1843  conn[cnt++] = prismpt[2];
1844 
1845  conn[cnt++] = prismpt[5];
1846  conn[cnt++] = prismpt[2];
1847  conn[cnt++] = prismpt[3];
1848  conn[cnt++] = prismpt[4];
1849 
1850  conn[cnt++] = prismpt[3];
1851  conn[cnt++] = prismpt[1];
1852  conn[cnt++] = prismpt[4];
1853  conn[cnt++] = prismpt[2];
1854 
1855  row += np-i;
1856  row1 += np-i-1;
1857  }
1858 
1859  }
1860  else
1861  { // lower diagonal along 0-4 on base
1862  for(int j = 0; j < np-1; ++j)
1863  {
1864  rowp1 += np-i;
1865  row1p1 += np-i-1;
1866  for(int k = 0; k < np-i-2; ++k)
1867  {
1868  // bottom prism block
1869  prismpt[rot[0]] = plane + row + k;
1870  prismpt[rot[1]] = plane + row + k+1;
1871  prismpt[rot[2]] = planep1 + row1 + k;
1872 
1873  prismpt[3+rot[0]] = plane + rowp1 + k;
1874  prismpt[3+rot[1]] = plane + rowp1 + k+1;
1875  prismpt[3+rot[2]] = planep1 + row1p1 + k;
1876 
1877  conn[cnt++] = prismpt[0];
1878  conn[cnt++] = prismpt[1];
1879  conn[cnt++] = prismpt[4];
1880  conn[cnt++] = prismpt[2];
1881 
1882  conn[cnt++] = prismpt[4];
1883  conn[cnt++] = prismpt[3];
1884  conn[cnt++] = prismpt[0];
1885  conn[cnt++] = prismpt[2];
1886 
1887  conn[cnt++] = prismpt[3];
1888  conn[cnt++] = prismpt[4];
1889  conn[cnt++] = prismpt[5];
1890  conn[cnt++] = prismpt[2];
1891 
1892  // upper prism block.
1893  prismpt[rot[0]] = planep1 + row1 + k+1;
1894  prismpt[rot[1]] = planep1 + row1 + k;
1895  prismpt[rot[2]] = plane + row + k+1;
1896 
1897  prismpt[3+rot[0]] = planep1 + row1p1 + k+1;
1898  prismpt[3+rot[1]] = planep1 + row1p1 + k;
1899  prismpt[3+rot[2]] = plane + rowp1 + k+1;
1900 
1901  conn[cnt++] = prismpt[0];
1902  conn[cnt++] = prismpt[2];
1903  conn[cnt++] = prismpt[1];
1904  conn[cnt++] = prismpt[5];
1905 
1906  conn[cnt++] = prismpt[3];
1907  conn[cnt++] = prismpt[5];
1908  conn[cnt++] = prismpt[0];
1909  conn[cnt++] = prismpt[1];
1910 
1911  conn[cnt++] = prismpt[5];
1912  conn[cnt++] = prismpt[3];
1913  conn[cnt++] = prismpt[4];
1914  conn[cnt++] = prismpt[1];
1915  }
1916 
1917  // bottom prism block
1918  prismpt[rot[0]] = plane + row + np-i-2;
1919  prismpt[rot[1]] = plane + row + np-i-1;
1920  prismpt[rot[2]] = planep1 + row1 + np-i-2;
1921 
1922  prismpt[3+rot[0]] = plane + rowp1 + np-i-2;
1923  prismpt[3+rot[1]] = plane + rowp1 + np-i-1;
1924  prismpt[3+rot[2]] = planep1 + row1p1 + np-i-2;
1925 
1926  conn[cnt++] = prismpt[0];
1927  conn[cnt++] = prismpt[1];
1928  conn[cnt++] = prismpt[4];
1929  conn[cnt++] = prismpt[2];
1930 
1931  conn[cnt++] = prismpt[4];
1932  conn[cnt++] = prismpt[3];
1933  conn[cnt++] = prismpt[0];
1934  conn[cnt++] = prismpt[2];
1935 
1936  conn[cnt++] = prismpt[3];
1937  conn[cnt++] = prismpt[4];
1938  conn[cnt++] = prismpt[5];
1939  conn[cnt++] = prismpt[2];
1940 
1941  row += np-i;
1942  row1 += np-i-1;
1943  }
1944 
1945  }
1946  plane += (np-i)*np;
1947  }
1948  }
1949 
1950  }//end of namespace
1951 }//end of namespace
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:1011
Nektar::LibUtilities::BasisKey::GetNumPoints
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:133
Nektar::StdRegions::StdMatrixKey::GetVarCoeffs
const VarCoeffMap & GetVarCoeffs() const
Definition: StdMatrixKey.h:161
Nektar::StdRegions::StdExpansion::GenMatrix
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:1060
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::LocalRegions::PrismExp::m_matrixManager
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: PrismExp.h:210
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:476
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:494
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::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
Nektar::LocalRegions::Expansion::BuildTransformationMatrix
DNekMatSharedPtr BuildTransformationMatrix(const DNekScalMatSharedPtr &r_bnd, const StdRegions::MatrixType matrixType)
Definition: Expansion.cpp:90
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:469
Nektar::StdRegions::StdExpansion::GetNumPoints
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:228
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:411
Nektar::DNekScalMatSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:69
Nektar::StdRegions::StdExpansion::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:945
Nektar::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::Expansion::GetGeom
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:167
Nektar::SpatialDomains::GeomFactorsSharedPtr
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:83
Nektar::DNekScalBlkMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:73
Nektar::LocalRegions::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:1019
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
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:537
Nektar::StdRegions::StdExpansion::LaplacianMatrixOp_MatFree
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:1508
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:488
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:445
Nektar::LibUtilities::eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
Nektar::LocalRegions::Expansion::m_metricinfo
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:128
std
STL namespace.
Vmath::Sdiv
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:274
Nektar::LocalRegions::PrismExp::v_GetSimplexEquiSpacedConnectivity
virtual void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true)
Definition: PrismExp.cpp:1703
Nektar::LocalRegions::Expansion::BuildVertexMatrix
DNekMatSharedPtr BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
Definition: Expansion.cpp:98
Nektar::LocalRegions::PrismExp::v_GetLocStaticCondMatrix
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: PrismExp.cpp:1114
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:985
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:244
Nektar::LocalRegions::Expansion::m_geom
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:127
Nektar::StdRegions::StdMatrixKey::GetMatrixType
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:81
Nektar::StdRegions::StdExpansion::GetLocStaticCondMatrix
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
Definition: StdExpansion.h:761
Nektar::LibUtilities::BasisKey::GetPointsKey
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:150
Nektar::LocalRegions::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_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::DNekBlkMatSharedPtr
std::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:71
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:714
Nektar::LocalRegions::PrismExp::v_GetStdExp
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
Definition: PrismExp.cpp:451
Nektar::LocalRegions::PrismExp::v_GetFacePhysMap
virtual void v_GetFacePhysMap(const int face, Array< OneD, int > &outarray)
Definition: PrismExp.cpp:600
Nektar::LocalRegions::PrismExp::CreateMatrix
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition: PrismExp.cpp:1124
Nektar::StdRegions::StdExpansion::GetPointsKeys
const LibUtilities::PointsKeyVector GetPointsKeys() const
Definition: StdExpansion.h:1330
SegGeom.h
Nektar::LibUtilities::Interp2D
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis...
Definition: Interp.cpp:115
Nektar::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::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: PrismExp.cpp:1098
Nektar::StdRegions::StdExpansionSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1926
Nektar::StdRegions::StdExpansion::GetStdStaticCondMatrix
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:719
Nektar::LocalRegions::PrismExp::v_StdPhysEvaluate
virtual NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
Definition: PrismExp.cpp:507
Nektar::StdRegions::StdExpansion3D::IProductWRTBase_SumFacKernel
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
Definition: StdExpansion3D.cpp:128
Nektar::LocalRegions::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:1547
Nektar::LocalRegions::PrismExp::m_staticCondMatrixManager
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: PrismExp.h:212
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
Nektar::LocalRegions::PrismExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: PrismExp.cpp:1075
PrismExp.h
Nektar::LibUtilities::eModified_B
Principle Modified Functions .
Definition: BasisType.h: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::Expansion::v_GetCoords
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: Expansion.cpp:231
Nektar::LocalRegions::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_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
Definition: PrismExp.cpp:1041
Nektar::LocalRegions::PrismExp::v_DropLocStaticCondMatrix
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: PrismExp.cpp:1119
Nektar::LocalRegions::PrismExp::v_ComputeFaceNormal
void v_ComputeFaceNormal(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:704
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:161
Nektar::StdRegions::StdExpansion::GetInteriorMap
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:817
Nektar::StdRegions::StdPrismExpSharedPtr
std::shared_ptr< StdPrismExp > StdPrismExpSharedPtr
Definition: StdPrismExp.h:257
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:217
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:59
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Blas::Dgemv
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = A x where A[m x n].
Definition: Blas.hpp:168
Interp.h
Nektar::StdRegions::StdExpansion3D::m_faceNormals
std::map< int, NormalVector > m_faceNormals
Definition: StdExpansion3D.h:203
Nektar::LocalRegions::PrismExp::v_GetLocMatrix
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: PrismExp.cpp:1109
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:121
Nektar::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:130
Nektar::LocalRegions::Expansion::GetLocMatrix
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:85
Blas::Dscal
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:125
Nektar::StdRegions::StdMatrixKey::GetConstFactors
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:135
Nektar::StdRegions::StdExpansion::GetPointsType
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:215
Nektar::LocalRegions::Expansion3D::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: Expansion3D.cpp:473
Vmath::Vvtvvtp
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:540
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:130
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:274
Nektar::SpatialDomains::eRegular
Geometry is straight-sided with constant geometric factors.
Definition: SpatialDomains.hpp:54
Nektar::StdRegions::StdExpansion::GetBasisType
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:164
Nektar::StdRegions::StdExpansion::DetFaceBasisKey
const LibUtilities::BasisKey DetFaceBasisKey(const int i, const int k) const
Definition: StdExpansion.h:323
Nektar::LocalRegions::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:515
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_GetLinStdExp
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
Definition: PrismExp.cpp:459
Vmath::Svtsvtp
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
vvtvvtp (scalar times vector plus scalar times vector):
Definition: Vmath.cpp:594
Nektar::StdRegions::StdExpansion::StdExpansion
StdExpansion()
Default Constructor.
Definition: StdExpansion.cpp:45
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:140
Nektar::SpatialDomains::PrismGeomSharedPtr
std::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:88
InterpCoeff.h
Nektar::SpatialDomains::GeomType
GeomType
Indicates the type of element geometry.
Definition: SpatialDomains.hpp:51
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:376
Nektar::StdRegions::StdExpansion3D::v_HelmholtzMatrixOp_MatFree
virtual void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: StdExpansion3D.cpp:249
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1429
Blas::Daxpy
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
Definition: Blas.hpp:110
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Nektar::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:993
Nektar::SpatialDomains::eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:56
Nektar::StdRegions::StdMatrixKey::GetShapeType
LibUtilities::ShapeType GetShapeType() const
Definition: StdMatrixKey.h:86
Nektar::LocalRegions::PrismExp::CreateStaticCondMatrix
DNekScalBlkMatSharedPtr CreateStaticCondMatrix(const MatrixKey &mkey)
Definition: PrismExp.cpp:1415
Nektar::StdRegions::StdExpansion::GetBoundaryMap
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:812
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:49
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186
Generated on Wed Jan 20 2021 03:33:40 for Nektar++ by   doxygen 1.8.13