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