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