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