Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: PrismExp routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
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  boost::bind(&PrismExp::CreateMatrix, this, _1),
64  std::string("PrismExpMatrix")),
65  m_staticCondMatrixManager(
66  boost::bind(&PrismExp::CreateStaticCondMatrix, this, _1),
67  std::string("PrismExpStaticCondMatrix"))
68  {
69  }
70 
71  PrismExp::PrismExp(const PrismExp &T):
72  StdExpansion(T),
73  StdExpansion3D(T),
74  StdRegions::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  int data_order0 = nummodes[mode_offset];
545  int fillorder0 = min(m_base[0]->GetNumModes(),data_order0);
546  int data_order1 = nummodes[mode_offset+1];
547  int order1 = m_base[1]->GetNumModes();
548  int fillorder1 = min(order1,data_order1);
549  int data_order2 = nummodes[mode_offset+2];
550  int order2 = m_base[2]->GetNumModes();
551  int fillorder2 = min(order2,data_order2);
552 
553  switch(m_base[0]->GetBasisType())
554  {
555  case LibUtilities::eModified_A:
556  {
557  int i,j;
558  int cnt = 0;
559  int cnt1 = 0;
560 
561  ASSERTL1(m_base[1]->GetBasisType() ==
562  LibUtilities::eModified_A,
563  "Extraction routine not set up for this basis");
564  ASSERTL1(m_base[2]->GetBasisType() ==
565  LibUtilities::eModified_B,
566  "Extraction routine not set up for this basis");
567 
568  Vmath::Zero(m_ncoeffs,coeffs,1);
569  for(j = 0; j < fillorder0; ++j)
570  {
571  for(i = 0; i < fillorder1; ++i)
572  {
573  Vmath::Vcopy(fillorder2-j, &data[cnt], 1,
574  &coeffs[cnt1], 1);
575  cnt += data_order2-j;
576  cnt1 += order2-j;
577  }
578 
579  // count out data for j iteration
580  for(i = fillorder1; i < data_order1; ++i)
581  {
582  cnt += data_order2-j;
583  }
584 
585  for(i = fillorder1; i < order1; ++i)
586  {
587  cnt1 += order2-j;
588  }
589  }
590  }
591  break;
592  default:
593  ASSERTL0(false, "basis is either not set up or not "
594  "hierarchicial");
595  }
596  }
597 
598  void PrismExp::v_GetFacePhysMap(const int face,
599  Array<OneD, int> &outarray)
600  {
601  int nquad0 = m_base[0]->GetNumPoints();
602  int nquad1 = m_base[1]->GetNumPoints();
603  int nquad2 = m_base[2]->GetNumPoints();
604  int nq0 = 0;
605  int nq1 = 0;
606 
607  switch(face)
608  {
609  case 0:
610  nq0 = nquad0;
611  nq1 = nquad1;
612  if(outarray.num_elements()!=nq0*nq1)
613  {
614  outarray = Array<OneD, int>(nq0*nq1);
615  }
616 
617  //Directions A and B positive
618  for(int i = 0; i < nquad0*nquad1; ++i)
619  {
620  outarray[i] = i;
621  }
622  break;
623  case 1:
624 
625  nq0 = nquad0;
626  nq1 = nquad2;
627  if(outarray.num_elements()!=nq0*nq1)
628  {
629  outarray = Array<OneD, int>(nq0*nq1);
630  }
631 
632  //Direction A and B positive
633  for (int k=0; k<nquad2; k++)
634  {
635  for(int i = 0; i < nquad0; ++i)
636  {
637  outarray[k*nquad0+i] = (nquad0*nquad1*k)+i;
638  }
639  }
640 
641  break;
642  case 2:
643 
644  nq0 = nquad1;
645  nq1 = nquad2;
646  if(outarray.num_elements()!=nq0*nq1)
647  {
648  outarray = Array<OneD, int>(nq0*nq1);
649  }
650 
651  //Directions A and B positive
652  for(int j = 0; j < nquad1*nquad2; ++j)
653  {
654  outarray[j] = nquad0-1 + j*nquad0;
655 
656  }
657  break;
658  case 3:
659  nq0 = nquad0;
660  nq1 = nquad2;
661  if(outarray.num_elements()!=nq0*nq1)
662  {
663  outarray = Array<OneD, int>(nq0*nq1);
664  }
665 
666  //Direction A and B positive
667  for (int k=0; k<nquad2; k++)
668  {
669  for(int i = 0; i < nquad0; ++i)
670  {
671  outarray[k*nquad0+i] = nquad0*(nquad1-1) + (nquad0*nquad1*k)+i;
672  }
673  }
674  break;
675  case 4:
676 
677  nq0 = nquad1;
678  nq1 = nquad2;
679  if(outarray.num_elements()!=nq0*nq1)
680  {
681  outarray = Array<OneD, int>(nq0*nq1);
682  }
683 
684  //Directions A and B positive
685  for(int j = 0; j < nquad1*nquad2; ++j)
686  {
687  outarray[j] = j*nquad0;
688 
689  }
690  break;
691  default:
692  ASSERTL0(false,"face value (> 4) is out of range");
693  break;
694  }
695 
696  }
697 
698  /** \brief Get the normals along specficied face
699  * Get the face normals interplated to a points0 x points 0
700  * type distribution
701  **/
702  void PrismExp::v_ComputeFaceNormal(const int face)
703  {
704  const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
705  GetGeom()->GetMetricInfo();
706  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
707  SpatialDomains::GeomType type = geomFactors->GetGtype();
708  const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys);
709  const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
710 
711  int nq0 = ptsKeys[0].GetNumPoints();
712  int nq1 = ptsKeys[1].GetNumPoints();
713  int nq2 = ptsKeys[2].GetNumPoints();
714  int nq01 = nq0*nq1;
715  int nqtot;
716 
717 
718  LibUtilities::BasisKey tobasis0 = DetFaceBasisKey(face,0);
719  LibUtilities::BasisKey tobasis1 = DetFaceBasisKey(face,1);
720 
721  // Number of quadrature points in face expansion.
722  int nq_face = tobasis0.GetNumPoints()*tobasis1.GetNumPoints();
723 
724  int vCoordDim = GetCoordim();
725  int i;
726 
727  m_faceNormals[face] = Array<OneD, Array<OneD, NekDouble> >(vCoordDim);
728  Array<OneD, Array<OneD, NekDouble> > &normal = m_faceNormals[face];
729  for (i = 0; i < vCoordDim; ++i)
730  {
731  normal[i] = Array<OneD, NekDouble>(nq_face);
732  }
733 
734  // Regular geometry case
735  if (type == SpatialDomains::eRegular ||
736  type == SpatialDomains::eMovingRegular)
737  {
738  NekDouble fac;
739  // Set up normals
740  switch(face)
741  {
742  case 0:
743  {
744  for(i = 0; i < vCoordDim; ++i)
745  {
746  normal[i][0] = -df[3*i+2][0];;
747  }
748  break;
749  }
750  case 1:
751  {
752  for(i = 0; i < vCoordDim; ++i)
753  {
754  normal[i][0] = -df[3*i+1][0];
755  }
756  break;
757  }
758  case 2:
759  {
760  for(i = 0; i < vCoordDim; ++i)
761  {
762  normal[i][0] = df[3*i][0]+df[3*i+2][0];
763  }
764  break;
765  }
766  case 3:
767  {
768  for(i = 0; i < vCoordDim; ++i)
769  {
770  normal[i][0] = df[3*i+1][0];
771  }
772  break;
773  }
774  case 4:
775  {
776  for(i = 0; i < vCoordDim; ++i)
777  {
778  normal[i][0] = -df[3*i][0];
779  }
780  break;
781  }
782  default:
783  ASSERTL0(false,"face is out of range (face < 4)");
784  }
785 
786  // Normalise resulting vector.
787  fac = 0.0;
788  for(i = 0; i < vCoordDim; ++i)
789  {
790  fac += normal[i][0]*normal[i][0];
791  }
792  fac = 1.0/sqrt(fac);
793  for (i = 0; i < vCoordDim; ++i)
794  {
795  Vmath::Fill(nq_face,fac*normal[i][0],normal[i],1);
796  }
797  }
798  else
799  {
800  // Set up deformed normals.
801  int j, k;
802 
803  // Determine number of quadrature points on the face of 3D elmt
804  if (face == 0)
805  {
806  nqtot = nq0*nq1;
807  }
808  else if (face == 1 || face == 3)
809  {
810  nqtot = nq0*nq2;
811  }
812  else
813  {
814  nqtot = nq1*nq2;
815  }
816 
817  LibUtilities::PointsKey points0;
818  LibUtilities::PointsKey points1;
819 
820  Array<OneD, NekDouble> faceJac(nqtot);
821  Array<OneD, NekDouble> normals(vCoordDim*nqtot,0.0);
822 
823  // Extract Jacobian along face and recover local derivatives
824  // (dx/dr) for polynomial interpolation by multiplying m_gmat by
825  // jacobian
826  switch(face)
827  {
828  case 0:
829  {
830  for(j = 0; j < nq01; ++j)
831  {
832  normals[j] = -df[2][j]*jac[j];
833  normals[nqtot+j] = -df[5][j]*jac[j];
834  normals[2*nqtot+j] = -df[8][j]*jac[j];
835  faceJac[j] = jac[j];
836  }
837 
838  points0 = ptsKeys[0];
839  points1 = ptsKeys[1];
840  break;
841  }
842 
843  case 1:
844  {
845  for (j = 0; j < nq0; ++j)
846  {
847  for(k = 0; k < nq2; ++k)
848  {
849  int tmp = j+nq01*k;
850  normals[j+k*nq0] =
851  -df[1][tmp]*jac[tmp];
852  normals[nqtot+j+k*nq0] =
853  -df[4][tmp]*jac[tmp];
854  normals[2*nqtot+j+k*nq0] =
855  -df[7][tmp]*jac[tmp];
856  faceJac[j+k*nq0] = jac[tmp];
857  }
858  }
859 
860  points0 = ptsKeys[0];
861  points1 = ptsKeys[2];
862  break;
863  }
864 
865  case 2:
866  {
867  for (j = 0; j < nq1; ++j)
868  {
869  for(k = 0; k < nq2; ++k)
870  {
871  int tmp = nq0-1+nq0*j+nq01*k;
872  normals[j+k*nq1] =
873  (df[0][tmp]+df[2][tmp])*jac[tmp];
874  normals[nqtot+j+k*nq1] =
875  (df[3][tmp]+df[5][tmp])*jac[tmp];
876  normals[2*nqtot+j+k*nq1] =
877  (df[6][tmp]+df[8][tmp])*jac[tmp];
878  faceJac[j+k*nq1] = jac[tmp];
879  }
880  }
881 
882  points0 = ptsKeys[1];
883  points1 = ptsKeys[2];
884  break;
885  }
886 
887  case 3:
888  {
889  for (j = 0; j < nq0; ++j)
890  {
891  for(k = 0; k < nq2; ++k)
892  {
893  int tmp = nq0*(nq1-1) + j + nq01*k;
894  normals[j+k*nq0] =
895  df[1][tmp]*jac[tmp];
896  normals[nqtot+j+k*nq0] =
897  df[4][tmp]*jac[tmp];
898  normals[2*nqtot+j+k*nq0] =
899  df[7][tmp]*jac[tmp];
900  faceJac[j+k*nq0] = jac[tmp];
901  }
902  }
903 
904  points0 = ptsKeys[0];
905  points1 = ptsKeys[2];
906  break;
907  }
908 
909  case 4:
910  {
911  for (j = 0; j < nq1; ++j)
912  {
913  for(k = 0; k < nq2; ++k)
914  {
915  int tmp = j*nq0+nq01*k;
916  normals[j+k*nq1] =
917  -df[0][tmp]*jac[tmp];
918  normals[nqtot+j+k*nq1] =
919  -df[3][tmp]*jac[tmp];
920  normals[2*nqtot+j+k*nq1] =
921  -df[6][tmp]*jac[tmp];
922  faceJac[j+k*nq1] = jac[tmp];
923  }
924  }
925 
926  points0 = ptsKeys[1];
927  points1 = ptsKeys[2];
928  break;
929  }
930 
931  default:
932  ASSERTL0(false,"face is out of range (face < 4)");
933  }
934 
935 
936  Array<OneD, NekDouble> work (nq_face, 0.0);
937  // Interpolate Jacobian and invert
938  LibUtilities::Interp2D(points0, points1, faceJac,
939  tobasis0.GetPointsKey(),
940  tobasis1.GetPointsKey(),
941  work);
942  Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
943 
944  // Interpolate normal and multiply by inverse Jacobian.
945  for(i = 0; i < vCoordDim; ++i)
946  {
947  LibUtilities::Interp2D(points0, points1,
948  &normals[i*nqtot],
949  tobasis0.GetPointsKey(),
950  tobasis1.GetPointsKey(),
951  &normal[i][0]);
952  Vmath::Vmul(nq_face,work,1,normal[i],1,normal[i],1);
953  }
954 
955  // Normalise to obtain unit normals.
956  Vmath::Zero(nq_face,work,1);
957  for(i = 0; i < GetCoordim(); ++i)
958  {
959  Vmath::Vvtvp(nq_face,normal[i],1,normal[i],1,work,1,work,1);
960  }
961 
962  Vmath::Vsqrt(nq_face,work,1,work,1);
963  Vmath::Sdiv (nq_face,1.0,work,1,work,1);
964 
965  for(i = 0; i < GetCoordim(); ++i)
966  {
967  Vmath::Vmul(nq_face,normal[i],1,work,1,normal[i],1);
968  }
969  }
970  }
971 
972  void PrismExp::v_MassMatrixOp(
973  const Array<OneD, const NekDouble> &inarray,
974  Array<OneD, NekDouble> &outarray,
975  const StdRegions::StdMatrixKey &mkey)
976  {
977  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
978  }
979 
980  void PrismExp::v_LaplacianMatrixOp(
981  const Array<OneD, const NekDouble> &inarray,
982  Array<OneD, NekDouble> &outarray,
983  const StdRegions::StdMatrixKey &mkey)
984  {
985  PrismExp::LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
986  }
987 
988  void PrismExp::v_LaplacianMatrixOp(
989  const int k1,
990  const int k2,
991  const Array<OneD, const NekDouble> &inarray,
992  Array<OneD, NekDouble> &outarray,
993  const StdRegions::StdMatrixKey &mkey)
994  {
995  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,mkey);
996  }
997 
998  void PrismExp::v_HelmholtzMatrixOp(
999  const Array<OneD, const NekDouble> &inarray,
1000  Array<OneD, NekDouble> &outarray,
1001  const StdRegions::StdMatrixKey &mkey)
1002  {
1003  PrismExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1004  }
1005 
1006  void PrismExp::v_GeneralMatrixOp_MatOp(
1007  const Array<OneD, const NekDouble> &inarray,
1008  Array<OneD, NekDouble> &outarray,
1009  const StdRegions::StdMatrixKey &mkey)
1010  {
1011  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1012 
1013  if(inarray.get() == outarray.get())
1014  {
1015  Array<OneD,NekDouble> tmp(m_ncoeffs);
1016  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1017 
1018  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1019  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1020  }
1021  else
1022  {
1023  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1024  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1025  }
1026  }
1027 
1028  void PrismExp::v_SVVLaplacianFilter(
1029  Array<OneD, NekDouble> &array,
1030  const StdRegions::StdMatrixKey &mkey)
1031  {
1032  int nq = GetTotPoints();
1033 
1034  // Calculate sqrt of the Jacobian
1035  Array<OneD, const NekDouble> jac =
1036  m_metricinfo->GetJac(GetPointsKeys());
1037  Array<OneD, NekDouble> sqrt_jac(nq);
1038  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1039  {
1040  Vmath::Vsqrt(nq,jac,1,sqrt_jac,1);
1041  }
1042  else
1043  {
1044  Vmath::Fill(nq,sqrt(jac[0]),sqrt_jac,1);
1045  }
1046 
1047  // Multiply array by sqrt(Jac)
1048  Vmath::Vmul(nq,sqrt_jac,1,array,1,array,1);
1049 
1050  // Apply std region filter
1051  StdPrismExp::v_SVVLaplacianFilter( array, mkey);
1052 
1053  // Divide by sqrt(Jac)
1054  Vmath::Vdiv(nq,array,1,sqrt_jac,1,array,1);
1055  }
1056 
1057 
1058  //---------------------------------------
1059  // Matrix creation functions
1060  //---------------------------------------
1061 
1062  DNekMatSharedPtr PrismExp::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
1063  {
1064  DNekMatSharedPtr returnval;
1065 
1066  switch(mkey.GetMatrixType())
1067  {
1068  case StdRegions::eHybridDGHelmholtz:
1069  case StdRegions::eHybridDGLamToU:
1070  case StdRegions::eHybridDGLamToQ0:
1071  case StdRegions::eHybridDGLamToQ1:
1072  case StdRegions::eHybridDGLamToQ2:
1073  case StdRegions::eHybridDGHelmBndLam:
1074  case StdRegions::eInvLaplacianWithUnityMean:
1075  returnval = Expansion3D::v_GenMatrix(mkey);
1076  break;
1077  default:
1078  returnval = StdPrismExp::v_GenMatrix(mkey);
1079  break;
1080  }
1081 
1082  return returnval;
1083  }
1084 
1085  DNekMatSharedPtr PrismExp::v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
1086  {
1087  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1088  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1089  LibUtilities::BasisKey bkey2 = m_base[2]->GetBasisKey();
1090  StdRegions::StdPrismExpSharedPtr tmp =
1091  MemoryManager<StdPrismExp>::AllocateSharedPtr(bkey0, bkey1, bkey2);
1092 
1093  return tmp->GetStdMatrix(mkey);
1094  }
1095 
1096  DNekScalMatSharedPtr PrismExp::v_GetLocMatrix(const MatrixKey &mkey)
1097  {
1098  return m_matrixManager[mkey];
1099  }
1100 
1101  DNekScalBlkMatSharedPtr PrismExp::v_GetLocStaticCondMatrix(const MatrixKey &mkey)
1102  {
1103  return m_staticCondMatrixManager[mkey];
1104  }
1105 
1106  void PrismExp::v_DropLocStaticCondMatrix(const MatrixKey &mkey)
1107  {
1108  m_staticCondMatrixManager.DeleteObject(mkey);
1109  }
1110 
1111  DNekScalMatSharedPtr PrismExp::CreateMatrix(const MatrixKey &mkey)
1112  {
1113  DNekScalMatSharedPtr returnval;
1114  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1115 
1116  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,
1117  "Geometric information is not set up");
1118 
1119  switch(mkey.GetMatrixType())
1120  {
1121  case StdRegions::eMass:
1122  {
1123  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1124  {
1125  NekDouble one = 1.0;
1126  DNekMatSharedPtr mat = GenMatrix(mkey);
1127  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1128  }
1129  else
1130  {
1131  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1132  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1133  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1134  }
1135  break;
1136  }
1137 
1138  case StdRegions::eInvMass:
1139  {
1140  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1141  {
1142  NekDouble one = 1.0;
1143  StdRegions::StdMatrixKey masskey(StdRegions::eMass,DetShapeType(),*this);
1144  DNekMatSharedPtr mat = GenMatrix(masskey);
1145  mat->Invert();
1146 
1147  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1148  }
1149  else
1150  {
1151  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1152  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1153  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(fac,mat);
1154  }
1155  break;
1156  }
1157 
1158  case StdRegions::eWeakDeriv0:
1159  case StdRegions::eWeakDeriv1:
1160  case StdRegions::eWeakDeriv2:
1161  {
1162  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1163  {
1164  NekDouble one = 1.0;
1165  DNekMatSharedPtr mat = GenMatrix(mkey);
1166 
1167  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1168  }
1169  else
1170  {
1171  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1172  Array<TwoD, const NekDouble> df =
1173  m_metricinfo->GetDerivFactors(ptsKeys);
1174  int dir = 0;
1175 
1176  switch(mkey.GetMatrixType())
1177  {
1178  case StdRegions::eWeakDeriv0:
1179  dir = 0;
1180  break;
1181  case StdRegions::eWeakDeriv1:
1182  dir = 1;
1183  break;
1184  case StdRegions::eWeakDeriv2:
1185  dir = 2;
1186  break;
1187  default:
1188  break;
1189  }
1190 
1191  MatrixKey deriv0key(StdRegions::eWeakDeriv0,
1192  mkey.GetShapeType(), *this);
1193  MatrixKey deriv1key(StdRegions::eWeakDeriv1,
1194  mkey.GetShapeType(), *this);
1195  MatrixKey deriv2key(StdRegions::eWeakDeriv2,
1196  mkey.GetShapeType(), *this);
1197 
1198  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
1199  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
1200  DNekMat &deriv2 = *GetStdMatrix(deriv2key);
1201 
1202  int rows = deriv0.GetRows();
1203  int cols = deriv1.GetColumns();
1204 
1205  DNekMatSharedPtr WeakDeriv = MemoryManager<DNekMat>
1206  ::AllocateSharedPtr(rows,cols);
1207 
1208  (*WeakDeriv) = df[3*dir ][0]*deriv0
1209  + df[3*dir+1][0]*deriv1
1210  + df[3*dir+2][0]*deriv2;
1211 
1212  returnval = MemoryManager<DNekScalMat>
1213  ::AllocateSharedPtr(jac,WeakDeriv);
1214  }
1215  break;
1216  }
1217 
1218  case StdRegions::eLaplacian:
1219  {
1220  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1221  mkey.GetNVarCoeff() > 0 ||
1222  mkey.ConstFactorExists(
1223  StdRegions::eFactorSVVCutoffRatio))
1224  {
1225  NekDouble one = 1.0;
1226  DNekMatSharedPtr mat = GenMatrix(mkey);
1227  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1228  }
1229  else
1230  {
1231  MatrixKey lap00key(StdRegions::eLaplacian00,
1232  mkey.GetShapeType(), *this);
1233  MatrixKey lap01key(StdRegions::eLaplacian01,
1234  mkey.GetShapeType(), *this);
1235  MatrixKey lap02key(StdRegions::eLaplacian02,
1236  mkey.GetShapeType(), *this);
1237  MatrixKey lap11key(StdRegions::eLaplacian11,
1238  mkey.GetShapeType(), *this);
1239  MatrixKey lap12key(StdRegions::eLaplacian12,
1240  mkey.GetShapeType(), *this);
1241  MatrixKey lap22key(StdRegions::eLaplacian22,
1242  mkey.GetShapeType(), *this);
1243 
1244  DNekMat &lap00 = *GetStdMatrix(lap00key);
1245  DNekMat &lap01 = *GetStdMatrix(lap01key);
1246  DNekMat &lap02 = *GetStdMatrix(lap02key);
1247  DNekMat &lap11 = *GetStdMatrix(lap11key);
1248  DNekMat &lap12 = *GetStdMatrix(lap12key);
1249  DNekMat &lap22 = *GetStdMatrix(lap22key);
1250 
1251  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1252  Array<TwoD, const NekDouble> gmat
1253  = m_metricinfo->GetGmat(ptsKeys);
1254 
1255  int rows = lap00.GetRows();
1256  int cols = lap00.GetColumns();
1257 
1258  DNekMatSharedPtr lap = MemoryManager<DNekMat>
1259  ::AllocateSharedPtr(rows,cols);
1260 
1261  (*lap) = gmat[0][0]*lap00
1262  + gmat[4][0]*lap11
1263  + gmat[8][0]*lap22
1264  + gmat[3][0]*(lap01 + Transpose(lap01))
1265  + gmat[6][0]*(lap02 + Transpose(lap02))
1266  + gmat[7][0]*(lap12 + Transpose(lap12));
1267 
1268  returnval = MemoryManager<DNekScalMat>
1269  ::AllocateSharedPtr(jac,lap);
1270  }
1271  break;
1272  }
1273 
1274  case StdRegions::eHelmholtz:
1275  {
1276  NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorLambda);
1277  MatrixKey masskey(StdRegions::eMass,
1278  mkey.GetShapeType(), *this);
1279  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1280  MatrixKey lapkey(StdRegions::eLaplacian,
1281  mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1282  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1283 
1284  int rows = LapMat.GetRows();
1285  int cols = LapMat.GetColumns();
1286 
1287  DNekMatSharedPtr helm = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1288 
1289  NekDouble one = 1.0;
1290  (*helm) = LapMat + factor*MassMat;
1291 
1292  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
1293  break;
1294  }
1295  case StdRegions::eHybridDGHelmholtz:
1296  case StdRegions::eHybridDGLamToU:
1297  case StdRegions::eHybridDGLamToQ0:
1298  case StdRegions::eHybridDGLamToQ1:
1299  case StdRegions::eHybridDGHelmBndLam:
1300  case StdRegions::eInvLaplacianWithUnityMean:
1301  {
1302  NekDouble one = 1.0;
1303 
1304  DNekMatSharedPtr mat = GenMatrix(mkey);
1305  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1306 
1307  break;
1308  }
1309 
1310  case StdRegions::eInvHybridDGHelmholtz:
1311  {
1312  NekDouble one = 1.0;
1313 
1314  MatrixKey hkey(StdRegions::eHybridDGHelmholtz, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1315 // StdRegions::StdMatrixKey hkey(StdRegions::eHybridDGHelmholtz,
1316 // DetExpansionType(),*this,
1317 // mkey.GetConstant(0),
1318 // mkey.GetConstant(1));
1319  DNekMatSharedPtr mat = GenMatrix(hkey);
1320 
1321  mat->Invert();
1322  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1323  break;
1324  }
1325 
1326  case StdRegions::eIProductWRTBase:
1327  {
1328  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1329  {
1330  NekDouble one = 1.0;
1331  DNekMatSharedPtr mat = GenMatrix(mkey);
1332  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1333  }
1334  else
1335  {
1336  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1337  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1338  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1339  }
1340  break;
1341  }
1342  case StdRegions::ePreconLinearSpace:
1343  {
1344  NekDouble one = 1.0;
1345  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1346  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1347  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1348  DNekMatSharedPtr R=BuildVertexMatrix(A);
1349 
1350  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1351  }
1352  break;
1353  case StdRegions::ePreconLinearSpaceMass:
1354  {
1355  NekDouble one = 1.0;
1356  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1357  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1358  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1359  DNekMatSharedPtr R=BuildVertexMatrix(A);
1360 
1361  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1362  }
1363  break;
1364  case StdRegions::ePreconR:
1365  {
1366  NekDouble one = 1.0;
1367  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this,mkey.GetConstFactors(), mkey.GetVarCoeffs());
1368  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1369  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1370 
1371  DNekScalMatSharedPtr Atmp;
1372  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1373 
1374  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1375  }
1376  break;
1377  case StdRegions::ePreconRT:
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  case StdRegions::ePreconRTMass:
1404  {
1405  NekDouble one = 1.0;
1406  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1407  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1408  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1409 
1410  DNekScalMatSharedPtr Atmp;
1411  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1412 
1413  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1414  }
1415  break;
1416  default:
1417  {
1418  NekDouble one = 1.0;
1419  DNekMatSharedPtr mat = GenMatrix(mkey);
1420 
1421  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1422  }
1423  }
1424 
1425  return returnval;
1426  }
1427 
1428  DNekScalBlkMatSharedPtr PrismExp::CreateStaticCondMatrix(const MatrixKey &mkey)
1429  {
1430  DNekScalBlkMatSharedPtr returnval;
1431 
1432  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
1433 
1434  // set up block matrix system
1435  unsigned int nbdry = NumBndryCoeffs();
1436  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
1437  unsigned int exp_size[] = {nbdry, nint};
1438  unsigned int nblks=2;
1439  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks, nblks, exp_size, exp_size);
1440  NekDouble factor = 1.0;
1441 
1442  switch(mkey.GetMatrixType())
1443  {
1444  case StdRegions::eLaplacian:
1445  case StdRegions::eHelmholtz: // special case since Helmholtz not defined in StdRegions
1446  // use Deformed case for both regular and deformed geometries
1447  factor = 1.0;
1448  goto UseLocRegionsMatrix;
1449  break;
1450  default:
1451  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1452  {
1453  factor = 1.0;
1454  goto UseLocRegionsMatrix;
1455  }
1456  else
1457  {
1458  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1459  factor = mat->Scale();
1460  goto UseStdRegionsMatrix;
1461  }
1462  break;
1463  UseStdRegionsMatrix:
1464  {
1465  NekDouble invfactor = 1.0/factor;
1466  NekDouble one = 1.0;
1467  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
1468  DNekScalMatSharedPtr Atmp;
1469  DNekMatSharedPtr Asubmat;
1470 
1471  //TODO: check below
1472  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
1473  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
1474  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
1475  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
1476  }
1477  break;
1478  UseLocRegionsMatrix:
1479  {
1480  int i,j;
1481  NekDouble invfactor = 1.0/factor;
1482  NekDouble one = 1.0;
1483  DNekScalMat &mat = *GetLocMatrix(mkey);
1484  DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
1485  DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
1486  DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
1487  DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
1488 
1489  Array<OneD,unsigned int> bmap(nbdry);
1490  Array<OneD,unsigned int> imap(nint);
1491  GetBoundaryMap(bmap);
1492  GetInteriorMap(imap);
1493 
1494  for(i = 0; i < nbdry; ++i)
1495  {
1496  for(j = 0; j < nbdry; ++j)
1497  {
1498  (*A)(i,j) = mat(bmap[i],bmap[j]);
1499  }
1500 
1501  for(j = 0; j < nint; ++j)
1502  {
1503  (*B)(i,j) = mat(bmap[i],imap[j]);
1504  }
1505  }
1506 
1507  for(i = 0; i < nint; ++i)
1508  {
1509  for(j = 0; j < nbdry; ++j)
1510  {
1511  (*C)(i,j) = mat(imap[i],bmap[j]);
1512  }
1513 
1514  for(j = 0; j < nint; ++j)
1515  {
1516  (*D)(i,j) = mat(imap[i],imap[j]);
1517  }
1518  }
1519 
1520  // Calculate static condensed system
1521  if(nint)
1522  {
1523  D->Invert();
1524  (*B) = (*B)*(*D);
1525  (*A) = (*A) - (*B)*(*C);
1526  }
1527 
1528  DNekScalMatSharedPtr Atmp;
1529 
1530  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
1531  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
1532  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
1533  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
1534 
1535  }
1536  break;
1537  }
1538  return returnval;
1539  }
1540 
1541 
1542  /**
1543  * @brief Calculate the Laplacian multiplication in a matrix-free
1544  * manner.
1545  *
1546  * This function is the kernel of the Laplacian matrix-free operator,
1547  * and is used in #v_HelmholtzMatrixOp_MatFree to determine the effect
1548  * of the Helmholtz operator in a similar fashion.
1549  *
1550  * The majority of the calculation is precisely the same as in the
1551  * hexahedral expansion; however the collapsed co-ordinate system must
1552  * be taken into account when constructing the geometric factors. How
1553  * this is done is detailed more exactly in the tetrahedral expansion.
1554  * On entry to this function, the input #inarray must be in its
1555  * backwards-transformed state (i.e. \f$\mathbf{u} =
1556  * \mathbf{B}\hat{\mathbf{u}}\f$). The output is in coefficient space.
1557  *
1558  * @see %TetExp::v_HelmholtzMatrixOp_MatFree
1559  */
1560  void PrismExp::v_LaplacianMatrixOp_MatFree_Kernel(
1561  const Array<OneD, const NekDouble> &inarray,
1562  Array<OneD, NekDouble> &outarray,
1563  Array<OneD, NekDouble> &wsp)
1564  {
1565  int nquad0 = m_base[0]->GetNumPoints();
1566  int nquad1 = m_base[1]->GetNumPoints();
1567  int nquad2 = m_base[2]->GetNumPoints();
1568  int nqtot = nquad0*nquad1*nquad2;
1569  int i;
1570 
1571  // Set up temporary storage.
1572  Array<OneD,NekDouble> alloc(11*nqtot,0.0);
1573  Array<OneD,NekDouble> wsp1 (alloc ); // TensorDeriv 1
1574  Array<OneD,NekDouble> wsp2 (alloc+ 1*nqtot); // TensorDeriv 2
1575  Array<OneD,NekDouble> wsp3 (alloc+ 2*nqtot); // TensorDeriv 3
1576  Array<OneD,NekDouble> g0 (alloc+ 3*nqtot); // g0
1577  Array<OneD,NekDouble> g1 (alloc+ 4*nqtot); // g1
1578  Array<OneD,NekDouble> g2 (alloc+ 5*nqtot); // g2
1579  Array<OneD,NekDouble> g3 (alloc+ 6*nqtot); // g3
1580  Array<OneD,NekDouble> g4 (alloc+ 7*nqtot); // g4
1581  Array<OneD,NekDouble> g5 (alloc+ 8*nqtot); // g5
1582  Array<OneD,NekDouble> h0 (alloc+ 3*nqtot); // h0 == g0
1583  Array<OneD,NekDouble> h1 (alloc+ 6*nqtot); // h1 == g3
1584  Array<OneD,NekDouble> wsp4 (alloc+ 4*nqtot); // wsp4 == g1
1585  Array<OneD,NekDouble> wsp5 (alloc+ 5*nqtot); // wsp5 == g2
1586  Array<OneD,NekDouble> wsp6 (alloc+ 8*nqtot); // wsp6 == g5
1587  Array<OneD,NekDouble> wsp7 (alloc+ 3*nqtot); // wsp7 == g0
1588  Array<OneD,NekDouble> wsp8 (alloc+ 9*nqtot); // wsp8
1589  Array<OneD,NekDouble> wsp9 (alloc+10*nqtot); // wsp9
1590 
1591  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
1592  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
1593  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata();
1594  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
1595  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
1596  const Array<OneD, const NekDouble>& dbase2 = m_base[2]->GetDbdata();
1597 
1598  // Step 1. LAPLACIAN MATRIX OPERATION
1599  // wsp1 = du_dxi1 = D_xi1 * wsp0 = D_xi1 * u
1600  // wsp2 = du_dxi2 = D_xi2 * wsp0 = D_xi2 * u
1601  // wsp3 = du_dxi3 = D_xi3 * wsp0 = D_xi3 * u
1602  StdExpansion3D::PhysTensorDeriv(inarray,wsp1,wsp2,wsp3);
1603 
1604  const Array<TwoD, const NekDouble>& df =
1605  m_metricinfo->GetDerivFactors(GetPointsKeys());
1606  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
1607  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1608 
1609  // Step 2. Calculate the metric terms of the collapsed
1610  // coordinate transformation (Spencer's book P152)
1611  for (i = 0; i < nquad2; ++i)
1612  {
1613  Vmath::Fill(nquad0*nquad1, 2.0/(1.0-z2[i]), &h0[0]+i*nquad0*nquad1,1);
1614  Vmath::Fill(nquad0*nquad1, 2.0/(1.0-z2[i]), &h1[0]+i*nquad0*nquad1,1);
1615  }
1616  for (i = 0; i < nquad0; i++)
1617  {
1618  Blas::Dscal(nquad1*nquad2, 0.5*(1+z0[i]), &h1[0]+i, nquad0);
1619  }
1620 
1621  // Step 3. Construct combined metric terms for physical space to
1622  // collapsed coordinate system. Order of construction optimised
1623  // to minimise temporary storage
1624  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1625  {
1626  // wsp4 = d eta_1/d x_1
1627  Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &df[2][0], 1, &h1[0], 1, &wsp4[0], 1);
1628  // wsp5 = d eta_2/d x_1
1629  Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &df[5][0], 1, &h1[0], 1, &wsp5[0], 1);
1630  // wsp6 = d eta_3/d x_1d
1631  Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &df[8][0], 1, &h1[0], 1, &wsp6[0], 1);
1632 
1633  // g0 (overwrites h0)
1634  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0], 1, &g0[0], 1);
1635  Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1636 
1637  // g3 (overwrites h1)
1638  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &wsp4[0], 1, &df[4][0], 1, &wsp5[0], 1, &g3[0], 1);
1639  Vmath::Vvtvp (nqtot, &df[7][0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1640 
1641  // g4
1642  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0], 1, &g4[0], 1);
1643  Vmath::Vvtvp (nqtot, &df[8][0], 1, &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1644 
1645  // Overwrite wsp4/5/6 with g1/2/5
1646  // g1
1647  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &df[1][0], 1, &df[4][0], 1, &df[4][0], 1, &g1[0], 1);
1648  Vmath::Vvtvp (nqtot, &df[7][0], 1, &df[7][0], 1, &g1[0], 1, &g1[0], 1);
1649 
1650  // g2
1651  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &df[2][0], 1, &df[5][0], 1, &df[5][0], 1, &g2[0], 1);
1652  Vmath::Vvtvp (nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1653 
1654  // g5
1655  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &df[2][0], 1, &df[4][0], 1, &df[5][0], 1, &g5[0], 1);
1656  Vmath::Vvtvp (nqtot, &df[7][0], 1, &df[8][0], 1, &g5[0], 1, &g5[0], 1);
1657  }
1658  else
1659  {
1660  // wsp4 = d eta_1/d x_1
1661  Vmath::Svtsvtp(nqtot, df[0][0], &h0[0], 1, df[2][0], &h1[0], 1, &wsp4[0], 1);
1662  // wsp5 = d eta_2/d x_1
1663  Vmath::Svtsvtp(nqtot, df[3][0], &h0[0], 1, df[5][0], &h1[0], 1, &wsp5[0], 1);
1664  // wsp6 = d eta_3/d x_1
1665  Vmath::Svtsvtp(nqtot, df[6][0], &h0[0], 1, df[8][0], &h1[0], 1, &wsp6[0], 1);
1666 
1667  // g0 (overwrites h0)
1668  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0], 1, &g0[0], 1);
1669  Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1670 
1671  // g3 (overwrites h1)
1672  Vmath::Svtsvtp(nqtot, df[1][0], &wsp4[0], 1, df[4][0], &wsp5[0], 1, &g3[0], 1);
1673  Vmath::Svtvp (nqtot, df[7][0], &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1674 
1675  // g4
1676  Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1, &g4[0], 1);
1677  Vmath::Svtvp (nqtot, df[8][0], &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1678 
1679  // Overwrite wsp4/5/6 with g1/2/5
1680  // g1
1681  Vmath::Fill(nqtot, df[1][0]*df[1][0] + df[4][0]*df[4][0] + df[7][0]*df[7][0], &g1[0], 1);
1682 
1683  // g2
1684  Vmath::Fill(nqtot, df[2][0]*df[2][0] + df[5][0]*df[5][0] + df[8][0]*df[8][0], &g2[0], 1);
1685 
1686  // g5
1687  Vmath::Fill(nqtot, df[1][0]*df[2][0] + df[4][0]*df[5][0] + df[7][0]*df[8][0], &g5[0], 1);
1688  }
1689  // Compute component derivatives into wsp7, 8, 9 (wsp7 overwrites
1690  // g0).
1691  Vmath::Vvtvvtp(nqtot,&g0[0],1,&wsp1[0],1,&g3[0],1,&wsp2[0],1,&wsp7[0],1);
1692  Vmath::Vvtvp (nqtot,&g4[0],1,&wsp3[0],1,&wsp7[0],1,&wsp7[0],1);
1693  Vmath::Vvtvvtp(nqtot,&g1[0],1,&wsp2[0],1,&g3[0],1,&wsp1[0],1,&wsp8[0],1);
1694  Vmath::Vvtvp (nqtot,&g5[0],1,&wsp3[0],1,&wsp8[0],1,&wsp8[0],1);
1695  Vmath::Vvtvvtp(nqtot,&g2[0],1,&wsp3[0],1,&g4[0],1,&wsp1[0],1,&wsp9[0],1);
1696  Vmath::Vvtvp (nqtot,&g5[0],1,&wsp2[0],1,&wsp9[0],1,&wsp9[0],1);
1697 
1698  // Step 4.
1699  // Multiply by quadrature metric
1700  MultiplyByQuadratureMetric(wsp7,wsp7);
1701  MultiplyByQuadratureMetric(wsp8,wsp8);
1702  MultiplyByQuadratureMetric(wsp9,wsp9);
1703 
1704  // Perform inner product w.r.t derivative bases.
1705  IProductWRTBase_SumFacKernel(dbase0,base1,base2,wsp7,wsp1, wsp,false,true,true);
1706  IProductWRTBase_SumFacKernel(base0,dbase1,base2,wsp8,wsp2, wsp,true,false,true);
1707  IProductWRTBase_SumFacKernel(base0,base1,dbase2,wsp9,outarray,wsp,true,true,false);
1708 
1709  // Step 5.
1710  // Sum contributions from wsp1, wsp2 and outarray.
1711  Vmath::Vadd(m_ncoeffs,wsp1.get(),1,outarray.get(),1,outarray.get(),1);
1712  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
1713  }
1714 
1715 
1716  void PrismExp::v_GetSimplexEquiSpacedConnectivity(
1717  Array<OneD, int> &conn,
1718  bool oldstandard)
1719  {
1720  int np0 = m_base[0]->GetNumPoints();
1721  int np1 = m_base[1]->GetNumPoints();
1722  int np2 = m_base[2]->GetNumPoints();
1723  int np = max(np0,max(np1,np2));
1724  Array<OneD, int> prismpt(6);
1725  bool standard = true;
1726 
1727  int vid0 = m_geom->GetVid(0);
1728  int vid1 = m_geom->GetVid(1);
1729  int vid2 = m_geom->GetVid(4);
1730  int rotate = 0;
1731 
1732  // sort out prism rotation according to
1733  if((vid2 < vid1)&&(vid2 < vid0)) // top triangle vertex is lowest id
1734  {
1735  rotate = 0;
1736  if(vid0 > vid1)
1737  {
1738  standard = false;// reverse base direction
1739  }
1740  }
1741  else if((vid1 < vid2)&&(vid1 < vid0))
1742  {
1743  rotate = 1;
1744  if(vid2 > vid0)
1745  {
1746  standard = false;// reverse base direction
1747  }
1748  }
1749  else if ((vid0 < vid2)&&(vid0 < vid1))
1750  {
1751  rotate = 2;
1752  if(vid1 > vid2)
1753  {
1754  standard = false; // reverse base direction
1755  }
1756  }
1757 
1758  conn = Array<OneD, int>(12*(np-1)*(np-1)*(np-1));
1759 
1760  int row = 0;
1761  int rowp1 = 0;
1762  int plane = 0;
1763  int row1 = 0;
1764  int row1p1 = 0;
1765  int planep1 = 0;
1766  int cnt = 0;
1767 
1768 
1769  Array<OneD, int> rot(3);
1770 
1771  rot[0] = (0+rotate)%3;
1772  rot[1] = (1+rotate)%3;
1773  rot[2] = (2+rotate)%3;
1774 
1775  // lower diagonal along 1-3 on base
1776  for(int i = 0; i < np-1; ++i)
1777  {
1778  planep1 += (np-i)*np;
1779  row = 0; // current plane row offset
1780  rowp1 = 0; // current plane row plus one offset
1781  row1 = 0; // next plane row offset
1782  row1p1 = 0; // nex plane row plus one offset
1783  if(standard == false)
1784  {
1785  for(int j = 0; j < np-1; ++j)
1786  {
1787  rowp1 += np-i;
1788  row1p1 += np-i-1;
1789  for(int k = 0; k < np-i-2; ++k)
1790  {
1791  // bottom prism block
1792  prismpt[rot[0]] = plane + row + k;
1793  prismpt[rot[1]] = plane + row + k+1;
1794  prismpt[rot[2]] = planep1 + row1 + k;
1795 
1796  prismpt[3+rot[0]] = plane + rowp1 + k;
1797  prismpt[3+rot[1]] = plane + rowp1 + k+1;
1798  prismpt[3+rot[2]] = planep1 + row1p1 + k;
1799 
1800  conn[cnt++] = prismpt[0];
1801  conn[cnt++] = prismpt[1];
1802  conn[cnt++] = prismpt[3];
1803  conn[cnt++] = prismpt[2];
1804 
1805  conn[cnt++] = prismpt[5];
1806  conn[cnt++] = prismpt[2];
1807  conn[cnt++] = prismpt[3];
1808  conn[cnt++] = prismpt[4];
1809 
1810  conn[cnt++] = prismpt[3];
1811  conn[cnt++] = prismpt[1];
1812  conn[cnt++] = prismpt[4];
1813  conn[cnt++] = prismpt[2];
1814 
1815  // upper prism block.
1816  prismpt[rot[0]] = planep1 + row1 + k+1;
1817  prismpt[rot[1]] = planep1 + row1 + k;
1818  prismpt[rot[2]] = plane + row + k+1;
1819 
1820  prismpt[3+rot[0]] = planep1 + row1p1 + k+1;
1821  prismpt[3+rot[1]] = planep1 + row1p1 + k;
1822  prismpt[3+rot[2]] = plane + rowp1 + k+1;
1823 
1824 
1825  conn[cnt++] = prismpt[0];
1826  conn[cnt++] = prismpt[1];
1827  conn[cnt++] = prismpt[2];
1828  conn[cnt++] = prismpt[5];
1829 
1830  conn[cnt++] = prismpt[5];
1831  conn[cnt++] = prismpt[0];
1832  conn[cnt++] = prismpt[4];
1833  conn[cnt++] = prismpt[1];
1834 
1835  conn[cnt++] = prismpt[3];
1836  conn[cnt++] = prismpt[4];
1837  conn[cnt++] = prismpt[0];
1838  conn[cnt++] = prismpt[5];
1839 
1840  }
1841 
1842  // bottom prism block
1843  prismpt[rot[0]] = plane + row + np-i-2;
1844  prismpt[rot[1]] = plane + row + np-i-1;
1845  prismpt[rot[2]] = planep1 + row1 + np-i-2;
1846 
1847  prismpt[3+rot[0]] = plane + rowp1 + np-i-2;
1848  prismpt[3+rot[1]] = plane + rowp1 + np-i-1;
1849  prismpt[3+rot[2]] = planep1 + row1p1 + np-i-2;
1850 
1851  conn[cnt++] = prismpt[0];
1852  conn[cnt++] = prismpt[1];
1853  conn[cnt++] = prismpt[3];
1854  conn[cnt++] = prismpt[2];
1855 
1856  conn[cnt++] = prismpt[5];
1857  conn[cnt++] = prismpt[2];
1858  conn[cnt++] = prismpt[3];
1859  conn[cnt++] = prismpt[4];
1860 
1861  conn[cnt++] = prismpt[3];
1862  conn[cnt++] = prismpt[1];
1863  conn[cnt++] = prismpt[4];
1864  conn[cnt++] = prismpt[2];
1865 
1866  row += np-i;
1867  row1 += np-i-1;
1868  }
1869 
1870  }
1871  else
1872  { // lower diagonal along 0-4 on base
1873  for(int j = 0; j < np-1; ++j)
1874  {
1875  rowp1 += np-i;
1876  row1p1 += np-i-1;
1877  for(int k = 0; k < np-i-2; ++k)
1878  {
1879  // bottom prism block
1880  prismpt[rot[0]] = plane + row + k;
1881  prismpt[rot[1]] = plane + row + k+1;
1882  prismpt[rot[2]] = planep1 + row1 + k;
1883 
1884  prismpt[3+rot[0]] = plane + rowp1 + k;
1885  prismpt[3+rot[1]] = plane + rowp1 + k+1;
1886  prismpt[3+rot[2]] = planep1 + row1p1 + k;
1887 
1888  conn[cnt++] = prismpt[0];
1889  conn[cnt++] = prismpt[1];
1890  conn[cnt++] = prismpt[4];
1891  conn[cnt++] = prismpt[2];
1892 
1893  conn[cnt++] = prismpt[4];
1894  conn[cnt++] = prismpt[3];
1895  conn[cnt++] = prismpt[0];
1896  conn[cnt++] = prismpt[2];
1897 
1898  conn[cnt++] = prismpt[3];
1899  conn[cnt++] = prismpt[4];
1900  conn[cnt++] = prismpt[5];
1901  conn[cnt++] = prismpt[2];
1902 
1903  // upper prism block.
1904  prismpt[rot[0]] = planep1 + row1 + k+1;
1905  prismpt[rot[1]] = planep1 + row1 + k;
1906  prismpt[rot[2]] = plane + row + k+1;
1907 
1908  prismpt[3+rot[0]] = planep1 + row1p1 + k+1;
1909  prismpt[3+rot[1]] = planep1 + row1p1 + k;
1910  prismpt[3+rot[2]] = plane + rowp1 + k+1;
1911 
1912  conn[cnt++] = prismpt[0];
1913  conn[cnt++] = prismpt[2];
1914  conn[cnt++] = prismpt[1];
1915  conn[cnt++] = prismpt[5];
1916 
1917  conn[cnt++] = prismpt[3];
1918  conn[cnt++] = prismpt[5];
1919  conn[cnt++] = prismpt[0];
1920  conn[cnt++] = prismpt[1];
1921 
1922  conn[cnt++] = prismpt[5];
1923  conn[cnt++] = prismpt[3];
1924  conn[cnt++] = prismpt[4];
1925  conn[cnt++] = prismpt[1];
1926  }
1927 
1928  // bottom prism block
1929  prismpt[rot[0]] = plane + row + np-i-2;
1930  prismpt[rot[1]] = plane + row + np-i-1;
1931  prismpt[rot[2]] = planep1 + row1 + np-i-2;
1932 
1933  prismpt[3+rot[0]] = plane + rowp1 + np-i-2;
1934  prismpt[3+rot[1]] = plane + rowp1 + np-i-1;
1935  prismpt[3+rot[2]] = planep1 + row1p1 + np-i-2;
1936 
1937  conn[cnt++] = prismpt[0];
1938  conn[cnt++] = prismpt[1];
1939  conn[cnt++] = prismpt[4];
1940  conn[cnt++] = prismpt[2];
1941 
1942  conn[cnt++] = prismpt[4];
1943  conn[cnt++] = prismpt[3];
1944  conn[cnt++] = prismpt[0];
1945  conn[cnt++] = prismpt[2];
1946 
1947  conn[cnt++] = prismpt[3];
1948  conn[cnt++] = prismpt[4];
1949  conn[cnt++] = prismpt[5];
1950  conn[cnt++] = prismpt[2];
1951 
1952  row += np-i;
1953  row1 += np-i-1;
1954  }
1955 
1956  }
1957  plane += (np-i)*np;
1958  }
1959  }
1960 
1961  }//end of namespace
1962 }//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:998
Nektar::StdRegions::StdExpansion::GetPointsKeys
const LibUtilities::PointsKeyVector GetPointsKeys() const
Definition: StdExpansion.h:1323
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:470
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
Nektar::StdRegions::StdExpansion::GenMatrix
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:1054
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Nektar::StdRegions::StdMatrixKey::GetConstFactors
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:142
Nektar::LocalRegions::PrismExp::m_matrixManager
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: PrismExp.h:211
Nektar::StdRegions::StdMatrixKey::GetVarCoeffs
const VarCoeffMap & GetVarCoeffs() const
Definition: StdMatrixKey.h:168
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:242
Nektar::LocalRegions::Expansion::BuildTransformationMatrix
DNekMatSharedPtr BuildTransformationMatrix(const DNekScalMatSharedPtr &r_bnd, const StdRegions::MatrixType matrixType)
Definition: Expansion.cpp:90
Nektar::StdRegions::StdMatrixKey::GetMatrixType
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:408
Nektar::StdRegions::StdExpansion::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:947
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::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
Nektar::StdRegions::StdPrismExpSharedPtr
boost::shared_ptr< StdPrismExp > StdPrismExpSharedPtr
Definition: StdPrismExp.h:258
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:1006
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
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:1485
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:485
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:442
Nektar::LibUtilities::eModified_A
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::LocalRegions::Expansion::m_metricinfo
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:131
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:271
Nektar::LocalRegions::PrismExp::v_GetSimplexEquiSpacedConnectivity
virtual void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true)
Definition: PrismExp.cpp:1716
Nektar::LocalRegions::Expansion::BuildVertexMatrix
DNekMatSharedPtr BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
Definition: Expansion.cpp:98
Nektar::StdRegions::StdMatrixKey::GetShapeType
LibUtilities::ShapeType GetShapeType() const
Definition: StdMatrixKey.h:87
Nektar::LocalRegions::PrismExp::v_GetLocStaticCondMatrix
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: PrismExp.cpp:1101
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:972
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:241
Nektar::LocalRegions::Expansion::m_geom
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:130
Nektar::StdRegions::StdExpansion::GetLocStaticCondMatrix
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
Definition: StdExpansion.h:753
Nektar::DNekMatSharedPtr
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
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::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:706
Nektar::DNekScalMatSharedPtr
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:73
Nektar::LocalRegions::PrismExp::v_GetFacePhysMap
virtual void v_GetFacePhysMap(const int face, Array< OneD, int > &outarray)
Definition: PrismExp.cpp:598
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:131
Nektar::LocalRegions::PrismExp::CreateMatrix
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition: PrismExp.cpp:1111
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
SegGeom.h
Nektar::LibUtilities::Interp2D
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis...
Definition: Interp.cpp:116
Nektar::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:1085
Nektar::StdRegions::StdExpansion::GetStdStaticCondMatrix
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:711
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:127
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:1560
Nektar::LocalRegions::PrismExp::m_staticCondMatrixManager
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: PrismExp.h:213
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:213
Nektar::LocalRegions::PrismExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: PrismExp.cpp:1062
PrismExp.h
Nektar::LibUtilities::eModified_B
Principle Modified Functions .
Definition: BasisType.h:50
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:224
Nektar::LibUtilities::BasisKey::GetNumPoints
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:128
Nektar::DNekScalBlkMatSharedPtr
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
Nektar::LocalRegions::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:1028
Nektar::LocalRegions::PrismExp::v_DropLocStaticCondMatrix
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: PrismExp.cpp:1106
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:702
Nektar::StdRegions::StdExpansion::GetInteriorMap
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:819
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:219
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Interp.h
Nektar::LocalRegions::PrismExp::v_GetLinStdExp
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
Definition: PrismExp.cpp:459
Nektar::StdRegions::StdExpansion3D::m_faceNormals
std::map< int, NormalVector > m_faceNormals
Definition: StdExpansion3D.h:204
Nektar::StdRegions::StdExpansion::GetBasisType
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:165
Nektar::LocalRegions::PrismExp::v_GetLocMatrix
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: PrismExp.cpp:1096
Nektar::DNekBlkMatSharedPtr
boost::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:72
Nektar::LocalRegions::Expansion::GetLocMatrix
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:85
Nektar::LocalRegions::Expansion3D::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: Expansion3D.cpp:383
Nektar::LibUtilities::BasisKey::GetPointsKey
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:145
Nektar::LocalRegions::Expansion::GetGeom
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:161
Nektar::SpatialDomains::GeomFactorsSharedPtr
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
Vmath::Vvtvvtp
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:537
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:250
Nektar::SpatialDomains::eRegular
Geometry is straight-sided with constant geometric factors.
Definition: SpatialDomains.hpp:55
Nektar::SpatialDomains::PrismGeomSharedPtr
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:109
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::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:131
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:591
Nektar::StdRegions::StdExpansion::DetFaceBasisKey
const LibUtilities::BasisKey DetFaceBasisKey(const int i, const int k) const
Definition: StdExpansion.h:324
InterpCoeff.h
Nektar::SpatialDomains::GeomType
GeomType
Indicates the type of element geometry.
Definition: SpatialDomains.hpp:52
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1878
Nektar::StdRegions::StdExpansion3D::v_HelmholtzMatrixOp_MatFree
virtual void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: StdExpansion3D.cpp:248
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1418
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
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:980
Nektar::SpatialDomains::eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:57
Nektar::LocalRegions::PrismExp::v_GetStdExp
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
Definition: PrismExp.cpp:451
Nektar::LocalRegions::PrismExp::CreateStaticCondMatrix
DNekScalBlkMatSharedPtr CreateStaticCondMatrix(const MatrixKey &mkey)
Definition: PrismExp.cpp:1428
Nektar::StdRegions::StdExpansion::GetBoundaryMap
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:814
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
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:299
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:183
Generated on Sat Nov 25 2017 16:18:15 for Nektar++ by   doxygen 1.8.8