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