Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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>
38 
39 namespace Nektar
40 {
41  namespace LocalRegions
42  {
43 
45  const LibUtilities::BasisKey &Bb,
46  const LibUtilities::BasisKey &Bc,
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 
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 
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  */
111  {
112  int nquad0 = m_base[0]->GetNumPoints();
113  int nquad1 = m_base[1]->GetNumPoints();
114  int nquad2 = m_base[2]->GetNumPoints();
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 
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();
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  */
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
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  */
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();
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 
305  {
307  ::AllocateSharedPtr(m_base[0]->GetBasisKey(),
308  m_base[1]->GetBasisKey(),
309  m_base[2]->GetBasisKey());
310  }
311 
313  {
315  2, m_base[0]->GetPointsKey());
317  2, m_base[1]->GetPointsKey());
319  2, m_base[2]->GetPointsKey());
320 
322  ::AllocateSharedPtr( bkey0, bkey1, bkey2);
323  }
324 
325  /*
326  * @brief Get the coordinates #coords at the local coordinates
327  * #Lcoords
328  */
330  Array<OneD, NekDouble>& coords)
331  {
332  int i;
333 
334  ASSERTL1(Lcoords[0] <= -1.0 && Lcoords[0] >= 1.0 &&
335  Lcoords[1] <= -1.0 && Lcoords[1] >= 1.0 &&
336  Lcoords[2] <= -1.0 && Lcoords[2] >= 1.0,
337  "Local coordinates are not in region [-1,1]");
338 
339  // m_geom->FillGeom(); // TODO: implement FillGeom()
340 
341  for(i = 0; i < m_geom->GetCoordim(); ++i)
342  {
343  coords[i] = m_geom->GetCoord(i,Lcoords);
344  }
345  }
346 
348  Array<OneD, NekDouble> &coords_1,
349  Array<OneD, NekDouble> &coords_2,
350  Array<OneD, NekDouble> &coords_3)
351  {
352  Expansion::v_GetCoords(coords_1, coords_2, coords_3);
353  }
354 
356  const Array<OneD, const NekDouble>& physvals)
357  {
358  Array<OneD,NekDouble> Lcoord(3);
359 
360  ASSERTL0(m_geom,"m_geom not defined");
361 
362  //TODO: check GetLocCoords()
363  m_geom->GetLocCoords(coord, Lcoord);
364 
365  return StdPyrExp::v_PhysEvaluate(Lcoord, physvals);
366  }
367 
368 
369  //---------------------------------------
370  // Helper functions
371  //---------------------------------------
372 
374  {
375  return m_geom->GetCoordim();
376  }
377 
378  void PyrExp::v_GetFacePhysMap(const int face,
379  Array<OneD, int> &outarray)
380  {
381  int nquad0 = m_base[0]->GetNumPoints();
382  int nquad1 = m_base[1]->GetNumPoints();
383  int nquad2 = m_base[2]->GetNumPoints();
384 
385  int nq0 = 0;
386  int nq1 = 0;
387 
388  switch(face)
389  {
390  case 0:
391  nq0 = nquad0;
392  nq1 = nquad1;
393  if(outarray.num_elements()!=nq0*nq1)
394  {
395  outarray = Array<OneD, int>(nq0*nq1);
396  }
397 
398  //Directions A and B positive
399  for(int i = 0; i < nquad0*nquad1; ++i)
400  {
401  outarray[i] = i;
402  }
403 
404  break;
405  case 1:
406  nq0 = nquad0;
407  nq1 = nquad2;
408  if(outarray.num_elements()!=nq0*nq1)
409  {
410  outarray = Array<OneD, int>(nq0*nq1);
411  }
412 
413  //Direction A and B positive
414  for (int k=0; k<nquad2; k++)
415  {
416  for(int i = 0; i < nquad0; ++i)
417  {
418  outarray[k*nquad0+i] = (nquad0*nquad1*k)+i;
419  }
420  }
421 
422  break;
423  case 2:
424  nq0 = nquad1;
425  nq1 = nquad2;
426  if(outarray.num_elements()!=nq0*nq1)
427  {
428  outarray = Array<OneD, int>(nq0*nq1);
429  }
430 
431  //Directions A and B positive
432  for(int j = 0; j < nquad1*nquad2; ++j)
433  {
434  outarray[j] = nquad0-1 + j*nquad0;
435 
436  }
437  break;
438  case 3:
439 
440  nq0 = nquad0;
441  nq1 = nquad2;
442  if(outarray.num_elements()!=nq0*nq1)
443  {
444  outarray = Array<OneD, int>(nq0*nq1);
445  }
446 
447  //Direction A and B positive
448  for (int k=0; k<nquad2; k++)
449  {
450  for(int i = 0; i < nquad0; ++i)
451  {
452  outarray[k*nquad0+i] = nquad0*(nquad1-1) + (nquad0*nquad1*k)+i;
453  }
454  }
455 
456  case 4:
457  nq0 = nquad1;
458  nq1 = nquad2;
459 
460  if(outarray.num_elements()!=nq0*nq1)
461  {
462  outarray = Array<OneD, int>(nq0*nq1);
463  }
464 
465  //Directions A and B positive
466  for(int j = 0; j < nquad1*nquad2; ++j)
467  {
468  outarray[j] = j*nquad0;
469 
470  }
471  break;
472  default:
473  ASSERTL0(false,"face value (> 4) is out of range");
474  break;
475  }
476  }
477 
478  void PyrExp::v_ComputeFaceNormal(const int face)
479  {
480  const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
481  GetGeom()->GetMetricInfo();
483  SpatialDomains::GeomType type = geomFactors->GetGtype();
484  const Array<TwoD, const NekDouble> &df = geomFactors->GetDerivFactors(ptsKeys);
485  const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
486 
487  LibUtilities::BasisKey tobasis0 = DetFaceBasisKey(face,0);
488  LibUtilities::BasisKey tobasis1 = DetFaceBasisKey(face,1);
489 
490  // Number of quadrature points in face expansion.
491  int nq_face = tobasis0.GetNumPoints()*tobasis1.GetNumPoints();
492 
493  int vCoordDim = GetCoordim();
494  int i;
495 
498  for (i = 0; i < vCoordDim; ++i)
499  {
500  normal[i] = Array<OneD, NekDouble>(nq_face);
501  }
502 
503  // Regular geometry case
504  if (type == SpatialDomains::eRegular ||
506  {
507  NekDouble fac;
508  // Set up normals
509  switch(face)
510  {
511  case 0:
512  {
513  for(i = 0; i < vCoordDim; ++i)
514  {
515  normal[i][0] = -df[3*i+2][0];
516  }
517  break;
518  }
519  case 1:
520  {
521  for(i = 0; i < vCoordDim; ++i)
522  {
523  normal[i][0] = -df[3*i+1][0];
524  }
525  break;
526  }
527  case 2:
528  {
529  for(i = 0; i < vCoordDim; ++i)
530  {
531  normal[i][0] = df[3*i][0]+df[3*i+2][0];
532  }
533  break;
534  }
535  case 3:
536  {
537  for(i = 0; i < vCoordDim; ++i)
538  {
539  normal[i][0] = df[3*i+1][0]+df[3*i+2][0];
540  }
541  break;
542  }
543  case 4:
544  {
545  for(i = 0; i < vCoordDim; ++i)
546  {
547  normal[i][0] = -df[3*i][0];
548  }
549  break;
550  }
551  default:
552  ASSERTL0(false,"face is out of range (face < 4)");
553  }
554 
555  // Normalise resulting vector.
556  fac = 0.0;
557  for(i = 0; i < vCoordDim; ++i)
558  {
559  fac += normal[i][0]*normal[i][0];
560  }
561  fac = 1.0/sqrt(fac);
562  for (i = 0; i < vCoordDim; ++i)
563  {
564  Vmath::Fill(nq_face,fac*normal[i][0],normal[i],1);
565  }
566  }
567  else
568  {
569  // Set up deformed normals.
570  int j, k;
571 
572  int nq0 = ptsKeys[0].GetNumPoints();
573  int nq1 = ptsKeys[1].GetNumPoints();
574  int nq2 = ptsKeys[2].GetNumPoints();
575  int nq01 = nq0*nq1;
576  int nqtot;
577 
578  // Determine number of quadrature points on the face.
579  if (face == 0)
580  {
581  nqtot = nq0*nq1;
582  }
583  else if (face == 1 || face == 3)
584  {
585  nqtot = nq0*nq2;
586  }
587  else
588  {
589  nqtot = nq1*nq2;
590  }
591 
592  LibUtilities::PointsKey points0;
593  LibUtilities::PointsKey points1;
594 
595  Array<OneD, NekDouble> faceJac(nqtot);
596  Array<OneD, NekDouble> normals(vCoordDim*nqtot,0.0);
597 
598  // Extract Jacobian along face and recover local derivatives
599  // (dx/dr) for polynomial interpolation by multiplying m_gmat by
600  // jacobian
601  switch(face)
602  {
603  case 0:
604  {
605  for(j = 0; j < nq01; ++j)
606  {
607  normals[j] = -df[2][j]*jac[j];
608  normals[nqtot+j] = -df[5][j]*jac[j];
609  normals[2*nqtot+j] = -df[8][j]*jac[j];
610  faceJac[j] = jac[j];
611  }
612 
613  points0 = ptsKeys[0];
614  points1 = ptsKeys[1];
615  break;
616  }
617 
618  case 1:
619  {
620  for (j = 0; j < nq0; ++j)
621  {
622  for(k = 0; k < nq2; ++k)
623  {
624  int tmp = j+nq01*k;
625  normals[j+k*nq0] =
626  -df[1][tmp]*jac[tmp];
627  normals[nqtot+j+k*nq0] =
628  -df[4][tmp]*jac[tmp];
629  normals[2*nqtot+j+k*nq0] =
630  -df[7][tmp]*jac[tmp];
631  faceJac[j+k*nq0] = jac[tmp];
632  }
633  }
634 
635  points0 = ptsKeys[0];
636  points1 = ptsKeys[2];
637  break;
638  }
639 
640  case 2:
641  {
642  for (j = 0; j < nq1; ++j)
643  {
644  for(k = 0; k < nq2; ++k)
645  {
646  int tmp = nq0-1+nq0*j+nq01*k;
647  normals[j+k*nq1] =
648  (df[0][tmp]+df[2][tmp])*jac[tmp];
649  normals[nqtot+j+k*nq1] =
650  (df[3][tmp]+df[5][tmp])*jac[tmp];
651  normals[2*nqtot+j+k*nq1] =
652  (df[6][tmp]+df[8][tmp])*jac[tmp];
653  faceJac[j+k*nq1] = jac[tmp];
654  }
655  }
656 
657  points0 = ptsKeys[1];
658  points1 = ptsKeys[2];
659  break;
660  }
661 
662  case 3:
663  {
664  for (j = 0; j < nq0; ++j)
665  {
666  for(k = 0; k < nq2; ++k)
667  {
668  int tmp = nq0*(nq1-1) + j + nq01*k;
669  normals[j+k*nq0] =
670  (df[1][tmp]+df[2][tmp])*jac[tmp];
671  normals[nqtot+j+k*nq0] =
672  (df[4][tmp]+df[5][tmp])*jac[tmp];
673  normals[2*nqtot+j+k*nq0] =
674  (df[7][tmp]+df[8][tmp])*jac[tmp];
675  faceJac[j+k*nq0] = jac[tmp];
676  }
677  }
678 
679  points0 = ptsKeys[0];
680  points1 = ptsKeys[2];
681  break;
682  }
683 
684  case 4:
685  {
686  for (j = 0; j < nq1; ++j)
687  {
688  for(k = 0; k < nq2; ++k)
689  {
690  int tmp = j*nq0+nq01*k;
691  normals[j+k*nq1] =
692  -df[0][tmp]*jac[tmp];
693  normals[nqtot+j+k*nq1] =
694  -df[3][tmp]*jac[tmp];
695  normals[2*nqtot+j+k*nq1] =
696  -df[6][tmp]*jac[tmp];
697  faceJac[j+k*nq1] = jac[tmp];
698  }
699  }
700 
701  points0 = ptsKeys[1];
702  points1 = ptsKeys[2];
703  break;
704  }
705 
706  default:
707  ASSERTL0(false,"face is out of range (face < 4)");
708  }
709 
710  Array<OneD, NekDouble> work (nq_face, 0.0);
711  // Interpolate Jacobian and invert
712  LibUtilities::Interp2D(points0, points1, faceJac,
713  tobasis0.GetPointsKey(),
714  tobasis1.GetPointsKey(),
715  work);
716  Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
717 
718  // Interpolate normal and multiply by inverse Jacobian.
719  for(i = 0; i < vCoordDim; ++i)
720  {
721  LibUtilities::Interp2D(points0, points1,
722  &normals[i*nqtot],
723  tobasis0.GetPointsKey(),
724  tobasis1.GetPointsKey(),
725  &normal[i][0]);
726  Vmath::Vmul(nq_face,work,1,normal[i],1,normal[i],1);
727  }
728 
729  // Normalise to obtain unit normals.
730  Vmath::Zero(nq_face,work,1);
731  for(i = 0; i < GetCoordim(); ++i)
732  {
733  Vmath::Vvtvp(nq_face,normal[i],1,normal[i],1,work,1,work,1);
734  }
735 
736  Vmath::Vsqrt(nq_face,work,1,work,1);
737  Vmath::Sdiv (nq_face,1.0,work,1,work,1);
738 
739  for(i = 0; i < GetCoordim(); ++i)
740  {
741  Vmath::Vmul(nq_face,normal[i],1,work,1,normal[i],1);
742  }
743  }
744  }
745 
746  //---------------------------------------
747  // Matrix creation functions
748  //---------------------------------------
749 
751  {
752  DNekMatSharedPtr returnval;
753 
754  switch(mkey.GetMatrixType())
755  {
762  returnval = Expansion3D::v_GenMatrix(mkey);
763  break;
764  default:
765  returnval = StdPyrExp::v_GenMatrix(mkey);
766  }
767 
768  return returnval;
769  }
770 
772  {
773  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
774  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
775  LibUtilities::BasisKey bkey2 = m_base[2]->GetBasisKey();
777  MemoryManager<StdPyrExp>::AllocateSharedPtr(bkey0, bkey1, bkey2);
778 
779  return tmp->GetStdMatrix(mkey);
780  }
781 
783  {
784  return m_matrixManager[mkey];
785  }
786 
788  {
789  return m_staticCondMatrixManager[mkey];
790  }
791 
793  {
794  m_staticCondMatrixManager.DeleteObject(mkey);
795  }
796 
798  {
799  DNekScalMatSharedPtr returnval;
801 
802  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
803 
804  switch(mkey.GetMatrixType())
805  {
806  case StdRegions::eMass:
807  {
808  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
809  {
810  NekDouble one = 1.0;
811  DNekMatSharedPtr mat = GenMatrix(mkey);
813  }
814  else
815  {
816  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
817  DNekMatSharedPtr mat = GetStdMatrix(mkey);
819  }
820  }
821  break;
823  {
824  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
825  {
826  NekDouble one = 1.0;
828  *this);
829  DNekMatSharedPtr mat = GenMatrix(masskey);
830  mat->Invert();
832  }
833  else
834  {
835  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
836  DNekMatSharedPtr mat = GetStdMatrix(mkey);
838  }
839  }
840  break;
842  {
843  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
844  mkey.GetNVarCoeff() > 0 ||
845  mkey.ConstFactorExists(
847  {
848  NekDouble one = 1.0;
849  DNekMatSharedPtr mat = GenMatrix(mkey);
850 
852  }
853  else
854  {
856  mkey.GetShapeType(), *this);
858  mkey.GetShapeType(), *this);
860  mkey.GetShapeType(), *this);
862  mkey.GetShapeType(), *this);
864  mkey.GetShapeType(), *this);
866  mkey.GetShapeType(), *this);
867 
868  DNekMat &lap00 = *GetStdMatrix(lap00key);
869  DNekMat &lap01 = *GetStdMatrix(lap01key);
870  DNekMat &lap02 = *GetStdMatrix(lap02key);
871  DNekMat &lap11 = *GetStdMatrix(lap11key);
872  DNekMat &lap12 = *GetStdMatrix(lap12key);
873  DNekMat &lap22 = *GetStdMatrix(lap22key);
874 
875  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
877  m_metricinfo->GetGmat(ptsKeys);
878 
879  int rows = lap00.GetRows();
880  int cols = lap00.GetColumns();
881 
883  ::AllocateSharedPtr(rows,cols);
884 
885  (*lap) = gmat[0][0]*lap00
886  + gmat[4][0]*lap11
887  + gmat[8][0]*lap22
888  + gmat[3][0]*(lap01 + Transpose(lap01))
889  + gmat[6][0]*(lap02 + Transpose(lap02))
890  + gmat[7][0]*(lap12 + Transpose(lap12));
891 
892  returnval = MemoryManager<DNekScalMat>
893  ::AllocateSharedPtr(jac, lap);
894  }
895  }
896  break;
898  {
900  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
901  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
902  MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
903  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
904 
905  int rows = LapMat.GetRows();
906  int cols = LapMat.GetColumns();
907 
909 
910  (*helm) = LapMat + factor*MassMat;
911 
912  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, helm);
913  }
914  break;
915  default:
916  NEKERROR(ErrorUtil::efatal, "Matrix creation not defined");
917  break;
918  }
919 
920  return returnval;
921  }
922 
924  {
925  DNekScalBlkMatSharedPtr returnval;
926 
927  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
928 
929  // set up block matrix system
930  unsigned int nbdry = NumBndryCoeffs();
931  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
932  unsigned int exp_size[] = {nbdry, nint};
933  unsigned int nblks = 2;
934  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks, nblks, exp_size, exp_size); //Really need a constructor which takes Arrays
935  NekDouble factor = 1.0;
936 
937  switch(mkey.GetMatrixType())
938  {
940  case StdRegions::eHelmholtz: // special case since Helmholtz not defined in StdRegions
941 
942  // use Deformed case for both regular and deformed geometries
943  factor = 1.0;
944  goto UseLocRegionsMatrix;
945  break;
946  default:
947  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
948  {
949  factor = 1.0;
950  goto UseLocRegionsMatrix;
951  }
952  else
953  {
955  factor = mat->Scale();
956  goto UseStdRegionsMatrix;
957  }
958  break;
959  UseStdRegionsMatrix:
960  {
961  NekDouble invfactor = 1.0/factor;
962  NekDouble one = 1.0;
965  DNekMatSharedPtr Asubmat;
966 
967  //TODO: check below
968  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
969  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
970  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
971  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
972  }
973  break;
974  UseLocRegionsMatrix:
975  {
976  int i,j;
977  NekDouble invfactor = 1.0/factor;
978  NekDouble one = 1.0;
979  DNekScalMat &mat = *GetLocMatrix(mkey);
984 
985  Array<OneD,unsigned int> bmap(nbdry);
986  Array<OneD,unsigned int> imap(nint);
987  GetBoundaryMap(bmap);
988  GetInteriorMap(imap);
989 
990  for(i = 0; i < nbdry; ++i)
991  {
992  for(j = 0; j < nbdry; ++j)
993  {
994  (*A)(i,j) = mat(bmap[i],bmap[j]);
995  }
996 
997  for(j = 0; j < nint; ++j)
998  {
999  (*B)(i,j) = mat(bmap[i],imap[j]);
1000  }
1001  }
1002 
1003  for(i = 0; i < nint; ++i)
1004  {
1005  for(j = 0; j < nbdry; ++j)
1006  {
1007  (*C)(i,j) = mat(imap[i],bmap[j]);
1008  }
1009 
1010  for(j = 0; j < nint; ++j)
1011  {
1012  (*D)(i,j) = mat(imap[i],imap[j]);
1013  }
1014  }
1015 
1016  // Calculate static condensed system
1017  if(nint)
1018  {
1019  D->Invert();
1020  (*B) = (*B)*(*D);
1021  (*A) = (*A) - (*B)*(*C);
1022  }
1023 
1024  DNekScalMatSharedPtr Atmp;
1025 
1026  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
1027  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
1028  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
1029  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
1030 
1031  }
1032  }
1033  return returnval;
1034  }
1035 
1037  {
1038  if (m_metrics.count(eMetricQuadrature) == 0)
1039  {
1041  }
1042 
1043  int i, j;
1044  const unsigned int nqtot = GetTotPoints();
1045  const unsigned int dim = 3;
1046  const MetricType m[3][3] = {
1050  };
1051 
1052  for (unsigned int i = 0; i < dim; ++i)
1053  {
1054  for (unsigned int j = i; j < dim; ++j)
1055  {
1056  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1057  }
1058  }
1059 
1060  // Define shorthand synonyms for m_metrics storage
1061  Array<OneD,NekDouble> g0 (m_metrics[m[0][0]]);
1062  Array<OneD,NekDouble> g1 (m_metrics[m[1][1]]);
1063  Array<OneD,NekDouble> g2 (m_metrics[m[2][2]]);
1064  Array<OneD,NekDouble> g3 (m_metrics[m[0][1]]);
1065  Array<OneD,NekDouble> g4 (m_metrics[m[0][2]]);
1066  Array<OneD,NekDouble> g5 (m_metrics[m[1][2]]);
1067 
1068  // Allocate temporary storage
1069  Array<OneD,NekDouble> alloc(9*nqtot,0.0);
1070  Array<OneD,NekDouble> h0 (nqtot, alloc);
1071  Array<OneD,NekDouble> h1 (nqtot, alloc+ 1*nqtot);
1072  Array<OneD,NekDouble> h2 (nqtot, alloc+ 2*nqtot);
1073  Array<OneD,NekDouble> wsp1 (nqtot, alloc+ 3*nqtot);
1074  Array<OneD,NekDouble> wsp2 (nqtot, alloc+ 4*nqtot);
1075  Array<OneD,NekDouble> wsp3 (nqtot, alloc+ 5*nqtot);
1076  Array<OneD,NekDouble> wsp4 (nqtot, alloc+ 6*nqtot);
1077  Array<OneD,NekDouble> wsp5 (nqtot, alloc+ 7*nqtot);
1078  Array<OneD,NekDouble> wsp6 (nqtot, alloc+ 8*nqtot);
1079 
1080  const Array<TwoD, const NekDouble>& df =
1081  m_metricinfo->GetDerivFactors(GetPointsKeys());
1082  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
1083  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
1084  const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
1085  const unsigned int nquad0 = m_base[0]->GetNumPoints();
1086  const unsigned int nquad1 = m_base[1]->GetNumPoints();
1087  const unsigned int nquad2 = m_base[2]->GetNumPoints();
1088 
1089  // Populate collapsed coordinate arrays h0, h1 and h2.
1090  for(j = 0; j < nquad2; ++j)
1091  {
1092  for(i = 0; i < nquad1; ++i)
1093  {
1094  Vmath::Fill(nquad0, 2.0/(1.0-z2[j]), &h0[0]+i*nquad0 + j*nquad0*nquad1,1);
1095  Vmath::Fill(nquad0, 1.0/(1.0-z2[j]), &h1[0]+i*nquad0 + j*nquad0*nquad1,1);
1096  Vmath::Fill(nquad0, (1.0+z1[i])/(1.0-z2[j]), &h2[0]+i*nquad0 + j*nquad0*nquad1,1);
1097  }
1098  }
1099  for(i = 0; i < nquad0; i++)
1100  {
1101  Blas::Dscal(nquad1*nquad2, 1+z0[i], &h1[0]+i, nquad0);
1102  }
1103 
1104  // Step 3. Construct combined metric terms for physical space to
1105  // collapsed coordinate system.
1106  // Order of construction optimised to minimise temporary storage
1107  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1108  {
1109  // f_{1k}
1110  Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &df[2][0], 1, &h1[0], 1, &wsp1[0], 1);
1111  Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &df[5][0], 1, &h1[0], 1, &wsp2[0], 1);
1112  Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &df[8][0], 1, &h1[0], 1, &wsp3[0], 1);
1113 
1114  // g0
1115  Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp1[0], 1, &wsp2[0], 1, &wsp2[0], 1, &g0[0], 1);
1116  Vmath::Vvtvp (nqtot, &wsp3[0], 1, &wsp3[0], 1, &g0[0], 1, &g0[0], 1);
1117 
1118  // g4
1119  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp1[0], 1, &df[5][0], 1, &wsp2[0], 1, &g4[0], 1);
1120  Vmath::Vvtvp (nqtot, &df[8][0], 1, &wsp3[0], 1, &g4[0], 1, &g4[0], 1);
1121 
1122  // f_{2k}
1123  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &h0[0], 1, &df[2][0], 1, &h2[0], 1, &wsp4[0], 1);
1124  Vmath::Vvtvvtp(nqtot, &df[4][0], 1, &h0[0], 1, &df[5][0], 1, &h2[0], 1, &wsp5[0], 1);
1125  Vmath::Vvtvvtp(nqtot, &df[7][0], 1, &h0[0], 1, &df[8][0], 1, &h2[0], 1, &wsp6[0], 1);
1126 
1127  // g1
1128  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0], 1, &g1[0], 1);
1129  Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp6[0], 1, &g1[0], 1, &g1[0], 1);
1130 
1131  // g3
1132  Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp4[0], 1, &wsp2[0], 1, &wsp5[0], 1, &g3[0], 1);
1133  Vmath::Vvtvp (nqtot, &wsp3[0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1134 
1135  // g5
1136  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0], 1, &g5[0], 1);
1137  Vmath::Vvtvp (nqtot, &df[8][0], 1, &wsp6[0], 1, &g5[0], 1, &g5[0], 1);
1138 
1139  // g2
1140  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &df[2][0], 1, &df[5][0], 1, &df[5][0], 1, &g2[0], 1);
1141  Vmath::Vvtvp (nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1142  }
1143  else
1144  {
1145  // f_{1k}
1146  Vmath::Svtsvtp(nqtot, df[0][0], &h0[0], 1, df[2][0], &h1[0], 1, &wsp1[0], 1);
1147  Vmath::Svtsvtp(nqtot, df[3][0], &h0[0], 1, df[5][0], &h1[0], 1, &wsp2[0], 1);
1148  Vmath::Svtsvtp(nqtot, df[6][0], &h0[0], 1, df[8][0], &h1[0], 1, &wsp3[0], 1);
1149 
1150  // g0
1151  Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp1[0], 1, &wsp2[0], 1, &wsp2[0], 1, &g0[0], 1);
1152  Vmath::Vvtvp (nqtot, &wsp3[0], 1, &wsp3[0], 1, &g0[0], 1, &g0[0], 1);
1153 
1154  // g4
1155  Vmath::Svtsvtp(nqtot, df[2][0], &wsp1[0], 1, df[5][0], &wsp2[0], 1, &g4[0], 1);
1156  Vmath::Svtvp (nqtot, df[8][0], &wsp3[0], 1, &g4[0], 1, &g4[0], 1);
1157 
1158  // f_{2k}
1159  Vmath::Svtsvtp(nqtot, df[1][0], &h0[0], 1, df[2][0], &h2[0], 1, &wsp4[0], 1);
1160  Vmath::Svtsvtp(nqtot, df[4][0], &h0[0], 1, df[5][0], &h2[0], 1, &wsp5[0], 1);
1161  Vmath::Svtsvtp(nqtot, df[7][0], &h0[0], 1, df[8][0], &h2[0], 1, &wsp6[0], 1);
1162 
1163  // g1
1164  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0], 1, &g1[0], 1);
1165  Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp6[0], 1, &g1[0], 1, &g1[0], 1);
1166 
1167  // g3
1168  Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp4[0], 1, &wsp2[0], 1, &wsp5[0], 1, &g3[0], 1);
1169  Vmath::Vvtvp (nqtot, &wsp3[0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1170 
1171  // g5
1172  Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1, &g5[0], 1);
1173  Vmath::Svtvp (nqtot, df[8][0], &wsp6[0], 1, &g5[0], 1, &g5[0], 1);
1174 
1175  // g2
1176  Vmath::Fill(nqtot, df[2][0]*df[2][0] + df[5][0]*df[5][0] + df[8][0]*df[8][0], &g2[0], 1);
1177  }
1178 
1179  for (unsigned int i = 0; i < dim; ++i)
1180  {
1181  for (unsigned int j = i; j < dim; ++j)
1182  {
1184  m_metrics[m[i][j]]);
1185 
1186  }
1187  }
1188  }
1189 
1191  const Array<OneD, const NekDouble> &inarray,
1192  Array<OneD, NekDouble> &outarray,
1194  {
1195  // This implementation is only valid when there are no coefficients
1196  // associated to the Laplacian operator
1197  if (m_metrics.count(eMetricLaplacian00) == 0)
1198  {
1200  }
1201 
1202  int nquad0 = m_base[0]->GetNumPoints();
1203  int nquad1 = m_base[1]->GetNumPoints();
1204  int nq2 = m_base[2]->GetNumPoints();
1205  int nqtot = nquad0*nquad1*nq2;
1206 
1207  ASSERTL1(wsp.num_elements() >= 6*nqtot,
1208  "Insufficient workspace size.");
1209  ASSERTL1(m_ncoeffs <= nqtot,
1210  "Workspace not set up for ncoeffs > nqtot");
1211 
1212  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
1213  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
1214  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata();
1215  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
1216  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
1217  const Array<OneD, const NekDouble>& dbase2 = m_base[2]->GetDbdata();
1224 
1225  // Allocate temporary storage
1226  Array<OneD,NekDouble> wsp0 (2*nqtot, wsp);
1227  Array<OneD,NekDouble> wsp1 ( nqtot, wsp+1*nqtot);
1228  Array<OneD,NekDouble> wsp2 ( nqtot, wsp+2*nqtot);
1229  Array<OneD,NekDouble> wsp3 ( nqtot, wsp+3*nqtot);
1230  Array<OneD,NekDouble> wsp4 ( nqtot, wsp+4*nqtot);
1231  Array<OneD,NekDouble> wsp5 ( nqtot, wsp+5*nqtot);
1232 
1233  // LAPLACIAN MATRIX OPERATION
1234  // wsp1 = du_dxi1 = D_xi1 * inarray = D_xi1 * u
1235  // wsp2 = du_dxi2 = D_xi2 * inarray = D_xi2 * u
1236  // wsp2 = du_dxi3 = D_xi3 * inarray = D_xi3 * u
1237  StdExpansion3D::PhysTensorDeriv(inarray,wsp0,wsp1,wsp2);
1238 
1239  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1240  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1241  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1242  // especially for this purpose
1243  Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp0[0],1,&metric01[0],1,&wsp1[0],1,&wsp3[0],1);
1244  Vmath::Vvtvp (nqtot,&metric02[0],1,&wsp2[0],1,&wsp3[0],1,&wsp3[0],1);
1245  Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp0[0],1,&metric11[0],1,&wsp1[0],1,&wsp4[0],1);
1246  Vmath::Vvtvp (nqtot,&metric12[0],1,&wsp2[0],1,&wsp4[0],1,&wsp4[0],1);
1247  Vmath::Vvtvvtp(nqtot,&metric02[0],1,&wsp0[0],1,&metric12[0],1,&wsp1[0],1,&wsp5[0],1);
1248  Vmath::Vvtvp (nqtot,&metric22[0],1,&wsp2[0],1,&wsp5[0],1,&wsp5[0],1);
1249 
1250  // outarray = m = (D_xi1 * B)^T * k
1251  // wsp1 = n = (D_xi2 * B)^T * l
1252  IProductWRTBase_SumFacKernel(dbase0,base1,base2,wsp3,outarray,wsp0,false,true,true);
1253  IProductWRTBase_SumFacKernel(base0,dbase1,base2,wsp4,wsp2, wsp0,true,false,true);
1254  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
1255  IProductWRTBase_SumFacKernel(base0,base1,dbase2,wsp5,wsp2, wsp0,true,true,false);
1256  Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
1257  }
1258  }//end of namespace
1259 }//end of namespace
const LibUtilities::PointsKeyVector GetPointsKeys() const
boost::shared_ptr< PyrGeom > PyrGeomSharedPtr
Definition: PyrGeom.h:84
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:470
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:142
DNekScalBlkMatSharedPtr CreateStaticCondMatrix(const MatrixKey &mkey)
Definition: PyrExp.cpp:923
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:191
const VarCoeffMap & GetVarCoeffs() const
Definition: StdMatrixKey.h:168
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
void v_ComputeFaceNormal(const int face)
Definition: PyrExp.cpp:478
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: PyrExp.h:155
virtual void v_ComputeLaplacianMetric()
Definition: PyrExp.cpp:1036
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:408
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:947
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition: PyrExp.cpp:797
virtual void v_GetFacePhysMap(const int face, Array< OneD, int > &outarray)
Definition: PyrExp.cpp:378
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: PyrExp.cpp:771
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Calculate the inner product of inarray with respect to the basis B=base0*base1*base2 and put into out...
Definition: PyrExp.cpp:276
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: PyrExp.cpp:750
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
Definition: PyrExp.cpp:304
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:485
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:442
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:131
virtual int v_GetCoordim()
Definition: PyrExp.cpp:373
STL namespace.
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:271
LibUtilities::ShapeType GetShapeType() const
Definition: StdMatrixKey.h:87
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:130
virtual void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
Definition: PyrExp.cpp:1190
boost::shared_ptr< StdPyrExp > StdPyrExpSharedPtr
Definition: StdPyrExp.h:263
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
Definition: PyrExp.cpp:329
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:706
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:131
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrate the physical point list inarray over pyramidic region and return the value.
Definition: PyrExp.cpp:110
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis...
Definition: Interp.cpp:116
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Forward transform from physical quadrature space stored in inarray and evaluate the expansion coeffic...
Definition: PyrExp.cpp:219
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:711
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:213
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: Expansion.cpp:224
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:128
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:819
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Calculate the derivative of the physical points.
Definition: PyrExp.cpp:137
std::map< int, NormalVector > m_faceNormals
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:165
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: PyrExp.cpp:792
boost::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:72
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:85
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:145
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:161
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:537
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:250
Geometry is straight-sided with constant geometric factors.
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:131
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
vvtvvtp (scalar times vector plus scalar times vector):
Definition: Vmath.cpp:591
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: PyrExp.cpp:782
const LibUtilities::BasisKey DetFaceBasisKey(const int i, const int k) const
Definition: StdExpansion.h:324
GeomType
Indicates the type of element geometry.
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: PyrExp.cpp:787
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: PyrExp.cpp:347
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
virtual NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coord, const Array< OneD, const NekDouble > &physvals)
This function evaluates the expansion at a single (arbitrary) point of the domain.
Definition: PyrExp.cpp:355
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
Array< OneD, LibUtilities::BasisSharedPtr > m_base
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
Geometry is curved or has non-constant factors.
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:814
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: PyrExp.h:154
Describes the specification for a Basis.
Definition: Basis.h:50
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
Definition: PyrExp.cpp:312
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:299
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:183
PyrExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc, const SpatialDomains::PyrGeomSharedPtr &geom)
Constructor using BasisKey class for quadrature points and order definition.
Definition: PyrExp.cpp:44