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