Nektar++
PrismExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File PrismExp.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // 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  ///Returns the physical values at the quadrature points of a face
583  void PrismExp::v_GetTracePhysVals(
584  const int face,
585  const StdRegions::StdExpansionSharedPtr &FaceExp,
586  const Array<OneD, const NekDouble> &inarray,
587  Array<OneD, NekDouble> &outarray,
588  StdRegions::Orientation orient)
589  {
590  v_GetFacePhysVals(face,FaceExp,inarray,outarray,orient);
591  }
592 
593  void PrismExp::v_GetFacePhysVals(
594  const int face,
595  const StdRegions::StdExpansionSharedPtr &FaceExp,
596  const Array<OneD, const NekDouble> &inarray,
597  Array<OneD, NekDouble> &outarray,
598  StdRegions::Orientation orient)
599  {
600  int nquad0 = m_base[0]->GetNumPoints();
601  int nquad1 = m_base[1]->GetNumPoints();
602  int nquad2 = m_base[2]->GetNumPoints();
603 
604  Array<OneD,NekDouble> o_tmp(GetFaceNumPoints(face));
605 
606  if (orient == StdRegions::eNoOrientation)
607  {
608  orient = GetForient(face);
609  }
610 
611  switch(face)
612  {
613  case 0:
614  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
615  {
616  //Directions A and B positive
617  Vmath::Vcopy(nquad0*nquad1,&(inarray[0]),1,&(o_tmp[0]),1);
618  }
619  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
620  {
621  //Direction A negative and B positive
622  for (int j=0; j<nquad1; j++)
623  {
624  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+j*nquad0,-1,&(o_tmp[0])+(j*nquad0),1);
625  }
626  }
627  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
628  {
629  //Direction A positive and B negative
630  for (int j=0; j<nquad1; j++)
631  {
632  Vmath::Vcopy(nquad0,&(inarray[0])+nquad0*(nquad1-1-j),1,&(o_tmp[0])+(j*nquad0),1);
633  }
634  }
635  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
636  {
637  //Direction A negative and B negative
638  for(int j=0; j<nquad1; j++)
639  {
640  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1-j*nquad0),-1,&(o_tmp[0])+(j*nquad0),1);
641  }
642  }
643  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
644  {
645  //Transposed, Direction A and B positive
646  for (int i=0; i<nquad0; i++)
647  {
648  Vmath::Vcopy(nquad1,&(inarray[0])+i,nquad0,&(o_tmp[0])+(i*nquad1),1);
649  }
650  }
651  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
652  {
653  //Transposed, Direction A positive and B negative
654  for (int i=0; i<nquad0; i++)
655  {
656  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0-1-i),nquad0,&(o_tmp[0])+(i*nquad1),1);
657  }
658  }
659  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
660  {
661  //Transposed, Direction A negative and B positive
662  for (int i=0; i<nquad0; i++)
663  {
664  Vmath::Vcopy(nquad1,&(inarray[0])+i+nquad0*(nquad1-1),-nquad0,&(o_tmp[0])+(i*nquad1),1);
665  }
666  }
667  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
668  {
669  //Transposed, Direction A and B negative
670  for (int i=0; i<nquad0; i++)
671  {
672  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0*nquad1-1-i),-nquad0,&(o_tmp[0])+(i*nquad1),1);
673  }
674  }
675  //interpolate
676  LibUtilities::Interp2D(m_base[0]->GetPointsKey(), m_base[1]->GetPointsKey(), o_tmp,
677  FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),outarray);
678  break;
679  case 1:
680  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
681  {
682  //Direction A and B positive
683  for (int k=0; k<nquad2; k++)
684  {
685  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1*k),1,&(o_tmp[0])+(k*nquad0),1);
686  }
687  }
688  else
689  {
690  //Direction A negative and B positive
691  for (int k=0; k<nquad2; k++)
692  {
693  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+(nquad0*nquad1*k),-1,&(o_tmp[0])+(k*nquad0),1);
694  }
695  }
696 
697  //interpolate
698  LibUtilities::Interp2D(m_base[0]->GetPointsKey(), m_base[2]->GetPointsKey(), o_tmp,
699  FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),outarray);
700  break;
701  case 2:
702  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
703  {
704  //Directions A and B positive
705  Vmath::Vcopy(nquad0*nquad2,&(inarray[0])+(nquad0-1),
706  nquad0,&(o_tmp[0]),1);
707  }
708  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
709  {
710  //Direction A negative and B positive
711  for (int k=0; k<nquad2; k++)
712  {
713  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1)+(k*nquad0*nquad1),
714  -nquad0,&(o_tmp[0])+(k*nquad0),1);
715  }
716  }
717  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
718  {
719  //Direction A positive and B negative
720  for (int k=0; k<nquad2; k++)
721  {
722  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0-1)+(nquad0*nquad1*(nquad2-1-k)),
723  nquad0,&(o_tmp[0])+(k*nquad0),1);
724  }
725  }
726  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
727  {
728  //Direction A negative and B negative
729  for (int k=0; k<nquad2; k++)
730  {
731  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1)+(nquad0*nquad1*(nquad2-1-k)),
732  -nquad0,&(o_tmp[0])+(k*nquad0),1);
733  }
734  }
735  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
736  {
737  //Transposed, Direction A and B positive
738  for (int j=0; j<nquad1; j++)
739  {
740  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0-1)+(j*nquad0),
741  nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
742  }
743  }
744  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
745  {
746  //Transposed, Direction A positive and B negative
747  for (int j=0; j<nquad0; j++)
748  {
749  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*nquad1-1-j*nquad0),
750  nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
751  }
752  }
753  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
754  {
755  //Transposed, Direction A negative and B positive
756  for (int j=0; j<nquad0; j++)
757  {
758  Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*nquad1*(nquad2-1)+nquad0+j*nquad0,
759  -nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
760  }
761  }
762  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
763  {
764  //Transposed, Direction A and B negative
765  for (int j=0; j<nquad0; j++)
766  {
767  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*nquad1*nquad2-1-j*nquad0),
768  -nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
769  }
770  }
771  //interpolate
772  LibUtilities::Interp2D(m_base[1]->GetPointsKey(), m_base[2]->GetPointsKey(), o_tmp,
773  FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),outarray);
774  break;
775  case 3:
776  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
777  {
778  //Directions A and B positive
779  for (int k=0; k<nquad2; k++)
780  {
781  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*(nquad1-1))+(k*nquad0*nquad1),
782  1,&(o_tmp[0])+(k*nquad0),1);
783  }
784  }
785  else
786  {
787  //Direction A negative and B positive
788  for (int k=0; k<nquad2; k++)
789  {
790  Vmath::Vcopy(nquad0,&(inarray[0])+(nquad0*nquad1-1)+(k*nquad0*nquad1),
791  -1,&(o_tmp[0])+(k*nquad0),1);
792  }
793  }
794  //interpolate
795  LibUtilities::Interp2D(m_base[0]->GetPointsKey(), m_base[2]->GetPointsKey(), o_tmp,
796  FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),outarray);
797 
798  break;
799  case 4:
800  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
801  {
802  //Directions A and B positive
803  Vmath::Vcopy(nquad1*nquad2,&(inarray[0]),nquad0,&(o_tmp[0]),1);
804  }
805  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
806  {
807  //Direction A negative and B positive
808  for (int k=0; k<nquad2; k++)
809  {
810  Vmath::Vcopy(nquad1,&(inarray[0])+nquad0*(nquad1-1)+(k*nquad0*nquad1),
811  -nquad0,&(o_tmp[0])+(k*nquad1),1);
812  }
813  }
814  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
815  {
816  //Direction A positive and B negative
817  for (int k=0; k<nquad2; k++)
818  {
819  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0*nquad1*(nquad2-1-k)),
820  nquad0,&(o_tmp[0])+(k*nquad1),1);
821  }
822  }
823  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
824  {
825  //Direction A negative and B negative
826  for (int k=0; k<nquad2; k++)
827  {
828  Vmath::Vcopy(nquad1,&(inarray[0])+nquad0*(nquad1-1)+(nquad0*nquad1*(nquad2-1-k)),
829  -nquad0,&(o_tmp[0])+(k*nquad1),1);
830  }
831  }
832  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
833  {
834  //Transposed, Direction A and B positive
835  for (int j=0; j<nquad1; j++)
836  {
837  Vmath::Vcopy(nquad2,&(inarray[0])+j*nquad0,nquad0*nquad1,
838  &(o_tmp[0])+(j*nquad2),1);
839  }
840  }
841  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
842  {
843  //Transposed, Direction A positive and B negative
844  for (int j=0; j<nquad1; j++)
845  {
846  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*(nquad1-1)-j*nquad0),
847  nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
848  }
849  }
850  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
851  {
852  //Transposed, Direction A negative and B positive
853  for (int j=0; j<nquad1; j++)
854  {
855  Vmath::Vcopy(nquad2,&(inarray[0])+nquad0*nquad1*(nquad2-1)+j*nquad0,
856  -nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
857  }
858  }
859  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
860  {
861  //Transposed, Direction A and B negative
862  for (int j=0; j<nquad1; j++)
863  {
864  Vmath::Vcopy(nquad2,&(inarray[0])+(nquad0*(nquad1*nquad2-1)-j*nquad0),
865  -nquad0*nquad1,&(o_tmp[0])+(j*nquad2),1);
866  }
867  }
868  //interpolate
869  LibUtilities::Interp2D(m_base[1]->GetPointsKey(), m_base[2]->GetPointsKey(), o_tmp,
870  FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),outarray);
871  break;
872  default:
873  ASSERTL0(false,"face value (> 4) is out of range");
874  break;
875  }
876  }
877 
878  /** \brief Get the normals along specficied face
879  * Get the face normals interplated to a points0 x points 0
880  * type distribution
881  **/
882  void PrismExp::v_ComputeFaceNormal(const int face)
883  {
884  const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
885  GetGeom()->GetMetricInfo();
886  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
887  SpatialDomains::GeomType type = geomFactors->GetGtype();
888  const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys);
889  const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
890 
891  int nq0 = ptsKeys[0].GetNumPoints();
892  int nq1 = ptsKeys[1].GetNumPoints();
893  int nq2 = ptsKeys[2].GetNumPoints();
894  int nq01 = nq0*nq1;
895  int nqtot;
896 
897 
898  LibUtilities::BasisKey tobasis0 = DetFaceBasisKey(face,0);
899  LibUtilities::BasisKey tobasis1 = DetFaceBasisKey(face,1);
900 
901  // Number of quadrature points in face expansion.
902  int nq_face = tobasis0.GetNumPoints()*tobasis1.GetNumPoints();
903 
904  int vCoordDim = GetCoordim();
905  int i;
906 
907  m_faceNormals[face] = Array<OneD, Array<OneD, NekDouble> >(vCoordDim);
908  Array<OneD, Array<OneD, NekDouble> > &normal = m_faceNormals[face];
909  for (i = 0; i < vCoordDim; ++i)
910  {
911  normal[i] = Array<OneD, NekDouble>(nq_face);
912  }
913 
914  // Regular geometry case
915  if (type == SpatialDomains::eRegular ||
916  type == SpatialDomains::eMovingRegular)
917  {
918  NekDouble fac;
919  // Set up normals
920  switch(face)
921  {
922  case 0:
923  {
924  for(i = 0; i < vCoordDim; ++i)
925  {
926  normal[i][0] = -df[3*i+2][0];;
927  }
928  break;
929  }
930  case 1:
931  {
932  for(i = 0; i < vCoordDim; ++i)
933  {
934  normal[i][0] = -df[3*i+1][0];
935  }
936  break;
937  }
938  case 2:
939  {
940  for(i = 0; i < vCoordDim; ++i)
941  {
942  normal[i][0] = df[3*i][0]+df[3*i+2][0];
943  }
944  break;
945  }
946  case 3:
947  {
948  for(i = 0; i < vCoordDim; ++i)
949  {
950  normal[i][0] = df[3*i+1][0];
951  }
952  break;
953  }
954  case 4:
955  {
956  for(i = 0; i < vCoordDim; ++i)
957  {
958  normal[i][0] = -df[3*i][0];
959  }
960  break;
961  }
962  default:
963  ASSERTL0(false,"face is out of range (face < 4)");
964  }
965 
966  // Normalise resulting vector.
967  fac = 0.0;
968  for(i = 0; i < vCoordDim; ++i)
969  {
970  fac += normal[i][0]*normal[i][0];
971  }
972  fac = 1.0/sqrt(fac);
973  for (i = 0; i < vCoordDim; ++i)
974  {
975  Vmath::Fill(nq_face,fac*normal[i][0],normal[i],1);
976  }
977  }
978  else
979  {
980  // Set up deformed normals.
981  int j, k;
982 
983  // Determine number of quadrature points on the face of 3D elmt
984  if (face == 0)
985  {
986  nqtot = nq0*nq1;
987  }
988  else if (face == 1 || face == 3)
989  {
990  nqtot = nq0*nq2;
991  }
992  else
993  {
994  nqtot = nq1*nq2;
995  }
996 
997  LibUtilities::PointsKey points0;
998  LibUtilities::PointsKey points1;
999 
1000  Array<OneD, NekDouble> faceJac(nqtot);
1001  Array<OneD, NekDouble> normals(vCoordDim*nqtot,0.0);
1002 
1003  // Extract Jacobian along face and recover local derivatives
1004  // (dx/dr) for polynomial interpolation by multiplying m_gmat by
1005  // jacobian
1006  switch(face)
1007  {
1008  case 0:
1009  {
1010  for(j = 0; j < nq01; ++j)
1011  {
1012  normals[j] = -df[2][j]*jac[j];
1013  normals[nqtot+j] = -df[5][j]*jac[j];
1014  normals[2*nqtot+j] = -df[8][j]*jac[j];
1015  faceJac[j] = jac[j];
1016  }
1017 
1018  points0 = ptsKeys[0];
1019  points1 = ptsKeys[1];
1020  break;
1021  }
1022 
1023  case 1:
1024  {
1025  for (j = 0; j < nq0; ++j)
1026  {
1027  for(k = 0; k < nq2; ++k)
1028  {
1029  int tmp = j+nq01*k;
1030  normals[j+k*nq0] =
1031  -df[1][tmp]*jac[tmp];
1032  normals[nqtot+j+k*nq0] =
1033  -df[4][tmp]*jac[tmp];
1034  normals[2*nqtot+j+k*nq0] =
1035  -df[7][tmp]*jac[tmp];
1036  faceJac[j+k*nq0] = jac[tmp];
1037  }
1038  }
1039 
1040  points0 = ptsKeys[0];
1041  points1 = ptsKeys[2];
1042  break;
1043  }
1044 
1045  case 2:
1046  {
1047  for (j = 0; j < nq1; ++j)
1048  {
1049  for(k = 0; k < nq2; ++k)
1050  {
1051  int tmp = nq0-1+nq0*j+nq01*k;
1052  normals[j+k*nq1] =
1053  (df[0][tmp]+df[2][tmp])*jac[tmp];
1054  normals[nqtot+j+k*nq1] =
1055  (df[3][tmp]+df[5][tmp])*jac[tmp];
1056  normals[2*nqtot+j+k*nq1] =
1057  (df[6][tmp]+df[8][tmp])*jac[tmp];
1058  faceJac[j+k*nq1] = jac[tmp];
1059  }
1060  }
1061 
1062  points0 = ptsKeys[1];
1063  points1 = ptsKeys[2];
1064  break;
1065  }
1066 
1067  case 3:
1068  {
1069  for (j = 0; j < nq0; ++j)
1070  {
1071  for(k = 0; k < nq2; ++k)
1072  {
1073  int tmp = nq0*(nq1-1) + j + nq01*k;
1074  normals[j+k*nq0] =
1075  df[1][tmp]*jac[tmp];
1076  normals[nqtot+j+k*nq0] =
1077  df[4][tmp]*jac[tmp];
1078  normals[2*nqtot+j+k*nq0] =
1079  df[7][tmp]*jac[tmp];
1080  faceJac[j+k*nq0] = jac[tmp];
1081  }
1082  }
1083 
1084  points0 = ptsKeys[0];
1085  points1 = ptsKeys[2];
1086  break;
1087  }
1088 
1089  case 4:
1090  {
1091  for (j = 0; j < nq1; ++j)
1092  {
1093  for(k = 0; k < nq2; ++k)
1094  {
1095  int tmp = j*nq0+nq01*k;
1096  normals[j+k*nq1] =
1097  -df[0][tmp]*jac[tmp];
1098  normals[nqtot+j+k*nq1] =
1099  -df[3][tmp]*jac[tmp];
1100  normals[2*nqtot+j+k*nq1] =
1101  -df[6][tmp]*jac[tmp];
1102  faceJac[j+k*nq1] = jac[tmp];
1103  }
1104  }
1105 
1106  points0 = ptsKeys[1];
1107  points1 = ptsKeys[2];
1108  break;
1109  }
1110 
1111  default:
1112  ASSERTL0(false,"face is out of range (face < 4)");
1113  }
1114 
1115 
1116  Array<OneD, NekDouble> work (nq_face, 0.0);
1117  // Interpolate Jacobian and invert
1118  LibUtilities::Interp2D(points0, points1, faceJac,
1119  tobasis0.GetPointsKey(),
1120  tobasis1.GetPointsKey(),
1121  work);
1122  Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
1123 
1124  // Interpolate normal and multiply by inverse Jacobian.
1125  for(i = 0; i < vCoordDim; ++i)
1126  {
1127  LibUtilities::Interp2D(points0, points1,
1128  &normals[i*nqtot],
1129  tobasis0.GetPointsKey(),
1130  tobasis1.GetPointsKey(),
1131  &normal[i][0]);
1132  Vmath::Vmul(nq_face,work,1,normal[i],1,normal[i],1);
1133  }
1134 
1135  // Normalise to obtain unit normals.
1136  Vmath::Zero(nq_face,work,1);
1137  for(i = 0; i < GetCoordim(); ++i)
1138  {
1139  Vmath::Vvtvp(nq_face,normal[i],1,normal[i],1,work,1,work,1);
1140  }
1141 
1142  Vmath::Vsqrt(nq_face,work,1,work,1);
1143  Vmath::Sdiv (nq_face,1.0,work,1,work,1);
1144 
1145  for(i = 0; i < GetCoordim(); ++i)
1146  {
1147  Vmath::Vmul(nq_face,normal[i],1,work,1,normal[i],1);
1148  }
1149  }
1150  }
1151 
1152  void PrismExp::v_MassMatrixOp(
1153  const Array<OneD, const NekDouble> &inarray,
1154  Array<OneD, NekDouble> &outarray,
1155  const StdRegions::StdMatrixKey &mkey)
1156  {
1157  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
1158  }
1159 
1160  void PrismExp::v_LaplacianMatrixOp(
1161  const Array<OneD, const NekDouble> &inarray,
1162  Array<OneD, NekDouble> &outarray,
1163  const StdRegions::StdMatrixKey &mkey)
1164  {
1165  PrismExp::LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1166  }
1167 
1168  void PrismExp::v_LaplacianMatrixOp(
1169  const int k1,
1170  const int k2,
1171  const Array<OneD, const NekDouble> &inarray,
1172  Array<OneD, NekDouble> &outarray,
1173  const StdRegions::StdMatrixKey &mkey)
1174  {
1175  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,mkey);
1176  }
1177 
1178  void PrismExp::v_HelmholtzMatrixOp(
1179  const Array<OneD, const NekDouble> &inarray,
1180  Array<OneD, NekDouble> &outarray,
1181  const StdRegions::StdMatrixKey &mkey)
1182  {
1183  PrismExp::v_HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1184  }
1185 
1186  void PrismExp::v_GeneralMatrixOp_MatOp(
1187  const Array<OneD, const NekDouble> &inarray,
1188  Array<OneD, NekDouble> &outarray,
1189  const StdRegions::StdMatrixKey &mkey)
1190  {
1191  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1192 
1193  if(inarray.get() == outarray.get())
1194  {
1195  Array<OneD,NekDouble> tmp(m_ncoeffs);
1196  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1197 
1198  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1199  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1200  }
1201  else
1202  {
1203  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1204  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1205  }
1206  }
1207 
1208  void PrismExp::v_SVVLaplacianFilter(
1209  Array<OneD, NekDouble> &array,
1210  const StdRegions::StdMatrixKey &mkey)
1211  {
1212  int nq = GetTotPoints();
1213 
1214  // Calculate sqrt of the Jacobian
1215  Array<OneD, const NekDouble> jac =
1216  m_metricinfo->GetJac(GetPointsKeys());
1217  Array<OneD, NekDouble> sqrt_jac(nq);
1218  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1219  {
1220  Vmath::Vsqrt(nq,jac,1,sqrt_jac,1);
1221  }
1222  else
1223  {
1224  Vmath::Fill(nq,sqrt(jac[0]),sqrt_jac,1);
1225  }
1226 
1227  // Multiply array by sqrt(Jac)
1228  Vmath::Vmul(nq,sqrt_jac,1,array,1,array,1);
1229 
1230  // Apply std region filter
1231  StdPrismExp::v_SVVLaplacianFilter( array, mkey);
1232 
1233  // Divide by sqrt(Jac)
1234  Vmath::Vdiv(nq,array,1,sqrt_jac,1,array,1);
1235  }
1236 
1237 
1238  //---------------------------------------
1239  // Matrix creation functions
1240  //---------------------------------------
1241 
1242  DNekMatSharedPtr PrismExp::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
1243  {
1244  DNekMatSharedPtr returnval;
1245 
1246  switch(mkey.GetMatrixType())
1247  {
1248  case StdRegions::eHybridDGHelmholtz:
1249  case StdRegions::eHybridDGLamToU:
1250  case StdRegions::eHybridDGLamToQ0:
1251  case StdRegions::eHybridDGLamToQ1:
1252  case StdRegions::eHybridDGLamToQ2:
1253  case StdRegions::eHybridDGHelmBndLam:
1254  case StdRegions::eInvLaplacianWithUnityMean:
1255  returnval = Expansion3D::v_GenMatrix(mkey);
1256  break;
1257  default:
1258  returnval = StdPrismExp::v_GenMatrix(mkey);
1259  break;
1260  }
1261 
1262  return returnval;
1263  }
1264 
1265  DNekMatSharedPtr PrismExp::v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
1266  {
1267  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1268  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1269  LibUtilities::BasisKey bkey2 = m_base[2]->GetBasisKey();
1270  StdRegions::StdPrismExpSharedPtr tmp =
1271  MemoryManager<StdPrismExp>::AllocateSharedPtr(bkey0, bkey1, bkey2);
1272 
1273  return tmp->GetStdMatrix(mkey);
1274  }
1275 
1276  DNekScalMatSharedPtr PrismExp::v_GetLocMatrix(const MatrixKey &mkey)
1277  {
1278  return m_matrixManager[mkey];
1279  }
1280 
1281  DNekScalBlkMatSharedPtr PrismExp::v_GetLocStaticCondMatrix(const MatrixKey &mkey)
1282  {
1283  return m_staticCondMatrixManager[mkey];
1284  }
1285 
1286  void PrismExp::v_DropLocStaticCondMatrix(const MatrixKey &mkey)
1287  {
1288  m_staticCondMatrixManager.DeleteObject(mkey);
1289  }
1290 
1291  DNekScalMatSharedPtr PrismExp::CreateMatrix(const MatrixKey &mkey)
1292  {
1293  DNekScalMatSharedPtr returnval;
1294  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1295 
1296  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,
1297  "Geometric information is not set up");
1298 
1299  switch(mkey.GetMatrixType())
1300  {
1301  case StdRegions::eMass:
1302  {
1303  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1304  {
1305  NekDouble one = 1.0;
1306  DNekMatSharedPtr mat = GenMatrix(mkey);
1307  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1308  }
1309  else
1310  {
1311  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1312  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1313  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1314  }
1315  break;
1316  }
1317 
1318  case StdRegions::eInvMass:
1319  {
1320  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1321  {
1322  NekDouble one = 1.0;
1323  StdRegions::StdMatrixKey masskey(StdRegions::eMass,DetShapeType(),*this);
1324  DNekMatSharedPtr mat = GenMatrix(masskey);
1325  mat->Invert();
1326 
1327  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1328  }
1329  else
1330  {
1331  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1332  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1333  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(fac,mat);
1334  }
1335  break;
1336  }
1337 
1338  case StdRegions::eWeakDeriv0:
1339  case StdRegions::eWeakDeriv1:
1340  case StdRegions::eWeakDeriv2:
1341  {
1342  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1343  {
1344  NekDouble one = 1.0;
1345  DNekMatSharedPtr mat = GenMatrix(mkey);
1346 
1347  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1348  }
1349  else
1350  {
1351  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1352  Array<TwoD, const NekDouble> df =
1353  m_metricinfo->GetDerivFactors(ptsKeys);
1354  int dir = 0;
1355 
1356  switch(mkey.GetMatrixType())
1357  {
1358  case StdRegions::eWeakDeriv0:
1359  dir = 0;
1360  break;
1361  case StdRegions::eWeakDeriv1:
1362  dir = 1;
1363  break;
1364  case StdRegions::eWeakDeriv2:
1365  dir = 2;
1366  break;
1367  default:
1368  break;
1369  }
1370 
1371  MatrixKey deriv0key(StdRegions::eWeakDeriv0,
1372  mkey.GetShapeType(), *this);
1373  MatrixKey deriv1key(StdRegions::eWeakDeriv1,
1374  mkey.GetShapeType(), *this);
1375  MatrixKey deriv2key(StdRegions::eWeakDeriv2,
1376  mkey.GetShapeType(), *this);
1377 
1378  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
1379  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
1380  DNekMat &deriv2 = *GetStdMatrix(deriv2key);
1381 
1382  int rows = deriv0.GetRows();
1383  int cols = deriv1.GetColumns();
1384 
1385  DNekMatSharedPtr WeakDeriv = MemoryManager<DNekMat>
1386  ::AllocateSharedPtr(rows,cols);
1387 
1388  (*WeakDeriv) = df[3*dir ][0]*deriv0
1389  + df[3*dir+1][0]*deriv1
1390  + df[3*dir+2][0]*deriv2;
1391 
1392  returnval = MemoryManager<DNekScalMat>
1393  ::AllocateSharedPtr(jac,WeakDeriv);
1394  }
1395  break;
1396  }
1397 
1398  case StdRegions::eLaplacian:
1399  {
1400  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1401  mkey.GetNVarCoeff() > 0 ||
1402  mkey.ConstFactorExists(
1403  StdRegions::eFactorSVVCutoffRatio))
1404  {
1405  NekDouble one = 1.0;
1406  DNekMatSharedPtr mat = GenMatrix(mkey);
1407  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1408  }
1409  else
1410  {
1411  MatrixKey lap00key(StdRegions::eLaplacian00,
1412  mkey.GetShapeType(), *this);
1413  MatrixKey lap01key(StdRegions::eLaplacian01,
1414  mkey.GetShapeType(), *this);
1415  MatrixKey lap02key(StdRegions::eLaplacian02,
1416  mkey.GetShapeType(), *this);
1417  MatrixKey lap11key(StdRegions::eLaplacian11,
1418  mkey.GetShapeType(), *this);
1419  MatrixKey lap12key(StdRegions::eLaplacian12,
1420  mkey.GetShapeType(), *this);
1421  MatrixKey lap22key(StdRegions::eLaplacian22,
1422  mkey.GetShapeType(), *this);
1423 
1424  DNekMat &lap00 = *GetStdMatrix(lap00key);
1425  DNekMat &lap01 = *GetStdMatrix(lap01key);
1426  DNekMat &lap02 = *GetStdMatrix(lap02key);
1427  DNekMat &lap11 = *GetStdMatrix(lap11key);
1428  DNekMat &lap12 = *GetStdMatrix(lap12key);
1429  DNekMat &lap22 = *GetStdMatrix(lap22key);
1430 
1431  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1432  Array<TwoD, const NekDouble> gmat
1433  = m_metricinfo->GetGmat(ptsKeys);
1434 
1435  int rows = lap00.GetRows();
1436  int cols = lap00.GetColumns();
1437 
1438  DNekMatSharedPtr lap = MemoryManager<DNekMat>
1439  ::AllocateSharedPtr(rows,cols);
1440 
1441  (*lap) = gmat[0][0]*lap00
1442  + gmat[4][0]*lap11
1443  + gmat[8][0]*lap22
1444  + gmat[3][0]*(lap01 + Transpose(lap01))
1445  + gmat[6][0]*(lap02 + Transpose(lap02))
1446  + gmat[7][0]*(lap12 + Transpose(lap12));
1447 
1448  returnval = MemoryManager<DNekScalMat>
1449  ::AllocateSharedPtr(jac,lap);
1450  }
1451  break;
1452  }
1453 
1454  case StdRegions::eHelmholtz:
1455  {
1456  NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorLambda);
1457  MatrixKey masskey(StdRegions::eMass,
1458  mkey.GetShapeType(), *this);
1459  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1460  MatrixKey lapkey(StdRegions::eLaplacian,
1461  mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1462  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1463 
1464  int rows = LapMat.GetRows();
1465  int cols = LapMat.GetColumns();
1466 
1467  DNekMatSharedPtr helm = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1468 
1469  NekDouble one = 1.0;
1470  (*helm) = LapMat + factor*MassMat;
1471 
1472  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
1473  break;
1474  }
1475  case StdRegions::eHybridDGHelmholtz:
1476  case StdRegions::eHybridDGLamToU:
1477  case StdRegions::eHybridDGLamToQ0:
1478  case StdRegions::eHybridDGLamToQ1:
1479  case StdRegions::eHybridDGHelmBndLam:
1480  case StdRegions::eInvLaplacianWithUnityMean:
1481  {
1482  NekDouble one = 1.0;
1483 
1484  DNekMatSharedPtr mat = GenMatrix(mkey);
1485  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1486 
1487  break;
1488  }
1489 
1490  case StdRegions::eInvHybridDGHelmholtz:
1491  {
1492  NekDouble one = 1.0;
1493 
1494  MatrixKey hkey(StdRegions::eHybridDGHelmholtz, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1495 // StdRegions::StdMatrixKey hkey(StdRegions::eHybridDGHelmholtz,
1496 // DetExpansionType(),*this,
1497 // mkey.GetConstant(0),
1498 // mkey.GetConstant(1));
1499  DNekMatSharedPtr mat = GenMatrix(hkey);
1500 
1501  mat->Invert();
1502  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1503  break;
1504  }
1505 
1506  case StdRegions::eIProductWRTBase:
1507  {
1508  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1509  {
1510  NekDouble one = 1.0;
1511  DNekMatSharedPtr mat = GenMatrix(mkey);
1512  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1513  }
1514  else
1515  {
1516  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1517  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1518  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1519  }
1520  break;
1521  }
1522  case StdRegions::ePreconLinearSpace:
1523  {
1524  NekDouble one = 1.0;
1525  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1526  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1527  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1528  DNekMatSharedPtr R=BuildVertexMatrix(A);
1529 
1530  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1531  }
1532  break;
1533  case StdRegions::ePreconLinearSpaceMass:
1534  {
1535  NekDouble one = 1.0;
1536  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1537  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1538  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1539  DNekMatSharedPtr R=BuildVertexMatrix(A);
1540 
1541  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1542  }
1543  break;
1544  case StdRegions::ePreconR:
1545  {
1546  NekDouble one = 1.0;
1547  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this,mkey.GetConstFactors(), mkey.GetVarCoeffs());
1548  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1549  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1550 
1551  DNekScalMatSharedPtr Atmp;
1552  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1553 
1554  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1555  }
1556  break;
1557  case StdRegions::ePreconRT:
1558  {
1559  NekDouble one = 1.0;
1560  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this,mkey.GetConstFactors(), mkey.GetVarCoeffs());
1561  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1562  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1563 
1564  DNekScalMatSharedPtr Atmp;
1565  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1566 
1567  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1568  }
1569  break;
1570  case StdRegions::ePreconRMass:
1571  {
1572  NekDouble one = 1.0;
1573  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1574  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1575  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1576 
1577  DNekScalMatSharedPtr Atmp;
1578  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1579 
1580  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1581  }
1582  break;
1583  case StdRegions::ePreconRTMass:
1584  {
1585  NekDouble one = 1.0;
1586  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1587  DNekScalBlkMatSharedPtr massStatCond = GetLocStaticCondMatrix(masskey);
1588  DNekScalMatSharedPtr A =massStatCond->GetBlock(0,0);
1589 
1590  DNekScalMatSharedPtr Atmp;
1591  DNekMatSharedPtr R=BuildTransformationMatrix(A,mkey.GetMatrixType());
1592 
1593  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1594  }
1595  break;
1596  default:
1597  {
1598  NekDouble one = 1.0;
1599  DNekMatSharedPtr mat = GenMatrix(mkey);
1600 
1601  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1602  }
1603  }
1604 
1605  return returnval;
1606  }
1607 
1608  DNekScalBlkMatSharedPtr PrismExp::CreateStaticCondMatrix(const MatrixKey &mkey)
1609  {
1610  DNekScalBlkMatSharedPtr returnval;
1611 
1612  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
1613 
1614  // set up block matrix system
1615  unsigned int nbdry = NumBndryCoeffs();
1616  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
1617  unsigned int exp_size[] = {nbdry, nint};
1618  unsigned int nblks=2;
1619  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks, nblks, exp_size, exp_size);
1620  NekDouble factor = 1.0;
1621 
1622  switch(mkey.GetMatrixType())
1623  {
1624  case StdRegions::eLaplacian:
1625  case StdRegions::eHelmholtz: // special case since Helmholtz not defined in StdRegions
1626  // use Deformed case for both regular and deformed geometries
1627  factor = 1.0;
1628  goto UseLocRegionsMatrix;
1629  break;
1630  default:
1631  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1632  {
1633  factor = 1.0;
1634  goto UseLocRegionsMatrix;
1635  }
1636  else
1637  {
1638  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1639  factor = mat->Scale();
1640  goto UseStdRegionsMatrix;
1641  }
1642  break;
1643  UseStdRegionsMatrix:
1644  {
1645  NekDouble invfactor = 1.0/factor;
1646  NekDouble one = 1.0;
1647  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
1648  DNekScalMatSharedPtr Atmp;
1649  DNekMatSharedPtr Asubmat;
1650 
1651  //TODO: check below
1652  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
1653  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
1654  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
1655  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
1656  }
1657  break;
1658  UseLocRegionsMatrix:
1659  {
1660  int i,j;
1661  NekDouble invfactor = 1.0/factor;
1662  NekDouble one = 1.0;
1663  DNekScalMat &mat = *GetLocMatrix(mkey);
1664  DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
1665  DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
1666  DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
1667  DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
1668 
1669  Array<OneD,unsigned int> bmap(nbdry);
1670  Array<OneD,unsigned int> imap(nint);
1671  GetBoundaryMap(bmap);
1672  GetInteriorMap(imap);
1673 
1674  for(i = 0; i < nbdry; ++i)
1675  {
1676  for(j = 0; j < nbdry; ++j)
1677  {
1678  (*A)(i,j) = mat(bmap[i],bmap[j]);
1679  }
1680 
1681  for(j = 0; j < nint; ++j)
1682  {
1683  (*B)(i,j) = mat(bmap[i],imap[j]);
1684  }
1685  }
1686 
1687  for(i = 0; i < nint; ++i)
1688  {
1689  for(j = 0; j < nbdry; ++j)
1690  {
1691  (*C)(i,j) = mat(imap[i],bmap[j]);
1692  }
1693 
1694  for(j = 0; j < nint; ++j)
1695  {
1696  (*D)(i,j) = mat(imap[i],imap[j]);
1697  }
1698  }
1699 
1700  // Calculate static condensed system
1701  if(nint)
1702  {
1703  D->Invert();
1704  (*B) = (*B)*(*D);
1705  (*A) = (*A) - (*B)*(*C);
1706  }
1707 
1708  DNekScalMatSharedPtr Atmp;
1709 
1710  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
1711  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
1712  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
1713  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
1714 
1715  }
1716  break;
1717  }
1718  return returnval;
1719  }
1720 
1721 
1722  /**
1723  * @brief Calculate the Laplacian multiplication in a matrix-free
1724  * manner.
1725  *
1726  * This function is the kernel of the Laplacian matrix-free operator,
1727  * and is used in #v_HelmholtzMatrixOp_MatFree to determine the effect
1728  * of the Helmholtz operator in a similar fashion.
1729  *
1730  * The majority of the calculation is precisely the same as in the
1731  * hexahedral expansion; however the collapsed co-ordinate system must
1732  * be taken into account when constructing the geometric factors. How
1733  * this is done is detailed more exactly in the tetrahedral expansion.
1734  * On entry to this function, the input #inarray must be in its
1735  * backwards-transformed state (i.e. \f$\mathbf{u} =
1736  * \mathbf{B}\hat{\mathbf{u}}\f$). The output is in coefficient space.
1737  *
1738  * @see %TetExp::v_HelmholtzMatrixOp_MatFree
1739  */
1740  void PrismExp::v_LaplacianMatrixOp_MatFree_Kernel(
1741  const Array<OneD, const NekDouble> &inarray,
1742  Array<OneD, NekDouble> &outarray,
1743  Array<OneD, NekDouble> &wsp)
1744  {
1745  int nquad0 = m_base[0]->GetNumPoints();
1746  int nquad1 = m_base[1]->GetNumPoints();
1747  int nquad2 = m_base[2]->GetNumPoints();
1748  int nqtot = nquad0*nquad1*nquad2;
1749  int i;
1750 
1751  // Set up temporary storage.
1752  Array<OneD,NekDouble> alloc(11*nqtot,0.0);
1753  Array<OneD,NekDouble> wsp1 (alloc ); // TensorDeriv 1
1754  Array<OneD,NekDouble> wsp2 (alloc+ 1*nqtot); // TensorDeriv 2
1755  Array<OneD,NekDouble> wsp3 (alloc+ 2*nqtot); // TensorDeriv 3
1756  Array<OneD,NekDouble> g0 (alloc+ 3*nqtot); // g0
1757  Array<OneD,NekDouble> g1 (alloc+ 4*nqtot); // g1
1758  Array<OneD,NekDouble> g2 (alloc+ 5*nqtot); // g2
1759  Array<OneD,NekDouble> g3 (alloc+ 6*nqtot); // g3
1760  Array<OneD,NekDouble> g4 (alloc+ 7*nqtot); // g4
1761  Array<OneD,NekDouble> g5 (alloc+ 8*nqtot); // g5
1762  Array<OneD,NekDouble> h0 (alloc+ 3*nqtot); // h0 == g0
1763  Array<OneD,NekDouble> h1 (alloc+ 6*nqtot); // h1 == g3
1764  Array<OneD,NekDouble> wsp4 (alloc+ 4*nqtot); // wsp4 == g1
1765  Array<OneD,NekDouble> wsp5 (alloc+ 5*nqtot); // wsp5 == g2
1766  Array<OneD,NekDouble> wsp6 (alloc+ 8*nqtot); // wsp6 == g5
1767  Array<OneD,NekDouble> wsp7 (alloc+ 3*nqtot); // wsp7 == g0
1768  Array<OneD,NekDouble> wsp8 (alloc+ 9*nqtot); // wsp8
1769  Array<OneD,NekDouble> wsp9 (alloc+10*nqtot); // wsp9
1770 
1771  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
1772  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
1773  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata();
1774  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
1775  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
1776  const Array<OneD, const NekDouble>& dbase2 = m_base[2]->GetDbdata();
1777 
1778  // Step 1. LAPLACIAN MATRIX OPERATION
1779  // wsp1 = du_dxi1 = D_xi1 * wsp0 = D_xi1 * u
1780  // wsp2 = du_dxi2 = D_xi2 * wsp0 = D_xi2 * u
1781  // wsp3 = du_dxi3 = D_xi3 * wsp0 = D_xi3 * u
1782  StdExpansion3D::PhysTensorDeriv(inarray,wsp1,wsp2,wsp3);
1783 
1784  const Array<TwoD, const NekDouble>& df =
1785  m_metricinfo->GetDerivFactors(GetPointsKeys());
1786  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
1787  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1788 
1789  // Step 2. Calculate the metric terms of the collapsed
1790  // coordinate transformation (Spencer's book P152)
1791  for (i = 0; i < nquad2; ++i)
1792  {
1793  Vmath::Fill(nquad0*nquad1, 2.0/(1.0-z2[i]), &h0[0]+i*nquad0*nquad1,1);
1794  Vmath::Fill(nquad0*nquad1, 2.0/(1.0-z2[i]), &h1[0]+i*nquad0*nquad1,1);
1795  }
1796  for (i = 0; i < nquad0; i++)
1797  {
1798  Blas::Dscal(nquad1*nquad2, 0.5*(1+z0[i]), &h1[0]+i, nquad0);
1799  }
1800 
1801  // Step 3. Construct combined metric terms for physical space to
1802  // collapsed coordinate system. Order of construction optimised
1803  // to minimise temporary storage
1804  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1805  {
1806  // wsp4 = d eta_1/d x_1
1807  Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &df[2][0], 1, &h1[0], 1, &wsp4[0], 1);
1808  // wsp5 = d eta_2/d x_1
1809  Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &df[5][0], 1, &h1[0], 1, &wsp5[0], 1);
1810  // wsp6 = d eta_3/d x_1d
1811  Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &df[8][0], 1, &h1[0], 1, &wsp6[0], 1);
1812 
1813  // g0 (overwrites h0)
1814  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0], 1, &g0[0], 1);
1815  Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1816 
1817  // g3 (overwrites h1)
1818  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &wsp4[0], 1, &df[4][0], 1, &wsp5[0], 1, &g3[0], 1);
1819  Vmath::Vvtvp (nqtot, &df[7][0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1820 
1821  // g4
1822  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0], 1, &g4[0], 1);
1823  Vmath::Vvtvp (nqtot, &df[8][0], 1, &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1824 
1825  // Overwrite wsp4/5/6 with g1/2/5
1826  // g1
1827  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &df[1][0], 1, &df[4][0], 1, &df[4][0], 1, &g1[0], 1);
1828  Vmath::Vvtvp (nqtot, &df[7][0], 1, &df[7][0], 1, &g1[0], 1, &g1[0], 1);
1829 
1830  // g2
1831  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &df[2][0], 1, &df[5][0], 1, &df[5][0], 1, &g2[0], 1);
1832  Vmath::Vvtvp (nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1833 
1834  // g5
1835  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &df[2][0], 1, &df[4][0], 1, &df[5][0], 1, &g5[0], 1);
1836  Vmath::Vvtvp (nqtot, &df[7][0], 1, &df[8][0], 1, &g5[0], 1, &g5[0], 1);
1837  }
1838  else
1839  {
1840  // wsp4 = d eta_1/d x_1
1841  Vmath::Svtsvtp(nqtot, df[0][0], &h0[0], 1, df[2][0], &h1[0], 1, &wsp4[0], 1);
1842  // wsp5 = d eta_2/d x_1
1843  Vmath::Svtsvtp(nqtot, df[3][0], &h0[0], 1, df[5][0], &h1[0], 1, &wsp5[0], 1);
1844  // wsp6 = d eta_3/d x_1
1845  Vmath::Svtsvtp(nqtot, df[6][0], &h0[0], 1, df[8][0], &h1[0], 1, &wsp6[0], 1);
1846 
1847  // g0 (overwrites h0)
1848  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0], 1, &g0[0], 1);
1849  Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1850 
1851  // g3 (overwrites h1)
1852  Vmath::Svtsvtp(nqtot, df[1][0], &wsp4[0], 1, df[4][0], &wsp5[0], 1, &g3[0], 1);
1853  Vmath::Svtvp (nqtot, df[7][0], &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1854 
1855  // g4
1856  Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1, &g4[0], 1);
1857  Vmath::Svtvp (nqtot, df[8][0], &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1858 
1859  // Overwrite wsp4/5/6 with g1/2/5
1860  // g1
1861  Vmath::Fill(nqtot, df[1][0]*df[1][0] + df[4][0]*df[4][0] + df[7][0]*df[7][0], &g1[0], 1);
1862 
1863  // g2
1864  Vmath::Fill(nqtot, df[2][0]*df[2][0] + df[5][0]*df[5][0] + df[8][0]*df[8][0], &g2[0], 1);
1865 
1866  // g5
1867  Vmath::Fill(nqtot, df[1][0]*df[2][0] + df[4][0]*df[5][0] + df[7][0]*df[8][0], &g5[0], 1);
1868  }
1869  // Compute component derivatives into wsp7, 8, 9 (wsp7 overwrites
1870  // g0).
1871  Vmath::Vvtvvtp(nqtot,&g0[0],1,&wsp1[0],1,&g3[0],1,&wsp2[0],1,&wsp7[0],1);
1872  Vmath::Vvtvp (nqtot,&g4[0],1,&wsp3[0],1,&wsp7[0],1,&wsp7[0],1);
1873  Vmath::Vvtvvtp(nqtot,&g1[0],1,&wsp2[0],1,&g3[0],1,&wsp1[0],1,&wsp8[0],1);
1874  Vmath::Vvtvp (nqtot,&g5[0],1,&wsp3[0],1,&wsp8[0],1,&wsp8[0],1);
1875  Vmath::Vvtvvtp(nqtot,&g2[0],1,&wsp3[0],1,&g4[0],1,&wsp1[0],1,&wsp9[0],1);
1876  Vmath::Vvtvp (nqtot,&g5[0],1,&wsp2[0],1,&wsp9[0],1,&wsp9[0],1);
1877 
1878  // Step 4.
1879  // Multiply by quadrature metric
1880  MultiplyByQuadratureMetric(wsp7,wsp7);
1881  MultiplyByQuadratureMetric(wsp8,wsp8);
1882  MultiplyByQuadratureMetric(wsp9,wsp9);
1883 
1884  // Perform inner product w.r.t derivative bases.
1885  IProductWRTBase_SumFacKernel(dbase0,base1,base2,wsp7,wsp1, wsp,false,true,true);
1886  IProductWRTBase_SumFacKernel(base0,dbase1,base2,wsp8,wsp2, wsp,true,false,true);
1887  IProductWRTBase_SumFacKernel(base0,base1,dbase2,wsp9,outarray,wsp,true,true,false);
1888 
1889  // Step 5.
1890  // Sum contributions from wsp1, wsp2 and outarray.
1891  Vmath::Vadd(m_ncoeffs,wsp1.get(),1,outarray.get(),1,outarray.get(),1);
1892  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
1893  }
1894 
1895 
1896  void PrismExp::v_GetSimplexEquiSpacedConnectivity(
1897  Array<OneD, int> &conn,
1898  bool oldstandard)
1899  {
1900  int np0 = m_base[0]->GetNumPoints();
1901  int np1 = m_base[1]->GetNumPoints();
1902  int np2 = m_base[2]->GetNumPoints();
1903  int np = max(np0,max(np1,np2));
1904  Array<OneD, int> prismpt(6);
1905  bool standard = true;
1906 
1907  int vid0 = m_geom->GetVid(0);
1908  int vid1 = m_geom->GetVid(1);
1909  int vid2 = m_geom->GetVid(4);
1910  int rotate = 0;
1911 
1912  // sort out prism rotation according to
1913  if((vid2 < vid1)&&(vid2 < vid0)) // top triangle vertex is lowest id
1914  {
1915  rotate = 0;
1916  if(vid0 > vid1)
1917  {
1918  standard = false;// reverse base direction
1919  }
1920  }
1921  else if((vid1 < vid2)&&(vid1 < vid0))
1922  {
1923  rotate = 1;
1924  if(vid2 > vid0)
1925  {
1926  standard = false;// reverse base direction
1927  }
1928  }
1929  else if ((vid0 < vid2)&&(vid0 < vid1))
1930  {
1931  rotate = 2;
1932  if(vid1 > vid2)
1933  {
1934  standard = false; // reverse base direction
1935  }
1936  }
1937 
1938  conn = Array<OneD, int>(12*(np-1)*(np-1)*(np-1));
1939 
1940  int row = 0;
1941  int rowp1 = 0;
1942  int plane = 0;
1943  int row1 = 0;
1944  int row1p1 = 0;
1945  int planep1 = 0;
1946  int cnt = 0;
1947 
1948 
1949  Array<OneD, int> rot(3);
1950 
1951  rot[0] = (0+rotate)%3;
1952  rot[1] = (1+rotate)%3;
1953  rot[2] = (2+rotate)%3;
1954 
1955  // lower diagonal along 1-3 on base
1956  for(int i = 0; i < np-1; ++i)
1957  {
1958  planep1 += (np-i)*np;
1959  row = 0; // current plane row offset
1960  rowp1 = 0; // current plane row plus one offset
1961  row1 = 0; // next plane row offset
1962  row1p1 = 0; // nex plane row plus one offset
1963  if(standard == false)
1964  {
1965  for(int j = 0; j < np-1; ++j)
1966  {
1967  rowp1 += np-i;
1968  row1p1 += np-i-1;
1969  for(int k = 0; k < np-i-2; ++k)
1970  {
1971  // bottom prism block
1972  prismpt[rot[0]] = plane + row + k;
1973  prismpt[rot[1]] = plane + row + k+1;
1974  prismpt[rot[2]] = planep1 + row1 + k;
1975 
1976  prismpt[3+rot[0]] = plane + rowp1 + k;
1977  prismpt[3+rot[1]] = plane + rowp1 + k+1;
1978  prismpt[3+rot[2]] = planep1 + row1p1 + k;
1979 
1980  conn[cnt++] = prismpt[0];
1981  conn[cnt++] = prismpt[1];
1982  conn[cnt++] = prismpt[3];
1983  conn[cnt++] = prismpt[2];
1984 
1985  conn[cnt++] = prismpt[5];
1986  conn[cnt++] = prismpt[2];
1987  conn[cnt++] = prismpt[3];
1988  conn[cnt++] = prismpt[4];
1989 
1990  conn[cnt++] = prismpt[3];
1991  conn[cnt++] = prismpt[1];
1992  conn[cnt++] = prismpt[4];
1993  conn[cnt++] = prismpt[2];
1994 
1995  // upper prism block.
1996  prismpt[rot[0]] = planep1 + row1 + k+1;
1997  prismpt[rot[1]] = planep1 + row1 + k;
1998  prismpt[rot[2]] = plane + row + k+1;
1999 
2000  prismpt[3+rot[0]] = planep1 + row1p1 + k+1;
2001  prismpt[3+rot[1]] = planep1 + row1p1 + k;
2002  prismpt[3+rot[2]] = plane + rowp1 + k+1;
2003 
2004 
2005  conn[cnt++] = prismpt[0];
2006  conn[cnt++] = prismpt[1];
2007  conn[cnt++] = prismpt[2];
2008  conn[cnt++] = prismpt[5];
2009 
2010  conn[cnt++] = prismpt[5];
2011  conn[cnt++] = prismpt[0];
2012  conn[cnt++] = prismpt[4];
2013  conn[cnt++] = prismpt[1];
2014 
2015  conn[cnt++] = prismpt[3];
2016  conn[cnt++] = prismpt[4];
2017  conn[cnt++] = prismpt[0];
2018  conn[cnt++] = prismpt[5];
2019 
2020  }
2021 
2022  // bottom prism block
2023  prismpt[rot[0]] = plane + row + np-i-2;
2024  prismpt[rot[1]] = plane + row + np-i-1;
2025  prismpt[rot[2]] = planep1 + row1 + np-i-2;
2026 
2027  prismpt[3+rot[0]] = plane + rowp1 + np-i-2;
2028  prismpt[3+rot[1]] = plane + rowp1 + np-i-1;
2029  prismpt[3+rot[2]] = planep1 + row1p1 + np-i-2;
2030 
2031  conn[cnt++] = prismpt[0];
2032  conn[cnt++] = prismpt[1];
2033  conn[cnt++] = prismpt[3];
2034  conn[cnt++] = prismpt[2];
2035 
2036  conn[cnt++] = prismpt[5];
2037  conn[cnt++] = prismpt[2];
2038  conn[cnt++] = prismpt[3];
2039  conn[cnt++] = prismpt[4];
2040 
2041  conn[cnt++] = prismpt[3];
2042  conn[cnt++] = prismpt[1];
2043  conn[cnt++] = prismpt[4];
2044  conn[cnt++] = prismpt[2];
2045 
2046  row += np-i;
2047  row1 += np-i-1;
2048  }
2049 
2050  }
2051  else
2052  { // lower diagonal along 0-4 on base
2053  for(int j = 0; j < np-1; ++j)
2054  {
2055  rowp1 += np-i;
2056  row1p1 += np-i-1;
2057  for(int k = 0; k < np-i-2; ++k)
2058  {
2059  // bottom prism block
2060  prismpt[rot[0]] = plane + row + k;
2061  prismpt[rot[1]] = plane + row + k+1;
2062  prismpt[rot[2]] = planep1 + row1 + k;
2063 
2064  prismpt[3+rot[0]] = plane + rowp1 + k;
2065  prismpt[3+rot[1]] = plane + rowp1 + k+1;
2066  prismpt[3+rot[2]] = planep1 + row1p1 + k;
2067 
2068  conn[cnt++] = prismpt[0];
2069  conn[cnt++] = prismpt[1];
2070  conn[cnt++] = prismpt[4];
2071  conn[cnt++] = prismpt[2];
2072 
2073  conn[cnt++] = prismpt[4];
2074  conn[cnt++] = prismpt[3];
2075  conn[cnt++] = prismpt[0];
2076  conn[cnt++] = prismpt[2];
2077 
2078  conn[cnt++] = prismpt[3];
2079  conn[cnt++] = prismpt[4];
2080  conn[cnt++] = prismpt[5];
2081  conn[cnt++] = prismpt[2];
2082 
2083  // upper prism block.
2084  prismpt[rot[0]] = planep1 + row1 + k+1;
2085  prismpt[rot[1]] = planep1 + row1 + k;
2086  prismpt[rot[2]] = plane + row + k+1;
2087 
2088  prismpt[3+rot[0]] = planep1 + row1p1 + k+1;
2089  prismpt[3+rot[1]] = planep1 + row1p1 + k;
2090  prismpt[3+rot[2]] = plane + rowp1 + k+1;
2091 
2092  conn[cnt++] = prismpt[0];
2093  conn[cnt++] = prismpt[2];
2094  conn[cnt++] = prismpt[1];
2095  conn[cnt++] = prismpt[5];
2096 
2097  conn[cnt++] = prismpt[3];
2098  conn[cnt++] = prismpt[5];
2099  conn[cnt++] = prismpt[0];
2100  conn[cnt++] = prismpt[1];
2101 
2102  conn[cnt++] = prismpt[5];
2103  conn[cnt++] = prismpt[3];
2104  conn[cnt++] = prismpt[4];
2105  conn[cnt++] = prismpt[1];
2106  }
2107 
2108  // bottom prism block
2109  prismpt[rot[0]] = plane + row + np-i-2;
2110  prismpt[rot[1]] = plane + row + np-i-1;
2111  prismpt[rot[2]] = planep1 + row1 + np-i-2;
2112 
2113  prismpt[3+rot[0]] = plane + rowp1 + np-i-2;
2114  prismpt[3+rot[1]] = plane + rowp1 + np-i-1;
2115  prismpt[3+rot[2]] = planep1 + row1p1 + np-i-2;
2116 
2117  conn[cnt++] = prismpt[0];
2118  conn[cnt++] = prismpt[1];
2119  conn[cnt++] = prismpt[4];
2120  conn[cnt++] = prismpt[2];
2121 
2122  conn[cnt++] = prismpt[4];
2123  conn[cnt++] = prismpt[3];
2124  conn[cnt++] = prismpt[0];
2125  conn[cnt++] = prismpt[2];
2126 
2127  conn[cnt++] = prismpt[3];
2128  conn[cnt++] = prismpt[4];
2129  conn[cnt++] = prismpt[5];
2130  conn[cnt++] = prismpt[2];
2131 
2132  row += np-i;
2133  row1 += np-i-1;
2134  }
2135 
2136  }
2137  plane += (np-i)*np;
2138  }
2139  }
2140 
2141  }//end of namespace
2142 }//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:1178
Nektar::StdRegions::StdExpansion::GetPointsKeys
const LibUtilities::PointsKeyVector GetPointsKeys() const
Definition: StdExpansion.h:1276
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:454
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
Nektar::LocalRegions::PrismExp::v_GetFacePhysVals
virtual void v_GetFacePhysVals(const int face, const StdRegions::StdExpansionSharedPtr &FaceExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient)
Definition: PrismExp.cpp:593
Nektar::StdRegions::StdExpansion::GenMatrix
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:1009
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::StdRegions::StdExpansion::GetFaceNumPoints
int GetFaceNumPoints(const int i) const
This function returns the number of quadrature points belonging to the i-th face. ...
Definition: StdExpansion.h:339
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:215
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:902
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:1186
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:1423
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:1896
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:1281
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:1152
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::GetForient
StdRegions::Orientation GetForient(int face)
Definition: StdExpansion.h:731
Nektar::StdRegions::StdExpansion::GetLocStaticCondMatrix
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
Definition: StdExpansion.h:721
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:684
Nektar::DNekScalMatSharedPtr
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:73
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:1291
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:1265
Nektar::StdRegions::StdExpansion::GetStdStaticCondMatrix
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:689
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:1740
Nektar::LocalRegions::PrismExp::m_staticCondMatrixManager
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: PrismExp.h:217
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:1242
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:1208
Nektar::LocalRegions::PrismExp::v_DropLocStaticCondMatrix
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: PrismExp.cpp:1286
Nektar::LocalRegions::PrismExp::v_GetTracePhysVals
virtual void v_GetTracePhysVals(const int face, const StdRegions::StdExpansionSharedPtr &FaceExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient)
Returns the physical values at the quadrature points of a face.
Definition: PrismExp.cpp:583
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:882
Nektar::StdRegions::StdExpansion::GetInteriorMap
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:795
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:43
Interp.h
Nektar::StdRegions::StdExpansion3D::m_faceNormals
std::map< int, NormalVector > m_faceNormals
Definition: StdExpansion3D.h:196
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:1276
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:380
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:1791
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:1356
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
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:1160
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:1608
Nektar::StdRegions::StdExpansion::GetBoundaryMap
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:790
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 Tue Dec 8 2015 23:45:29 for Nektar++ by   doxygen 1.8.9.1