Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PyrExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File PyrExp.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: PyrExp routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <LocalRegions/PyrExp.h>
37 #include <LibUtilities/Foundations/Interp.h>
38 
39 namespace Nektar
40 {
41  namespace LocalRegions
42  {
43 
44  PyrExp::PyrExp(const LibUtilities::BasisKey &Ba,
45  const LibUtilities::BasisKey &Bb,
46  const LibUtilities::BasisKey &Bc,
47  const SpatialDomains::PyrGeomSharedPtr &geom):
48  StdExpansion (LibUtilities::StdPyrData::getNumberOfCoefficients(
49  Ba.GetNumModes(),
50  Bb.GetNumModes(),
51  Bc.GetNumModes()),
52  3, Ba, Bb, Bc),
53  StdExpansion3D(LibUtilities::StdPyrData::getNumberOfCoefficients(
54  Ba.GetNumModes(),
55  Bb.GetNumModes(),
56  Bc.GetNumModes()),
57  Ba, Bb, Bc),
58  StdPyrExp (Ba,Bb,Bc),
59  Expansion (geom),
60  Expansion3D (geom),
61  m_matrixManager(
62  boost::bind(&PyrExp::CreateMatrix, this, _1),
63  std::string("PyrExpMatrix")),
64  m_staticCondMatrixManager(
65  boost::bind(&PyrExp::CreateStaticCondMatrix, this, _1),
66  std::string("PyrExpStaticCondMatrix"))
67  {
68  }
69 
70  PyrExp::PyrExp(const PyrExp &T):
71  StdExpansion (T),
72  StdExpansion3D(T),
73  StdPyrExp (T),
74  Expansion (T),
75  Expansion3D (T),
76  m_matrixManager(T.m_matrixManager),
77  m_staticCondMatrixManager(T.m_staticCondMatrixManager)
78  {
79  }
80 
81  PyrExp::~PyrExp()
82  {
83  }
84 
85 
86  //----------------------------
87  // Integration Methods
88  //----------------------------
89 
90  /**
91  * \brief Integrate the physical point list \a inarray over pyramidic
92  * region and return the value.
93  *
94  * Inputs:\n
95  *
96  * - \a inarray: definition of function to be returned at quadrature
97  * point of expansion.
98  *
99  * Outputs:\n
100  *
101  * - returns \f$\int^1_{-1}\int^1_{-1}\int^1_{-1} u(\bar \eta_1,
102  * \eta_2, \eta_3) J[i,j,k] d \bar \eta_1 d \eta_2 d \eta_3\f$ \n \f$=
103  * \sum_{i=0}^{Q_1 - 1} \sum_{j=0}^{Q_2 - 1} \sum_{k=0}^{Q_3 - 1}
104  * u(\bar \eta_{1i}^{0,0}, \eta_{2j}^{0,0},\eta_{3k}^{2,0})w_{i}^{0,0}
105  * w_{j}^{0,0} \hat w_{k}^{2,0} \f$ \n where \f$inarray[i,j, k] =
106  * u(\bar \eta_{1i},\eta_{2j}, \eta_{3k}) \f$, \n \f$\hat w_{k}^{2,0}
107  * = \frac {w^{2,0}} {2} \f$ \n and \f$ J[i,j,k] \f$ is the Jacobian
108  * evaluated at the quadrature point.
109  */
110  NekDouble PyrExp::v_Integral(const Array<OneD, const NekDouble> &inarray)
111  {
112  int nquad0 = m_base[0]->GetNumPoints();
113  int nquad1 = m_base[1]->GetNumPoints();
114  int nquad2 = m_base[2]->GetNumPoints();
115  Array<OneD, const NekDouble> jac = m_metricinfo->GetJac(GetPointsKeys());
116  Array<OneD, NekDouble> tmp(nquad0*nquad1*nquad2);
117 
118  // multiply inarray with Jacobian
119  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
120  {
121  Vmath::Vmul(nquad0*nquad1*nquad2,&jac[0],1,(NekDouble*)&inarray[0],1, &tmp[0],1);
122  }
123  else
124  {
125  Vmath::Smul(nquad0*nquad1*nquad2,(NekDouble) jac[0], (NekDouble*)&inarray[0],1,&tmp[0],1);
126  }
127 
128  // call StdPyrExp version;
129  return StdPyrExp::v_Integral(tmp);
130  }
131 
132 
133  //----------------------------
134  // Differentiation Methods
135  //----------------------------
136 
137  void PyrExp::v_PhysDeriv(const Array<OneD, const NekDouble>& inarray,
138  Array<OneD, NekDouble>& out_d0,
139  Array<OneD, NekDouble>& out_d1,
140  Array<OneD, NekDouble>& out_d2)
141  {
142  int nquad0 = m_base[0]->GetNumPoints();
143  int nquad1 = m_base[1]->GetNumPoints();
144  int nquad2 = m_base[2]->GetNumPoints();
145  Array<TwoD, const NekDouble> gmat =
146  m_metricinfo->GetDerivFactors(GetPointsKeys());
147  Array<OneD,NekDouble> diff0(nquad0*nquad1*nquad2);
148  Array<OneD,NekDouble> diff1(nquad0*nquad1*nquad2);
149  Array<OneD,NekDouble> diff2(nquad0*nquad1*nquad2);
150 
151  StdPyrExp::v_PhysDeriv(inarray, diff0, diff1, diff2);
152 
153  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
154  {
155  if(out_d0.num_elements())
156  {
157  Vmath::Vmul (nquad0*nquad1*nquad2,&gmat[0][0],1,&diff0[0],1, &out_d0[0], 1);
158  Vmath::Vvtvp (nquad0*nquad1*nquad2,&gmat[1][0],1,&diff1[0],1, &out_d0[0], 1,&out_d0[0],1);
159  Vmath::Vvtvp (nquad0*nquad1*nquad2,&gmat[2][0],1,&diff2[0],1, &out_d0[0], 1,&out_d0[0],1);
160  }
161 
162  if(out_d1.num_elements())
163  {
164  Vmath::Vmul (nquad0*nquad1*nquad2,&gmat[3][0],1,&diff0[0],1, &out_d1[0], 1);
165  Vmath::Vvtvp (nquad0*nquad1*nquad2,&gmat[4][0],1,&diff1[0],1, &out_d1[0], 1,&out_d1[0],1);
166  Vmath::Vvtvp (nquad0*nquad1*nquad2,&gmat[5][0],1,&diff2[0],1, &out_d1[0], 1,&out_d1[0],1);
167  }
168 
169  if(out_d2.num_elements())
170  {
171  Vmath::Vmul (nquad0*nquad1*nquad2,&gmat[6][0],1,&diff0[0],1, &out_d2[0], 1);
172  Vmath::Vvtvp (nquad0*nquad1*nquad2,&gmat[7][0],1,&diff1[0],1, &out_d2[0], 1, &out_d2[0],1);
173  Vmath::Vvtvp (nquad0*nquad1*nquad2,&gmat[8][0],1,&diff2[0],1, &out_d2[0], 1, &out_d2[0],1);
174  }
175  }
176  else // regular geometry
177  {
178  if(out_d0.num_elements())
179  {
180  Vmath::Smul (nquad0*nquad1*nquad2,gmat[0][0],&diff0[0],1, &out_d0[0], 1);
181  Blas::Daxpy (nquad0*nquad1*nquad2,gmat[1][0],&diff1[0],1, &out_d0[0], 1);
182  Blas::Daxpy (nquad0*nquad1*nquad2,gmat[2][0],&diff2[0],1, &out_d0[0], 1);
183  }
184 
185  if(out_d1.num_elements())
186  {
187  Vmath::Smul (nquad0*nquad1*nquad2,gmat[3][0],&diff0[0],1, &out_d1[0], 1);
188  Blas::Daxpy (nquad0*nquad1*nquad2,gmat[4][0],&diff1[0],1, &out_d1[0], 1);
189  Blas::Daxpy (nquad0*nquad1*nquad2,gmat[5][0],&diff2[0],1, &out_d1[0], 1);
190  }
191 
192  if(out_d2.num_elements())
193  {
194  Vmath::Smul (nquad0*nquad1*nquad2,gmat[6][0],&diff0[0],1, &out_d2[0], 1);
195  Blas::Daxpy (nquad0*nquad1*nquad2,gmat[7][0],&diff1[0],1, &out_d2[0], 1);
196  Blas::Daxpy (nquad0*nquad1*nquad2,gmat[8][0],&diff2[0],1, &out_d2[0], 1);
197  }
198  }
199  }
200 
201 
202  //---------------------------------------
203  // Transforms
204  //---------------------------------------
205 
206  /**
207  * \brief Forward transform from physical quadrature space stored in
208  * \a inarray and evaluate the expansion coefficients and store in \a
209  * (this)->m_coeffs
210  *
211  * Inputs:\n
212  *
213  * - \a inarray: array of physical quadrature points to be transformed
214  *
215  * Outputs:\n
216  *
217  * - (this)->_coeffs: updated array of expansion coefficients.
218  */
219  void PyrExp::v_FwdTrans(const Array<OneD, const NekDouble>& inarray,
220  Array<OneD, NekDouble>& outarray)
221  {
222  if(m_base[0]->Collocation() &&
223  m_base[1]->Collocation() &&
224  m_base[2]->Collocation())
225  {
226  Vmath::Vcopy(GetNcoeffs(),&inarray[0],1,&outarray[0],1);
227  }
228  else
229  {
230  v_IProductWRTBase(inarray,outarray);
231 
232  // get Mass matrix inverse
233  MatrixKey masskey(StdRegions::eInvMass,
234  DetShapeType(),*this);
235  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
236 
237  // copy inarray in case inarray == outarray
238  DNekVec in (m_ncoeffs,outarray);
239  DNekVec out(m_ncoeffs,outarray,eWrapper);
240 
241  out = (*matsys)*in;
242  }
243  }
244 
245 
246  //---------------------------------------
247  // Inner product functions
248  //---------------------------------------
249 
250  /**
251  * \brief Calculate the inner product of inarray with respect to the
252  * basis B=base0*base1*base2 and put into outarray:
253  *
254  * \f$ \begin{array}{rcl} I_{pqr} = (\phi_{pqr}, u)_{\delta} & = &
255  * \sum_{i=0}^{nq_0} \sum_{j=0}^{nq_1} \sum_{k=0}^{nq_2} \psi_{p}^{a}
256  * (\bar \eta_{1i}) \psi_{q}^{a} (\eta_{2j}) \psi_{pqr}^{c}
257  * (\eta_{3k}) w_i w_j w_k u(\bar \eta_{1,i} \eta_{2,j} \eta_{3,k})
258  * J_{i,j,k}\\ & = & \sum_{i=0}^{nq_0} \psi_p^a(\bar \eta_{1,i})
259  * \sum_{j=0}^{nq_1} \psi_{q}^a(\eta_{2,j}) \sum_{k=0}^{nq_2}
260  * \psi_{pqr}^c u(\bar \eta_{1i},\eta_{2j},\eta_{3k}) J_{i,j,k}
261  * \end{array} \f$ \n
262  *
263  * where
264  *
265  * \f$\phi_{pqr} (\xi_1 , \xi_2 , \xi_3) = \psi_p^a (\bar \eta_1)
266  * \psi_{q}^a (\eta_2) \psi_{pqr}^c (\eta_3) \f$ \n
267  *
268  * which can be implemented as \n \f$f_{pqr} (\xi_{3k}) =
269  * \sum_{k=0}^{nq_3} \psi_{pqr}^c u(\bar
270  * \eta_{1i},\eta_{2j},\eta_{3k}) J_{i,j,k} = {\bf B_3 U} \f$ \n \f$
271  * g_{pq} (\xi_{3k}) = \sum_{j=0}^{nq_1} \psi_{q}^a (\xi_{2j}) f_{pqr}
272  * (\xi_{3k}) = {\bf B_2 F} \f$ \n \f$ (\phi_{pqr}, u)_{\delta} =
273  * \sum_{k=0}^{nq_0} \psi_{p}^a (\xi_{3k}) g_{pq} (\xi_{3k}) = {\bf
274  * B_1 G} \f$
275  */
276  void PyrExp::v_IProductWRTBase(
277  const Array<OneD, const NekDouble> &inarray,
278  Array<OneD, NekDouble> &outarray)
279  {
280  int nquad0 = m_base[0]->GetNumPoints();
281  int nquad1 = m_base[1]->GetNumPoints();
282  int nquad2 = m_base[2]->GetNumPoints();
283  Array<OneD, const NekDouble> jac = m_metricinfo->GetJac(GetPointsKeys());
284  Array<OneD, NekDouble> tmp(nquad0*nquad1*nquad2);
285 
286  // multiply inarray with Jacobian
287  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
288  {
289  Vmath::Vmul(nquad0*nquad1*nquad2,&jac[0],1,(NekDouble*)&inarray[0],1,&tmp[0],1);
290  }
291  else
292  {
293  Vmath::Smul(nquad0*nquad1*nquad2,jac[0],(NekDouble*)&inarray[0],1,&tmp[0],1);
294  }
295 
296  StdPyrExp::v_IProductWRTBase(tmp,outarray);
297  }
298 
299 
300  //---------------------------------------
301  // Evaluation functions
302  //---------------------------------------
303 
304  /*
305  * @brief Get the coordinates #coords at the local coordinates
306  * #Lcoords
307  */
308  void PyrExp::v_GetCoord(const Array<OneD, const NekDouble>& Lcoords,
309  Array<OneD, NekDouble>& coords)
310  {
311  int i;
312 
313  ASSERTL1(Lcoords[0] <= -1.0 && Lcoords[0] >= 1.0 &&
314  Lcoords[1] <= -1.0 && Lcoords[1] >= 1.0 &&
315  Lcoords[2] <= -1.0 && Lcoords[2] >= 1.0,
316  "Local coordinates are not in region [-1,1]");
317 
318  // m_geom->FillGeom(); // TODO: implement FillGeom()
319 
320  for(i = 0; i < m_geom->GetCoordim(); ++i)
321  {
322  coords[i] = m_geom->GetCoord(i,Lcoords);
323  }
324  }
325 
326  void PyrExp::v_GetCoords(
327  Array<OneD, NekDouble> &coords_1,
328  Array<OneD, NekDouble> &coords_2,
329  Array<OneD, NekDouble> &coords_3)
330  {
331  Expansion::v_GetCoords(coords_1, coords_2, coords_3);
332  }
333 
334  NekDouble PyrExp::v_PhysEvaluate(const Array<OneD, const NekDouble>& coord,
335  const Array<OneD, const NekDouble>& physvals)
336  {
337  Array<OneD,NekDouble> Lcoord(3);
338 
339  ASSERTL0(m_geom,"m_geom not defined");
340 
341  //TODO: check GetLocCoords()
342  m_geom->GetLocCoords(coord, Lcoord);
343 
344  return StdPyrExp::v_PhysEvaluate(Lcoord, physvals);
345  }
346 
347 
348  //---------------------------------------
349  // Helper functions
350  //---------------------------------------
351 
352  int PyrExp::v_GetCoordim()
353  {
354  return m_geom->GetCoordim();
355  }
356 
357  StdRegions::Orientation PyrExp::v_GetFaceOrient(int face)
358  {
359  return GetGeom3D()->GetFaceOrient(face);
360  }
361 
362  void PyrExp::v_GetFacePhysVals(
363  const int face,
364  const StdRegions::StdExpansionSharedPtr &FaceExp,
365  const Array<OneD, const NekDouble> &inarray,
366  Array<OneD, NekDouble> &outarray,
367  StdRegions::Orientation orient)
368  {
369  int nq0 = m_base[0]->GetNumPoints();
370  int nq1 = m_base[1]->GetNumPoints();
371  int nq2 = m_base[2]->GetNumPoints();
372 
373  Array<OneD,NekDouble> o_tmp(nq0*nq1*nq2);
374 
375  if (orient == StdRegions::eNoOrientation)
376  {
377  orient = GetFaceOrient(face);
378  }
379 
380  switch(face)
381  {
382  case 0:
383  if(orient == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
384  {
385  //Directions A and B positive
386  Vmath::Vcopy(nq0*nq1,&(inarray[0]),1,&(o_tmp[0]),1);
387  }
388  else if(orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
389  {
390  //Direction A negative and B positive
391  for (int j=0; j<nq1; j++)
392  {
393  Vmath::Vcopy(nq0,&(inarray[0])+(nq0-1)+j*nq0,-1,&(o_tmp[0])+(j*nq0),1);
394  }
395  }
396  else if(orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2)
397  {
398  //Direction A positive and B negative
399  for (int j=0; j<nq1; j++)
400  {
401  Vmath::Vcopy(nq0,&(inarray[0])+nq0*(nq1-1-j),1,&(o_tmp[0])+(j*nq0),1);
402  }
403  }
404  else if(orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2)
405  {
406  //Direction A negative and B negative
407  for(int j=0; j<nq1; j++)
408  {
409  Vmath::Vcopy(nq0,&(inarray[0])+(nq0*nq1-1-j*nq0),-1,&(o_tmp[0])+(j*nq0),1);
410  }
411  }
412  else if(orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
413  {
414  //Transposed, Direction A and B positive
415  for (int i=0; i<nq0; i++)
416  {
417  Vmath::Vcopy(nq1,&(inarray[0])+i,nq0,&(o_tmp[0])+(i*nq1),1);
418  }
419  }
420  else if(orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1)
421  {
422  //Transposed, Direction A positive and B negative
423  for (int i=0; i<nq0; i++)
424  {
425  Vmath::Vcopy(nq1,&(inarray[0])+(nq0-1-i),nq0,&(o_tmp[0])+(i*nq1),1);
426  }
427  }
428  else if(orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1)
429  {
430  //Transposed, Direction A negative and B positive
431  for (int i=0; i<nq0; i++)
432  {
433  Vmath::Vcopy(nq1,&(inarray[0])+i+nq0*(nq1-1),-nq0,&(o_tmp[0])+(i*nq1),1);
434  }
435  }
436  else if(orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
437  {
438  //Transposed, Direction A and B negative
439  for (int i=0; i<nq0; i++)
440  {
441  Vmath::Vcopy(nq1,&(inarray[0])+(nq0*nq1-1-i),-nq0,&(o_tmp[0])+(i*nq1),1);
442  }
443  }
444  LibUtilities::Interp2D(m_base[0]->GetPointsKey(), m_base[1]->GetPointsKey(), o_tmp,
445  FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),outarray);
446  break;
447 
448  case 1:
449  {
450  for (int k = 0; k < nq2; k++)
451  {
452  Vmath::Vcopy(nq0,inarray.get()+nq0*nq1*k,1,outarray.get()+k*nq0,1);
453  }
454  LibUtilities::Interp2D(m_base[0]->GetPointsKey(), m_base[2]->GetPointsKey(), outarray.get(),
455  FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),o_tmp.get());
456  break;
457  }
458 
459  case 2:
460  {
461  Vmath::Vcopy(nq1*nq2,inarray.get()+(nq0-1),nq0,outarray.get(),1);
462  LibUtilities::Interp2D(m_base[1]->GetPointsKey(), m_base[2]->GetPointsKey(), outarray.get(),
463  FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),o_tmp.get());
464  break;
465  }
466 
467  case 3:
468  {
469  for (int k = 0; k < nq2; k++)
470  {
471  Vmath::Vcopy(nq0,inarray.get()+nq0*(nq1-1)+nq0*nq1*k,1,outarray.get()+(k*nq0),1);
472  }
473  LibUtilities::Interp2D(m_base[0]->GetPointsKey(), m_base[2]->GetPointsKey(), outarray.get(),
474  FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),o_tmp.get());
475  }
476 
477  case 4:
478  {
479  Vmath::Vcopy(nq1*nq2,inarray.get(),nq0,outarray.get(),1);
480  LibUtilities::Interp2D(m_base[1]->GetPointsKey(), m_base[2]->GetPointsKey(), outarray.get(),
481  FaceExp->GetBasis(0)->GetPointsKey(),FaceExp->GetBasis(1)->GetPointsKey(),o_tmp.get());
482  break;
483  }
484 
485  default:
486  ASSERTL0(false,"face value (> 4) is out of range");
487  break;
488  }
489 
490  if (face > 0)
491  {
492  int fnq1 = FaceExp->GetNumPoints(0);
493  int fnq2 = FaceExp->GetNumPoints(1);
494 
495  if ((int)orient == 7)
496  {
497  for (int j = 0; j < fnq2; ++j)
498  {
499  Vmath::Vcopy(fnq1, o_tmp.get()+((j+1)*fnq1-1), -1, outarray.get()+j*fnq1, 1);
500  }
501  }
502  else
503  {
504  Vmath::Vcopy(fnq1*fnq2, o_tmp.get(), 1, outarray.get(), 1);
505  }
506  }
507  }
508 
509  void PyrExp::v_ComputeFaceNormal(const int face)
510  {
511  const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
512  GetGeom()->GetMetricInfo();
513  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
514  SpatialDomains::GeomType type = geomFactors->GetGtype();
515  const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys);
516  const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
517 
518  // Number of quadrature points in face expansion.
519  int nq = m_base[0]->GetNumPoints()*m_base[0]->GetNumPoints();
520  int vCoordDim = GetCoordim();
521  int i;
522 
523  m_faceNormals[face] = Array<OneD, Array<OneD, NekDouble> >(vCoordDim);
524  Array<OneD, Array<OneD, NekDouble> > &normal = m_faceNormals[face];
525  for (i = 0; i < vCoordDim; ++i)
526  {
527  normal[i] = Array<OneD, NekDouble>(nq);
528  }
529 
530  // Regular geometry case
531  if (type == SpatialDomains::eRegular ||
532  type == SpatialDomains::eMovingRegular)
533  {
534  NekDouble fac;
535  // Set up normals
536  switch(face)
537  {
538  case 0:
539  {
540  for(i = 0; i < vCoordDim; ++i)
541  {
542  Vmath::Fill(nq,-df[3*i+2][0],normal[i],1);
543  }
544  break;
545  }
546  case 1:
547  {
548  for(i = 0; i < vCoordDim; ++i)
549  {
550  Vmath::Fill(nq,-df[3*i+1][0],normal[i],1);
551  }
552  break;
553  }
554  case 2:
555  {
556  for(i = 0; i < vCoordDim; ++i)
557  {
558  Vmath::Fill(nq,df[3*i][0]+df[3*i+2][0],normal[i],1);
559  }
560  break;
561  }
562  case 3:
563  {
564  for(i = 0; i < vCoordDim; ++i)
565  {
566  Vmath::Fill(nq,df[3*i+1][0]+df[3*i+2][0],normal[i],1);
567  }
568  break;
569  }
570  case 4:
571  {
572  for(i = 0; i < vCoordDim; ++i)
573  {
574  Vmath::Fill(nq,-df[3*i][0],normal[i],1);
575  }
576  break;
577  }
578  default:
579  ASSERTL0(false,"face is out of range (face < 4)");
580  }
581 
582  // Normalise resulting vector.
583  fac = 0.0;
584  for(i = 0; i < vCoordDim; ++i)
585  {
586  fac += normal[i][0]*normal[i][0];
587  }
588  fac = 1.0/sqrt(fac);
589  for (i = 0; i < vCoordDim; ++i)
590  {
591  Vmath::Smul(nq,fac,normal[i],1,normal[i],1);
592  }
593  }
594  else
595  {
596  // Set up deformed normals.
597  int j, k;
598 
599  int nq0 = ptsKeys[0].GetNumPoints();
600  int nq1 = ptsKeys[1].GetNumPoints();
601  int nq2 = ptsKeys[2].GetNumPoints();
602  int nq01 = nq0*nq1;
603  int nqtot;
604 
605  // Determine number of quadrature points on the face.
606  if (face == 0)
607  {
608  nqtot = nq0*nq1;
609  }
610  else if (face == 1 || face == 3)
611  {
612  nqtot = nq0*nq2;
613  }
614  else
615  {
616  nqtot = nq1*nq2;
617  }
618 
619  LibUtilities::PointsKey points0;
620  LibUtilities::PointsKey points1;
621 
622  Array<OneD, NekDouble> work (nq, 0.0);
623  Array<OneD, NekDouble> normals(vCoordDim*nqtot,0.0);
624 
625  // Extract Jacobian along face and recover local derivatives
626  // (dx/dr) for polynomial interpolation by multiplying m_gmat by
627  // jacobian
628  switch(face)
629  {
630  case 0:
631  {
632  for(j = 0; j < nq01; ++j)
633  {
634  normals[j] = -df[2][j]*jac[j];
635  normals[nqtot+j] = -df[5][j]*jac[j];
636  normals[2*nqtot+j] = -df[8][j]*jac[j];
637  }
638 
639  points0 = ptsKeys[0];
640  points1 = ptsKeys[1];
641  break;
642  }
643 
644  case 1:
645  {
646  for (j = 0; j < nq0; ++j)
647  {
648  for(k = 0; k < nq2; ++k)
649  {
650  int tmp = j+nq01*k;
651  normals[j+k*nq0] =
652  -df[1][tmp]*jac[tmp];
653  normals[nqtot+j+k*nq0] =
654  -df[4][tmp]*jac[tmp];
655  normals[2*nqtot+j+k*nq0] =
656  -df[7][tmp]*jac[tmp];
657  }
658  }
659 
660  points0 = ptsKeys[0];
661  points1 = ptsKeys[2];
662  break;
663  }
664 
665  case 2:
666  {
667  for (j = 0; j < nq1; ++j)
668  {
669  for(k = 0; k < nq2; ++k)
670  {
671  int tmp = nq0-1+nq0*j+nq01*k;
672  normals[j+k*nq1] =
673  (df[0][tmp]+df[2][tmp])*jac[tmp];
674  normals[nqtot+j+k*nq1] =
675  (df[3][tmp]+df[5][tmp])*jac[tmp];
676  normals[2*nqtot+j+k*nq1] =
677  (df[6][tmp]+df[8][tmp])*jac[tmp];
678  }
679  }
680 
681  points0 = ptsKeys[1];
682  points1 = ptsKeys[2];
683  break;
684  }
685 
686  case 3:
687  {
688  for (j = 0; j < nq0; ++j)
689  {
690  for(k = 0; k < nq2; ++k)
691  {
692  int tmp = nq0*(nq1-1) + j + nq01*k;
693  normals[j+k*nq0] =
694  (df[1][tmp]+df[2][tmp])*jac[tmp];
695  normals[nqtot+j+k*nq0] =
696  (df[4][tmp]+df[5][tmp])*jac[tmp];
697  normals[2*nqtot+j+k*nq0] =
698  (df[7][tmp]+df[8][tmp])*jac[tmp];
699  }
700  }
701 
702  points0 = ptsKeys[0];
703  points1 = ptsKeys[2];
704  break;
705  }
706 
707  case 4:
708  {
709  for (j = 0; j < nq1; ++j)
710  {
711  for(k = 0; k < nq2; ++k)
712  {
713  int tmp = j*nq0+nq01*k;
714  normals[j+k*nq1] =
715  -df[0][tmp]*jac[tmp];
716  normals[nqtot+j+k*nq1] =
717  -df[3][tmp]*jac[tmp];
718  normals[2*nqtot+j+k*nq1] =
719  -df[6][tmp]*jac[tmp];
720  }
721  }
722 
723  points0 = ptsKeys[1];
724  points1 = ptsKeys[2];
725  break;
726  }
727 
728  default:
729  ASSERTL0(false,"face is out of range (face < 4)");
730  }
731 
732  // Interpolate Jacobian and invert
733  LibUtilities::Interp2D(points0, points1, jac,
734  m_base[0]->GetPointsKey(),
735  m_base[0]->GetPointsKey(),
736  work);
737  Vmath::Sdiv(nq, 1.0, &work[0], 1, &work[0], 1);
738 
739  // Interpolate normal and multiply by inverse Jacobian.
740  for(i = 0; i < vCoordDim; ++i)
741  {
742  LibUtilities::Interp2D(points0, points1,
743  &normals[i*nqtot],
744  m_base[0]->GetPointsKey(),
745  m_base[0]->GetPointsKey(),
746  &normal[i][0]);
747  Vmath::Vmul(nq,work,1,normal[i],1,normal[i],1);
748  }
749 
750  // Normalise to obtain unit normals.
751  Vmath::Zero(nq,work,1);
752  for(i = 0; i < GetCoordim(); ++i)
753  {
754  Vmath::Vvtvp(nq,normal[i],1,normal[i],1,work,1,work,1);
755  }
756 
757  Vmath::Vsqrt(nq,work,1,work,1);
758  Vmath::Sdiv (nq,1.0,work,1,work,1);
759 
760  for(i = 0; i < GetCoordim(); ++i)
761  {
762  Vmath::Vmul(nq,normal[i],1,work,1,normal[i],1);
763  }
764  }
765  }
766 
767  //---------------------------------------
768  // Matrix creation functions
769  //---------------------------------------
770 
771  DNekMatSharedPtr PyrExp::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
772  {
773  DNekMatSharedPtr returnval;
774 
775  switch(mkey.GetMatrixType())
776  {
777  case StdRegions::eHybridDGHelmholtz:
778  case StdRegions::eHybridDGLamToU:
779  case StdRegions::eHybridDGLamToQ0:
780  case StdRegions::eHybridDGLamToQ1:
781  case StdRegions::eHybridDGLamToQ2:
782  case StdRegions::eHybridDGHelmBndLam:
783  returnval = Expansion3D::v_GenMatrix(mkey);
784  break;
785  default:
786  returnval = StdPyrExp::v_GenMatrix(mkey);
787  }
788 
789  return returnval;
790  }
791 
792  DNekMatSharedPtr PyrExp::v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
793  {
794  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
795  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
796  LibUtilities::BasisKey bkey2 = m_base[2]->GetBasisKey();
797  StdRegions::StdPyrExpSharedPtr tmp =
798  MemoryManager<StdPyrExp>::AllocateSharedPtr(bkey0, bkey1, bkey2);
799 
800  return tmp->GetStdMatrix(mkey);
801  }
802 
803  DNekScalMatSharedPtr PyrExp::v_GetLocMatrix(const MatrixKey &mkey)
804  {
805  return m_matrixManager[mkey];
806  }
807 
808  DNekScalBlkMatSharedPtr PyrExp::v_GetLocStaticCondMatrix(const MatrixKey &mkey)
809  {
810  return m_staticCondMatrixManager[mkey];
811  }
812 
813  void PyrExp::v_DropLocStaticCondMatrix(const MatrixKey &mkey)
814  {
815  m_staticCondMatrixManager.DeleteObject(mkey);
816  }
817 
818  DNekScalMatSharedPtr PyrExp::CreateMatrix(const MatrixKey &mkey)
819  {
820  DNekScalMatSharedPtr returnval;
821  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
822 
823  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
824 
825  switch(mkey.GetMatrixType())
826  {
827  case StdRegions::eMass:
828  {
829  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
830  {
831  NekDouble one = 1.0;
832  DNekMatSharedPtr mat = GenMatrix(mkey);
833  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
834  }
835  else
836  {
837  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
838  DNekMatSharedPtr mat = GetStdMatrix(mkey);
839  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
840  }
841  }
842  break;
843  case StdRegions::eInvMass:
844  {
845  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
846  {
847  NekDouble one = 1.0;
848  StdRegions::StdMatrixKey masskey(StdRegions::eMass,DetShapeType(),
849  *this);
850  DNekMatSharedPtr mat = GenMatrix(masskey);
851  mat->Invert();
852  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
853  }
854  else
855  {
856  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
857  DNekMatSharedPtr mat = GetStdMatrix(mkey);
858  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(fac,mat);
859  }
860  }
861  break;
862  case StdRegions::eLaplacian:
863  {
864  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
865  mkey.GetNVarCoeff() > 0 ||
866  mkey.ConstFactorExists(
867  StdRegions::eFactorSVVCutoffRatio))
868  {
869  NekDouble one = 1.0;
870  DNekMatSharedPtr mat = GenMatrix(mkey);
871 
872  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
873  }
874  else
875  {
876  MatrixKey lap00key(StdRegions::eLaplacian00,
877  mkey.GetShapeType(), *this);
878  MatrixKey lap01key(StdRegions::eLaplacian01,
879  mkey.GetShapeType(), *this);
880  MatrixKey lap02key(StdRegions::eLaplacian02,
881  mkey.GetShapeType(), *this);
882  MatrixKey lap11key(StdRegions::eLaplacian11,
883  mkey.GetShapeType(), *this);
884  MatrixKey lap12key(StdRegions::eLaplacian12,
885  mkey.GetShapeType(), *this);
886  MatrixKey lap22key(StdRegions::eLaplacian22,
887  mkey.GetShapeType(), *this);
888 
889  DNekMat &lap00 = *GetStdMatrix(lap00key);
890  DNekMat &lap01 = *GetStdMatrix(lap01key);
891  DNekMat &lap02 = *GetStdMatrix(lap02key);
892  DNekMat &lap11 = *GetStdMatrix(lap11key);
893  DNekMat &lap12 = *GetStdMatrix(lap12key);
894  DNekMat &lap22 = *GetStdMatrix(lap22key);
895 
896  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
897  Array<TwoD, const NekDouble> gmat =
898  m_metricinfo->GetGmat(ptsKeys);
899 
900  int rows = lap00.GetRows();
901  int cols = lap00.GetColumns();
902 
903  DNekMatSharedPtr lap = MemoryManager<DNekMat>
904  ::AllocateSharedPtr(rows,cols);
905 
906  (*lap) = gmat[0][0]*lap00
907  + gmat[4][0]*lap11
908  + gmat[8][0]*lap22
909  + gmat[3][0]*(lap01 + Transpose(lap01))
910  + gmat[6][0]*(lap02 + Transpose(lap02))
911  + gmat[7][0]*(lap12 + Transpose(lap12));
912 
913  returnval = MemoryManager<DNekScalMat>
914  ::AllocateSharedPtr(jac, lap);
915  }
916  }
917  break;
918  case StdRegions::eHelmholtz:
919  {
920  NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorLambda);
921  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
922  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
923  MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
924  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
925 
926  int rows = LapMat.GetRows();
927  int cols = LapMat.GetColumns();
928 
929  DNekMatSharedPtr helm = MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
930 
931  (*helm) = LapMat + factor*MassMat;
932 
933  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, helm);
934  }
935  break;
936  default:
937  NEKERROR(ErrorUtil::efatal, "Matrix creation not defined");
938  break;
939  }
940 
941  return returnval;
942  }
943 
944  DNekScalBlkMatSharedPtr PyrExp::CreateStaticCondMatrix(const MatrixKey &mkey)
945  {
946  DNekScalBlkMatSharedPtr returnval;
947 
948  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
949 
950  // set up block matrix system
951  unsigned int nbdry = NumBndryCoeffs();
952  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
953  unsigned int exp_size[] = {nbdry, nint};
954  unsigned int nblks = 2;
955  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks, nblks, exp_size, exp_size); //Really need a constructor which takes Arrays
956  NekDouble factor = 1.0;
957 
958  switch(mkey.GetMatrixType())
959  {
960  case StdRegions::eLaplacian:
961  case StdRegions::eHelmholtz: // special case since Helmholtz not defined in StdRegions
962 
963  // use Deformed case for both regular and deformed geometries
964  factor = 1.0;
965  goto UseLocRegionsMatrix;
966  break;
967  default:
968  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
969  {
970  factor = 1.0;
971  goto UseLocRegionsMatrix;
972  }
973  else
974  {
975  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
976  factor = mat->Scale();
977  goto UseStdRegionsMatrix;
978  }
979  break;
980  UseStdRegionsMatrix:
981  {
982  NekDouble invfactor = 1.0/factor;
983  NekDouble one = 1.0;
984  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
985  DNekScalMatSharedPtr Atmp;
986  DNekMatSharedPtr Asubmat;
987 
988  //TODO: check below
989  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
990  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
991  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
992  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
993  }
994  break;
995  UseLocRegionsMatrix:
996  {
997  int i,j;
998  NekDouble invfactor = 1.0/factor;
999  NekDouble one = 1.0;
1000  DNekScalMat &mat = *GetLocMatrix(mkey);
1001  DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
1002  DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
1003  DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
1004  DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
1005 
1006  Array<OneD,unsigned int> bmap(nbdry);
1007  Array<OneD,unsigned int> imap(nint);
1008  GetBoundaryMap(bmap);
1009  GetInteriorMap(imap);
1010 
1011  for(i = 0; i < nbdry; ++i)
1012  {
1013  for(j = 0; j < nbdry; ++j)
1014  {
1015  (*A)(i,j) = mat(bmap[i],bmap[j]);
1016  }
1017 
1018  for(j = 0; j < nint; ++j)
1019  {
1020  (*B)(i,j) = mat(bmap[i],imap[j]);
1021  }
1022  }
1023 
1024  for(i = 0; i < nint; ++i)
1025  {
1026  for(j = 0; j < nbdry; ++j)
1027  {
1028  (*C)(i,j) = mat(imap[i],bmap[j]);
1029  }
1030 
1031  for(j = 0; j < nint; ++j)
1032  {
1033  (*D)(i,j) = mat(imap[i],imap[j]);
1034  }
1035  }
1036 
1037  // Calculate static condensed system
1038  if(nint)
1039  {
1040  D->Invert();
1041  (*B) = (*B)*(*D);
1042  (*A) = (*A) - (*B)*(*C);
1043  }
1044 
1045  DNekScalMatSharedPtr Atmp;
1046 
1047  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
1048  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
1049  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
1050  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
1051 
1052  }
1053  }
1054  return returnval;
1055  }
1056 
1057  void PyrExp::v_ComputeLaplacianMetric()
1058  {
1059  if (m_metrics.count(MetricQuadrature) == 0)
1060  {
1061  ComputeQuadratureMetric();
1062  }
1063 
1064  int i, j;
1065  const unsigned int nqtot = GetTotPoints();
1066  const unsigned int dim = 3;
1067  const MetricType m[3][3] = {
1068  { MetricLaplacian00, MetricLaplacian01, MetricLaplacian02 },
1069  { MetricLaplacian01, MetricLaplacian11, MetricLaplacian12 },
1070  { MetricLaplacian02, MetricLaplacian12, MetricLaplacian22 }
1071  };
1072 
1073  for (unsigned int i = 0; i < dim; ++i)
1074  {
1075  for (unsigned int j = i; j < dim; ++j)
1076  {
1077  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1078  }
1079  }
1080 
1081  // Define shorthand synonyms for m_metrics storage
1082  Array<OneD,NekDouble> g0 (m_metrics[m[0][0]]);
1083  Array<OneD,NekDouble> g1 (m_metrics[m[1][1]]);
1084  Array<OneD,NekDouble> g2 (m_metrics[m[2][2]]);
1085  Array<OneD,NekDouble> g3 (m_metrics[m[0][1]]);
1086  Array<OneD,NekDouble> g4 (m_metrics[m[0][2]]);
1087  Array<OneD,NekDouble> g5 (m_metrics[m[1][2]]);
1088 
1089  // Allocate temporary storage
1090  Array<OneD,NekDouble> alloc(9*nqtot,0.0);
1091  Array<OneD,NekDouble> h0 (nqtot, alloc);
1092  Array<OneD,NekDouble> h1 (nqtot, alloc+ 1*nqtot);
1093  Array<OneD,NekDouble> h2 (nqtot, alloc+ 2*nqtot);
1094  Array<OneD,NekDouble> wsp1 (nqtot, alloc+ 3*nqtot);
1095  Array<OneD,NekDouble> wsp2 (nqtot, alloc+ 4*nqtot);
1096  Array<OneD,NekDouble> wsp3 (nqtot, alloc+ 5*nqtot);
1097  Array<OneD,NekDouble> wsp4 (nqtot, alloc+ 6*nqtot);
1098  Array<OneD,NekDouble> wsp5 (nqtot, alloc+ 7*nqtot);
1099  Array<OneD,NekDouble> wsp6 (nqtot, alloc+ 8*nqtot);
1100 
1101  const Array<TwoD, const NekDouble>& df =
1102  m_metricinfo->GetDerivFactors(GetPointsKeys());
1103  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
1104  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
1105  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1106  const unsigned int nquad0 = m_base[0]->GetNumPoints();
1107  const unsigned int nquad1 = m_base[1]->GetNumPoints();
1108  const unsigned int nquad2 = m_base[2]->GetNumPoints();
1109 
1110  // Populate collapsed coordinate arrays h0, h1 and h2.
1111  for(j = 0; j < nquad2; ++j)
1112  {
1113  for(i = 0; i < nquad1; ++i)
1114  {
1115  Vmath::Fill(nquad0, 2.0/(1.0-z2[j]), &h0[0]+i*nquad0 + j*nquad0*nquad1,1);
1116  Vmath::Fill(nquad0, 1.0/(1.0-z2[j]), &h1[0]+i*nquad0 + j*nquad0*nquad1,1);
1117  Vmath::Fill(nquad0, (1.0+z1[i])/(1.0-z2[j]), &h2[0]+i*nquad0 + j*nquad0*nquad1,1);
1118  }
1119  }
1120  for(i = 0; i < nquad0; i++)
1121  {
1122  Blas::Dscal(nquad1*nquad2, 1+z0[i], &h1[0]+i, nquad0);
1123  }
1124 
1125  // Step 3. Construct combined metric terms for physical space to
1126  // collapsed coordinate system.
1127  // Order of construction optimised to minimise temporary storage
1128  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1129  {
1130  // f_{1k}
1131  Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &df[2][0], 1, &h1[0], 1, &wsp1[0], 1);
1132  Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &df[5][0], 1, &h1[0], 1, &wsp2[0], 1);
1133  Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &df[8][0], 1, &h1[0], 1, &wsp3[0], 1);
1134 
1135  // g0
1136  Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp1[0], 1, &wsp2[0], 1, &wsp2[0], 1, &g0[0], 1);
1137  Vmath::Vvtvp (nqtot, &wsp3[0], 1, &wsp3[0], 1, &g0[0], 1, &g0[0], 1);
1138 
1139  // g4
1140  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp1[0], 1, &df[5][0], 1, &wsp2[0], 1, &g4[0], 1);
1141  Vmath::Vvtvp (nqtot, &df[8][0], 1, &wsp3[0], 1, &g4[0], 1, &g4[0], 1);
1142 
1143  // f_{2k}
1144  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &h0[0], 1, &df[2][0], 1, &h2[0], 1, &wsp4[0], 1);
1145  Vmath::Vvtvvtp(nqtot, &df[4][0], 1, &h0[0], 1, &df[5][0], 1, &h2[0], 1, &wsp5[0], 1);
1146  Vmath::Vvtvvtp(nqtot, &df[7][0], 1, &h0[0], 1, &df[8][0], 1, &h2[0], 1, &wsp6[0], 1);
1147 
1148  // g1
1149  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0], 1, &g1[0], 1);
1150  Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp6[0], 1, &g1[0], 1, &g1[0], 1);
1151 
1152  // g3
1153  Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp4[0], 1, &wsp2[0], 1, &wsp5[0], 1, &g3[0], 1);
1154  Vmath::Vvtvp (nqtot, &wsp3[0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1155 
1156  // g5
1157  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0], 1, &g5[0], 1);
1158  Vmath::Vvtvp (nqtot, &df[8][0], 1, &wsp6[0], 1, &g5[0], 1, &g5[0], 1);
1159 
1160  // g2
1161  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &df[2][0], 1, &df[5][0], 1, &df[5][0], 1, &g2[0], 1);
1162  Vmath::Vvtvp (nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1163  }
1164  else
1165  {
1166  // f_{1k}
1167  Vmath::Svtsvtp(nqtot, df[0][0], &h0[0], 1, df[2][0], &h1[0], 1, &wsp1[0], 1);
1168  Vmath::Svtsvtp(nqtot, df[3][0], &h0[0], 1, df[5][0], &h1[0], 1, &wsp2[0], 1);
1169  Vmath::Svtsvtp(nqtot, df[6][0], &h0[0], 1, df[8][0], &h1[0], 1, &wsp3[0], 1);
1170 
1171  // g0
1172  Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp1[0], 1, &wsp2[0], 1, &wsp2[0], 1, &g0[0], 1);
1173  Vmath::Vvtvp (nqtot, &wsp3[0], 1, &wsp3[0], 1, &g0[0], 1, &g0[0], 1);
1174 
1175  // g4
1176  Vmath::Svtsvtp(nqtot, df[2][0], &wsp1[0], 1, df[5][0], &wsp2[0], 1, &g4[0], 1);
1177  Vmath::Svtvp (nqtot, df[8][0], &wsp3[0], 1, &g4[0], 1, &g4[0], 1);
1178 
1179  // f_{2k}
1180  Vmath::Svtsvtp(nqtot, df[1][0], &h0[0], 1, df[2][0], &h2[0], 1, &wsp4[0], 1);
1181  Vmath::Svtsvtp(nqtot, df[4][0], &h0[0], 1, df[5][0], &h2[0], 1, &wsp5[0], 1);
1182  Vmath::Svtsvtp(nqtot, df[7][0], &h0[0], 1, df[8][0], &h2[0], 1, &wsp6[0], 1);
1183 
1184  // g1
1185  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0], 1, &g1[0], 1);
1186  Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp6[0], 1, &g1[0], 1, &g1[0], 1);
1187 
1188  // g3
1189  Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp4[0], 1, &wsp2[0], 1, &wsp5[0], 1, &g3[0], 1);
1190  Vmath::Vvtvp (nqtot, &wsp3[0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1191 
1192  // g5
1193  Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1, &g5[0], 1);
1194  Vmath::Svtvp (nqtot, df[8][0], &wsp6[0], 1, &g5[0], 1, &g5[0], 1);
1195 
1196  // g2
1197  Vmath::Fill(nqtot, df[2][0]*df[2][0] + df[5][0]*df[5][0] + df[8][0]*df[8][0], &g2[0], 1);
1198  }
1199 
1200  for (unsigned int i = 0; i < dim; ++i)
1201  {
1202  for (unsigned int j = i; j < dim; ++j)
1203  {
1204  MultiplyByQuadratureMetric(m_metrics[m[i][j]],
1205  m_metrics[m[i][j]]);
1206 
1207  }
1208  }
1209  }
1210 
1211  void PyrExp::v_LaplacianMatrixOp_MatFree_Kernel(
1212  const Array<OneD, const NekDouble> &inarray,
1213  Array<OneD, NekDouble> &outarray,
1214  Array<OneD, NekDouble> &wsp)
1215  {
1216  // This implementation is only valid when there are no coefficients
1217  // associated to the Laplacian operator
1218  if (m_metrics.count(MetricLaplacian00) == 0)
1219  {
1220  ComputeLaplacianMetric();
1221  }
1222 
1223  int nquad0 = m_base[0]->GetNumPoints();
1224  int nquad1 = m_base[1]->GetNumPoints();
1225  int nq2 = m_base[2]->GetNumPoints();
1226  int nqtot = nquad0*nquad1*nq2;
1227 
1228  ASSERTL1(wsp.num_elements() >= 6*nqtot,
1229  "Insufficient workspace size.");
1230  ASSERTL1(m_ncoeffs <= nqtot,
1231  "Workspace not set up for ncoeffs > nqtot");
1232 
1233  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
1234  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
1235  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata();
1236  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
1237  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
1238  const Array<OneD, const NekDouble>& dbase2 = m_base[2]->GetDbdata();
1239  const Array<OneD, const NekDouble>& metric00 = m_metrics[MetricLaplacian00];
1240  const Array<OneD, const NekDouble>& metric01 = m_metrics[MetricLaplacian01];
1241  const Array<OneD, const NekDouble>& metric02 = m_metrics[MetricLaplacian02];
1242  const Array<OneD, const NekDouble>& metric11 = m_metrics[MetricLaplacian11];
1243  const Array<OneD, const NekDouble>& metric12 = m_metrics[MetricLaplacian12];
1244  const Array<OneD, const NekDouble>& metric22 = m_metrics[MetricLaplacian22];
1245 
1246  // Allocate temporary storage
1247  Array<OneD,NekDouble> wsp0 (2*nqtot, wsp);
1248  Array<OneD,NekDouble> wsp1 ( nqtot, wsp+1*nqtot);
1249  Array<OneD,NekDouble> wsp2 ( nqtot, wsp+2*nqtot);
1250  Array<OneD,NekDouble> wsp3 ( nqtot, wsp+3*nqtot);
1251  Array<OneD,NekDouble> wsp4 ( nqtot, wsp+4*nqtot);
1252  Array<OneD,NekDouble> wsp5 ( nqtot, wsp+5*nqtot);
1253 
1254  // LAPLACIAN MATRIX OPERATION
1255  // wsp1 = du_dxi1 = D_xi1 * inarray = D_xi1 * u
1256  // wsp2 = du_dxi2 = D_xi2 * inarray = D_xi2 * u
1257  // wsp2 = du_dxi3 = D_xi3 * inarray = D_xi3 * u
1258  StdExpansion3D::PhysTensorDeriv(inarray,wsp0,wsp1,wsp2);
1259 
1260  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1261  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1262  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1263  // especially for this purpose
1264  Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp0[0],1,&metric01[0],1,&wsp1[0],1,&wsp3[0],1);
1265  Vmath::Vvtvp (nqtot,&metric02[0],1,&wsp2[0],1,&wsp3[0],1,&wsp3[0],1);
1266  Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp0[0],1,&metric11[0],1,&wsp1[0],1,&wsp4[0],1);
1267  Vmath::Vvtvp (nqtot,&metric12[0],1,&wsp2[0],1,&wsp4[0],1,&wsp4[0],1);
1268  Vmath::Vvtvvtp(nqtot,&metric02[0],1,&wsp0[0],1,&metric12[0],1,&wsp1[0],1,&wsp5[0],1);
1269  Vmath::Vvtvp (nqtot,&metric22[0],1,&wsp2[0],1,&wsp5[0],1,&wsp5[0],1);
1270 
1271  // outarray = m = (D_xi1 * B)^T * k
1272  // wsp1 = n = (D_xi2 * B)^T * l
1273  IProductWRTBase_SumFacKernel(dbase0,base1,base2,wsp3,outarray,wsp0,false,true,true);
1274  IProductWRTBase_SumFacKernel(base0,dbase1,base2,wsp4,wsp2, wsp0,true,false,true);
1275  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
1276  IProductWRTBase_SumFacKernel(base0,base1,dbase2,wsp5,wsp2, wsp0,true,true,false);
1277  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
1278  }
1279  }//end of namespace
1280 }//end of namespace
Generated on Sun Mar 15 2015 00:11:34 for Nektar++ by   doxygen 1.8.1.2