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