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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: PyrExp routines
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
36 #include <LocalRegions/PyrExp.h>
37 
38 using namespace std;
39 
40 namespace Nektar
41 {
42 namespace LocalRegions
43 {
44 
45 PyrExp::PyrExp(const LibUtilities::BasisKey &Ba,
46  const LibUtilities::BasisKey &Bb,
47  const LibUtilities::BasisKey &Bc,
49  : StdExpansion(LibUtilities::StdPyrData::getNumberOfCoefficients(
50  Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
51  3, Ba, Bb, Bc),
52  StdExpansion3D(LibUtilities::StdPyrData::getNumberOfCoefficients(
53  Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
54  Ba, Bb, Bc),
55  StdPyrExp(Ba, Bb, Bc), Expansion(geom), Expansion3D(geom),
56  m_matrixManager(
57  std::bind(&Expansion3D::CreateMatrix, this, std::placeholders::_1),
58  std::string("PyrExpMatrix")),
59  m_staticCondMatrixManager(std::bind(&Expansion::CreateStaticCondMatrix,
60  this, std::placeholders::_1),
61  std::string("PyrExpStaticCondMatrix"))
62 {
63 }
64 
66  : StdExpansion(T), StdExpansion3D(T), StdPyrExp(T), Expansion(T),
67  Expansion3D(T), m_matrixManager(T.m_matrixManager),
68  m_staticCondMatrixManager(T.m_staticCondMatrixManager)
69 {
70 }
71 
73 {
74 }
75 
76 //----------------------------
77 // Integration Methods
78 //----------------------------
79 
80 /**
81  * \brief Integrate the physical point list \a inarray over pyramidic
82  * region and return the value.
83  *
84  * Inputs:\n
85  *
86  * - \a inarray: definition of function to be returned at quadrature
87  * point of expansion.
88  *
89  * Outputs:\n
90  *
91  * - returns \f$\int^1_{-1}\int^1_{-1}\int^1_{-1} u(\bar \eta_1,
92  * \eta_2, \eta_3) J[i,j,k] d \bar \eta_1 d \eta_2 d \eta_3\f$ \n \f$=
93  * \sum_{i=0}^{Q_1 - 1} \sum_{j=0}^{Q_2 - 1} \sum_{k=0}^{Q_3 - 1}
94  * u(\bar \eta_{1i}^{0,0}, \eta_{2j}^{0,0},\eta_{3k}^{2,0})w_{i}^{0,0}
95  * w_{j}^{0,0} \hat w_{k}^{2,0} \f$ \n where \f$inarray[i,j, k] =
96  * u(\bar \eta_{1i},\eta_{2j}, \eta_{3k}) \f$, \n \f$\hat w_{k}^{2,0}
97  * = \frac {w^{2,0}} {2} \f$ \n and \f$ J[i,j,k] \f$ is the Jacobian
98  * evaluated at the quadrature point.
99  */
101 {
102  int nquad0 = m_base[0]->GetNumPoints();
103  int nquad1 = m_base[1]->GetNumPoints();
104  int nquad2 = m_base[2]->GetNumPoints();
106  Array<OneD, NekDouble> tmp(nquad0 * nquad1 * nquad2);
107 
108  // multiply inarray with Jacobian
109  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
110  {
111  Vmath::Vmul(nquad0 * nquad1 * nquad2, &jac[0], 1,
112  (NekDouble *)&inarray[0], 1, &tmp[0], 1);
113  }
114  else
115  {
116  Vmath::Smul(nquad0 * nquad1 * nquad2, (NekDouble)jac[0],
117  (NekDouble *)&inarray[0], 1, &tmp[0], 1);
118  }
119 
120  // call StdPyrExp version;
121  return StdPyrExp::v_Integral(tmp);
122 }
123 
124 //----------------------------
125 // Differentiation Methods
126 //----------------------------
127 
129  Array<OneD, NekDouble> &out_d0,
130  Array<OneD, NekDouble> &out_d1,
131  Array<OneD, NekDouble> &out_d2)
132 {
133  int nquad0 = m_base[0]->GetNumPoints();
134  int nquad1 = m_base[1]->GetNumPoints();
135  int nquad2 = m_base[2]->GetNumPoints();
137  m_metricinfo->GetDerivFactors(GetPointsKeys());
138  Array<OneD, NekDouble> diff0(nquad0 * nquad1 * nquad2);
139  Array<OneD, NekDouble> diff1(nquad0 * nquad1 * nquad2);
140  Array<OneD, NekDouble> diff2(nquad0 * nquad1 * nquad2);
141 
142  StdPyrExp::v_PhysDeriv(inarray, diff0, diff1, diff2);
143 
144  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
145  {
146  if (out_d0.size())
147  {
148  Vmath::Vmul(nquad0 * nquad1 * nquad2, &gmat[0][0], 1, &diff0[0], 1,
149  &out_d0[0], 1);
150  Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[1][0], 1, &diff1[0], 1,
151  &out_d0[0], 1, &out_d0[0], 1);
152  Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[2][0], 1, &diff2[0], 1,
153  &out_d0[0], 1, &out_d0[0], 1);
154  }
155 
156  if (out_d1.size())
157  {
158  Vmath::Vmul(nquad0 * nquad1 * nquad2, &gmat[3][0], 1, &diff0[0], 1,
159  &out_d1[0], 1);
160  Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[4][0], 1, &diff1[0], 1,
161  &out_d1[0], 1, &out_d1[0], 1);
162  Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[5][0], 1, &diff2[0], 1,
163  &out_d1[0], 1, &out_d1[0], 1);
164  }
165 
166  if (out_d2.size())
167  {
168  Vmath::Vmul(nquad0 * nquad1 * nquad2, &gmat[6][0], 1, &diff0[0], 1,
169  &out_d2[0], 1);
170  Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[7][0], 1, &diff1[0], 1,
171  &out_d2[0], 1, &out_d2[0], 1);
172  Vmath::Vvtvp(nquad0 * nquad1 * nquad2, &gmat[8][0], 1, &diff2[0], 1,
173  &out_d2[0], 1, &out_d2[0], 1);
174  }
175  }
176  else // regular geometry
177  {
178  if (out_d0.size())
179  {
180  Vmath::Smul(nquad0 * nquad1 * nquad2, gmat[0][0], &diff0[0], 1,
181  &out_d0[0], 1);
182  Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[1][0], &diff1[0], 1,
183  &out_d0[0], 1);
184  Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[2][0], &diff2[0], 1,
185  &out_d0[0], 1);
186  }
187 
188  if (out_d1.size())
189  {
190  Vmath::Smul(nquad0 * nquad1 * nquad2, gmat[3][0], &diff0[0], 1,
191  &out_d1[0], 1);
192  Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[4][0], &diff1[0], 1,
193  &out_d1[0], 1);
194  Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[5][0], &diff2[0], 1,
195  &out_d1[0], 1);
196  }
197 
198  if (out_d2.size())
199  {
200  Vmath::Smul(nquad0 * nquad1 * nquad2, gmat[6][0], &diff0[0], 1,
201  &out_d2[0], 1);
202  Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[7][0], &diff1[0], 1,
203  &out_d2[0], 1);
204  Blas::Daxpy(nquad0 * nquad1 * nquad2, gmat[8][0], &diff2[0], 1,
205  &out_d2[0], 1);
206  }
207  }
208 }
209 
210 //---------------------------------------
211 // Transforms
212 //---------------------------------------
213 
214 /**
215  * \brief Forward transform from physical quadrature space stored in
216  * \a inarray and evaluate the expansion coefficients and store in \a
217  * (this)->m_coeffs
218  *
219  * Inputs:\n
220  *
221  * - \a inarray: array of physical quadrature points to be transformed
222  *
223  * Outputs:\n
224  *
225  * - (this)->_coeffs: updated array of expansion coefficients.
226  */
228  Array<OneD, NekDouble> &outarray)
229 {
230  if (m_base[0]->Collocation() && m_base[1]->Collocation() &&
231  m_base[2]->Collocation())
232  {
233  Vmath::Vcopy(GetNcoeffs(), &inarray[0], 1, &outarray[0], 1);
234  }
235  else
236  {
237  v_IProductWRTBase(inarray, outarray);
238 
239  // get Mass matrix inverse
240  MatrixKey masskey(StdRegions::eInvMass, DetShapeType(), *this);
241  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
242 
243  // copy inarray in case inarray == outarray
244  DNekVec in(m_ncoeffs, outarray);
245  DNekVec out(m_ncoeffs, outarray, eWrapper);
246 
247  out = (*matsys) * in;
248  }
249 }
250 
251 //---------------------------------------
252 // Inner product functions
253 //---------------------------------------
254 
255 /**
256  * \brief Calculate the inner product of inarray with respect to the
257  * basis B=base0*base1*base2 and put into outarray:
258  *
259  * \f$ \begin{array}{rcl} I_{pqr} = (\phi_{pqr}, u)_{\delta} & = &
260  * \sum_{i=0}^{nq_0} \sum_{j=0}^{nq_1} \sum_{k=0}^{nq_2} \psi_{p}^{a}
261  * (\bar \eta_{1i}) \psi_{q}^{a} (\eta_{2j}) \psi_{pqr}^{c}
262  * (\eta_{3k}) w_i w_j w_k u(\bar \eta_{1,i} \eta_{2,j} \eta_{3,k})
263  * J_{i,j,k}\\ & = & \sum_{i=0}^{nq_0} \psi_p^a(\bar \eta_{1,i})
264  * \sum_{j=0}^{nq_1} \psi_{q}^a(\eta_{2,j}) \sum_{k=0}^{nq_2}
265  * \psi_{pqr}^c u(\bar \eta_{1i},\eta_{2j},\eta_{3k}) J_{i,j,k}
266  * \end{array} \f$ \n
267  *
268  * where
269  *
270  * \f$\phi_{pqr} (\xi_1 , \xi_2 , \xi_3) = \psi_p^a (\bar \eta_1)
271  * \psi_{q}^a (\eta_2) \psi_{pqr}^c (\eta_3) \f$ \n
272  *
273  * which can be implemented as \n \f$f_{pqr} (\xi_{3k}) =
274  * \sum_{k=0}^{nq_3} \psi_{pqr}^c u(\bar
275  * \eta_{1i},\eta_{2j},\eta_{3k}) J_{i,j,k} = {\bf B_3 U} \f$ \n \f$
276  * g_{pq} (\xi_{3k}) = \sum_{j=0}^{nq_1} \psi_{q}^a (\xi_{2j}) f_{pqr}
277  * (\xi_{3k}) = {\bf B_2 F} \f$ \n \f$ (\phi_{pqr}, u)_{\delta} =
278  * \sum_{k=0}^{nq_0} \psi_{p}^a (\xi_{3k}) g_{pq} (\xi_{3k}) = {\bf
279  * B_1 G} \f$
280  */
281 
283  Array<OneD, NekDouble> &outarray)
284 {
285  v_IProductWRTBase_SumFac(inarray, outarray);
286 }
287 
289  const Array<OneD, const NekDouble> &inarray,
290  Array<OneD, NekDouble> &outarray, bool multiplybyweights)
291 {
292  const int nquad0 = m_base[0]->GetNumPoints();
293  const int nquad1 = m_base[1]->GetNumPoints();
294  const int nquad2 = m_base[2]->GetNumPoints();
295  const int order0 = m_base[0]->GetNumModes();
296  const int order1 = m_base[1]->GetNumModes();
297 
298  Array<OneD, NekDouble> wsp(order0 * nquad2 * (nquad1 + order1));
299 
300  if (multiplybyweights)
301  {
302  Array<OneD, NekDouble> tmp(nquad0 * nquad1 * nquad2);
303 
304  MultiplyByQuadratureMetric(inarray, tmp);
305 
307  m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
308  tmp, outarray, wsp, true, true, true);
309  }
310  else
311  {
313  m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
314  inarray, outarray, wsp, true, true, true);
315  }
316 }
317 
318 /**
319  * @brief Calculates the inner product \f$ I_{pqr} = (u,
320  * \partial_{x_i} \phi_{pqr}) \f$.
321  *
322  * The derivative of the basis functions is performed using the chain
323  * rule in order to incorporate the geometric factors. Assuming that
324  * the basis functions are a tensor product
325  * \f$\phi_{pqr}(\eta_1,\eta_2,\eta_3) =
326  * \phi_1(\eta_1)\phi_2(\eta_2)\phi_3(\eta_3)\f$, this yields the
327  * result
328  *
329  * \f[
330  * I_{pqr} = \sum_{j=1}^3 \left(u, \frac{\partial u}{\partial \eta_j}
331  * \frac{\partial \eta_j}{\partial x_i}\right)
332  * \f]
333  *
334  * In the pyramid element, we must also incorporate a second set
335  * of geometric factors which incorporate the collapsed co-ordinate
336  * system, so that
337  *
338  * \f[ \frac{\partial\eta_j}{\partial x_i} = \sum_{k=1}^3
339  * \frac{\partial\eta_j}{\partial\xi_k}\frac{\partial\xi_k}{\partial
340  * x_i} \f]
341  *
342  * These derivatives can be found on p152 of Sherwin & Karniadakis.
343  *
344  * @param dir Direction in which to take the derivative.
345  * @param inarray The function \f$ u \f$.
346  * @param outarray Value of the inner product.
347  */
349  const Array<OneD, const NekDouble> &inarray,
350  Array<OneD, NekDouble> &outarray)
351 {
352  v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
353 }
354 
356  const int dir, const Array<OneD, const NekDouble> &inarray,
357  Array<OneD, NekDouble> &outarray)
358 {
359  const int nquad0 = m_base[0]->GetNumPoints();
360  const int nquad1 = m_base[1]->GetNumPoints();
361  const int nquad2 = m_base[2]->GetNumPoints();
362  const int order0 = m_base[0]->GetNumModes();
363  const int order1 = m_base[1]->GetNumModes();
364  const int nqtot = nquad0 * nquad1 * nquad2;
365 
366  Array<OneD, NekDouble> tmp1(nqtot);
367  Array<OneD, NekDouble> tmp2(nqtot);
368  Array<OneD, NekDouble> tmp3(nqtot);
369  Array<OneD, NekDouble> tmp4(nqtot);
372  std::max(nqtot, order0 * nquad2 * (nquad1 + order1)));
373 
374  MultiplyByQuadratureMetric(inarray, tmp1);
375 
377  tmp2D[0] = tmp2;
378  tmp2D[1] = tmp3;
379  tmp2D[2] = tmp4;
380 
381  PyrExp::v_AlignVectorToCollapsedDir(dir, tmp1, tmp2D);
382 
383  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
384  m_base[2]->GetBdata(), tmp2, outarray, wsp,
385  false, true, true);
386 
387  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
388  m_base[2]->GetBdata(), tmp3, tmp6, wsp, true,
389  false, true);
390 
391  Vmath::Vadd(m_ncoeffs, tmp6, 1, outarray, 1, outarray, 1);
392 
393  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(),
394  m_base[2]->GetDbdata(), tmp4, tmp6, wsp, true,
395  true, false);
396 
397  Vmath::Vadd(m_ncoeffs, tmp6, 1, outarray, 1, outarray, 1);
398 }
399 
401  const int dir, const Array<OneD, const NekDouble> &inarray,
402  Array<OneD, Array<OneD, NekDouble>> &outarray)
403 {
404  const int nquad0 = m_base[0]->GetNumPoints();
405  const int nquad1 = m_base[1]->GetNumPoints();
406  const int nquad2 = m_base[2]->GetNumPoints();
407  const int order0 = m_base[0]->GetNumModes();
408  const int order1 = m_base[1]->GetNumModes();
409  const int nqtot = nquad0 * nquad1 * nquad2;
410 
411  const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
412  const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
413  const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
414 
415  Array<OneD, NekDouble> gfac0(nquad0);
416  Array<OneD, NekDouble> gfac1(nquad1);
417  Array<OneD, NekDouble> gfac2(nquad2);
418  Array<OneD, NekDouble> tmp5(nqtot);
420  std::max(nqtot, order0 * nquad2 * (nquad1 + order1)));
421 
422  Array<OneD, NekDouble> tmp2 = outarray[0];
423  Array<OneD, NekDouble> tmp3 = outarray[1];
424  Array<OneD, NekDouble> tmp4 = outarray[2];
425 
426  const Array<TwoD, const NekDouble> &df =
427  m_metricinfo->GetDerivFactors(GetPointsKeys());
428 
430  tmp1 = inarray;
431 
432  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
433  {
434  Vmath::Vmul(nqtot, &df[3 * dir][0], 1, tmp1.get(), 1, tmp2.get(), 1);
435  Vmath::Vmul(nqtot, &df[3 * dir + 1][0], 1, tmp1.get(), 1, tmp3.get(),
436  1);
437  Vmath::Vmul(nqtot, &df[3 * dir + 2][0], 1, tmp1.get(), 1, tmp4.get(),
438  1);
439  }
440  else
441  {
442  Vmath::Smul(nqtot, df[3 * dir][0], tmp1.get(), 1, tmp2.get(), 1);
443  Vmath::Smul(nqtot, df[3 * dir + 1][0], tmp1.get(), 1, tmp3.get(), 1);
444  Vmath::Smul(nqtot, df[3 * dir + 2][0], tmp1.get(), 1, tmp4.get(), 1);
445  }
446 
447  // set up geometric factor: (1+z0)/2
448  for (int i = 0; i < nquad0; ++i)
449  {
450  gfac0[i] = 0.5 * (1 + z0[i]);
451  }
452 
453  // set up geometric factor: (1+z1)/2
454  for (int i = 0; i < nquad1; ++i)
455  {
456  gfac1[i] = 0.5 * (1 + z1[i]);
457  }
458 
459  // Set up geometric factor: 2/(1-z2)
460  for (int i = 0; i < nquad2; ++i)
461  {
462  gfac2[i] = 2.0 / (1 - z2[i]);
463  }
464 
465  const int nq01 = nquad0 * nquad1;
466 
467  for (int i = 0; i < nquad2; ++i)
468  {
469  Vmath::Smul(nq01, gfac2[i], &tmp2[0] + i * nq01, 1, &tmp2[0] + i * nq01,
470  1); // 2/(1-z2) for d/dxi_0
471  Vmath::Smul(nq01, gfac2[i], &tmp3[0] + i * nq01, 1, &tmp3[0] + i * nq01,
472  1); // 2/(1-z2) for d/dxi_1
473  Vmath::Smul(nq01, gfac2[i], &tmp4[0] + i * nq01, 1, &tmp5[0] + i * nq01,
474  1); // 2/(1-z2) for d/dxi_2
475  }
476 
477  // (1+z0)/(1-z2) for d/d eta_0
478  for (int i = 0; i < nquad1 * nquad2; ++i)
479  {
480  Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp5[0] + i * nquad0, 1,
481  &wsp[0] + i * nquad0, 1);
482  }
483 
484  Vmath::Vadd(nqtot, &tmp2[0], 1, &wsp[0], 1, &tmp2[0], 1);
485 
486  // (1+z1)/(1-z2) for d/d eta_1
487  for (int i = 0; i < nquad1 * nquad2; ++i)
488  {
489  Vmath::Smul(nquad0, gfac1[i % nquad1], &tmp5[0] + i * nquad0, 1,
490  &tmp5[0] + i * nquad0, 1);
491  }
492  Vmath::Vadd(nqtot, &tmp3[0], 1, &tmp5[0], 1, &tmp3[0], 1);
493 }
494 
495 //---------------------------------------
496 // Evaluation functions
497 //---------------------------------------
498 
500 {
502  m_base[0]->GetBasisKey(), m_base[1]->GetBasisKey(),
503  m_base[2]->GetBasisKey());
504 }
505 
507 {
509  m_base[0]->GetPointsKey());
511  m_base[1]->GetPointsKey());
513  m_base[2]->GetPointsKey());
514 
516  bkey2);
517 }
518 
519 /*
520  * @brief Get the coordinates #coords at the local coordinates
521  * #Lcoords
522  */
524  Array<OneD, NekDouble> &coords)
525 {
526  int i;
527 
528  ASSERTL1(Lcoords[0] <= -1.0 && Lcoords[0] >= 1.0 && Lcoords[1] <= -1.0 &&
529  Lcoords[1] >= 1.0 && Lcoords[2] <= -1.0 && Lcoords[2] >= 1.0,
530  "Local coordinates are not in region [-1,1]");
531 
532  // m_geom->FillGeom(); // TODO: implement FillGeom()
533 
534  for (i = 0; i < m_geom->GetCoordim(); ++i)
535  {
536  coords[i] = m_geom->GetCoord(i, Lcoords);
537  }
538 }
539 
541  Array<OneD, NekDouble> &coords_2,
542  Array<OneD, NekDouble> &coords_3)
543 {
544  Expansion::v_GetCoords(coords_1, coords_2, coords_3);
545 }
546 
548  const NekDouble *data, const std::vector<unsigned int> &nummodes,
549  const int mode_offset, NekDouble *coeffs,
550  std::vector<LibUtilities::BasisType> &fromType)
551 {
552  int data_order0 = nummodes[mode_offset];
553  int fillorder0 = min(m_base[0]->GetNumModes(), data_order0);
554  int data_order1 = nummodes[mode_offset + 1];
555  int order1 = m_base[1]->GetNumModes();
556  int fillorder1 = min(order1, data_order1);
557  int data_order2 = nummodes[mode_offset + 2];
558  int order2 = m_base[2]->GetNumModes();
559  int fillorder2 = min(order2, data_order2);
560 
561  // Check if not same order or basis and if not make temp
562  // element to read in data
563  if (fromType[0] != m_base[0]->GetBasisType() ||
564  fromType[1] != m_base[1]->GetBasisType() ||
565  fromType[2] != m_base[2]->GetBasisType() || data_order0 != fillorder0 ||
566  data_order1 != fillorder1 || data_order2 != fillorder2)
567  {
568  // Construct a pyr with the appropriate basis type at our
569  // quadrature points, and one more to do a forwards
570  // transform. We can then copy the output to coeffs.
571  StdRegions::StdPyrExp tmpPyr(
572  LibUtilities::BasisKey(fromType[0], data_order0,
573  m_base[0]->GetPointsKey()),
574  LibUtilities::BasisKey(fromType[1], data_order1,
575  m_base[1]->GetPointsKey()),
576  LibUtilities::BasisKey(fromType[2], data_order2,
577  m_base[2]->GetPointsKey()));
578 
579  StdRegions::StdPyrExp tmpPyr2(m_base[0]->GetBasisKey(),
580  m_base[1]->GetBasisKey(),
581  m_base[2]->GetBasisKey());
582 
583  Array<OneD, const NekDouble> tmpData(tmpPyr.GetNcoeffs(), data);
584  Array<OneD, NekDouble> tmpBwd(tmpPyr2.GetTotPoints());
585  Array<OneD, NekDouble> tmpOut(tmpPyr2.GetNcoeffs());
586 
587  tmpPyr.BwdTrans(tmpData, tmpBwd);
588  tmpPyr2.FwdTrans(tmpBwd, tmpOut);
589  Vmath::Vcopy(tmpOut.size(), &tmpOut[0], 1, coeffs, 1);
590  }
591  else
592  {
593  Vmath::Vcopy(m_ncoeffs, &data[0], 1, coeffs, 1);
594  }
595 }
596 
597 /**
598  * Given the local cartesian coordinate \a Lcoord evaluate the
599  * value of physvals at this point by calling through to the
600  * StdExpansion method
601  */
603  const Array<OneD, const NekDouble> &Lcoord,
604  const Array<OneD, const NekDouble> &physvals)
605 {
606  // Evaluate point in local coordinates.
607  return StdPyrExp::v_PhysEvaluate(Lcoord, physvals);
608 }
609 
611  const Array<OneD, const NekDouble> &physvals)
612 {
613  Array<OneD, NekDouble> Lcoord(3);
614 
615  ASSERTL0(m_geom, "m_geom not defined");
616 
617  // TODO: check GetLocCoords()
618  m_geom->GetLocCoords(coord, Lcoord);
619 
620  return StdPyrExp::v_PhysEvaluate(Lcoord, physvals);
621 }
622 
623 //---------------------------------------
624 // Helper functions
625 //---------------------------------------
626 
628 {
629  return m_geom->GetCoordim();
630 }
631 
632 void PyrExp::v_GetTracePhysMap(const int face, Array<OneD, int> &outarray)
633 {
634  int nquad0 = m_base[0]->GetNumPoints();
635  int nquad1 = m_base[1]->GetNumPoints();
636  int nquad2 = m_base[2]->GetNumPoints();
637 
638  int nq0 = 0;
639  int nq1 = 0;
640 
641  switch (face)
642  {
643  case 0:
644  nq0 = nquad0;
645  nq1 = nquad1;
646  if (outarray.size() != nq0 * nq1)
647  {
648  outarray = Array<OneD, int>(nq0 * nq1);
649  }
650 
651  // Directions A and B positive
652  for (int i = 0; i < nquad0 * nquad1; ++i)
653  {
654  outarray[i] = i;
655  }
656 
657  break;
658  case 1:
659  nq0 = nquad0;
660  nq1 = nquad2;
661  if (outarray.size() != nq0 * nq1)
662  {
663  outarray = Array<OneD, int>(nq0 * nq1);
664  }
665 
666  // Direction A and B positive
667  for (int k = 0; k < nquad2; k++)
668  {
669  for (int i = 0; i < nquad0; ++i)
670  {
671  outarray[k * nquad0 + i] = (nquad0 * nquad1 * k) + i;
672  }
673  }
674 
675  break;
676  case 2:
677  nq0 = nquad1;
678  nq1 = nquad2;
679  if (outarray.size() != nq0 * nq1)
680  {
681  outarray = Array<OneD, int>(nq0 * nq1);
682  }
683 
684  // Directions A and B positive
685  for (int j = 0; j < nquad1 * nquad2; ++j)
686  {
687  outarray[j] = nquad0 - 1 + j * nquad0;
688  }
689  break;
690  case 3:
691 
692  nq0 = nquad0;
693  nq1 = nquad2;
694  if (outarray.size() != nq0 * nq1)
695  {
696  outarray = Array<OneD, int>(nq0 * nq1);
697  }
698 
699  // Direction A and B positive
700  for (int k = 0; k < nquad2; k++)
701  {
702  for (int i = 0; i < nquad0; ++i)
703  {
704  outarray[k * nquad0 + i] =
705  nquad0 * (nquad1 - 1) + (nquad0 * nquad1 * k) + i;
706  }
707  }
708  break;
709  case 4:
710  nq0 = nquad1;
711  nq1 = nquad2;
712 
713  if (outarray.size() != nq0 * nq1)
714  {
715  outarray = Array<OneD, int>(nq0 * nq1);
716  }
717 
718  // Directions A and B positive
719  for (int j = 0; j < nquad1 * nquad2; ++j)
720  {
721  outarray[j] = j * nquad0;
722  }
723  break;
724  default:
725  ASSERTL0(false, "face value (> 4) is out of range");
726  break;
727  }
728 }
729 
730 void PyrExp::v_ComputeTraceNormal(const int face)
731 {
732  const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
733  GetGeom()->GetMetricInfo();
734 
736  for (int i = 0; i < ptsKeys.size(); ++i)
737  {
738  // Need at least 2 points for computing normals
739  if (ptsKeys[i].GetNumPoints() == 1)
740  {
741  LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
742  ptsKeys[i] = pKey;
743  }
744  }
745 
746  SpatialDomains::GeomType type = geomFactors->GetGtype();
747  const Array<TwoD, const NekDouble> &df =
748  geomFactors->GetDerivFactors(ptsKeys);
749  const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
750 
751  LibUtilities::BasisKey tobasis0 = GetTraceBasisKey(face, 0);
752  LibUtilities::BasisKey tobasis1 = GetTraceBasisKey(face, 1);
753 
754  // Number of quadrature points in face expansion.
755  int nq_face = tobasis0.GetNumPoints() * tobasis1.GetNumPoints();
756 
757  int vCoordDim = GetCoordim();
758  int i;
759 
762  for (i = 0; i < vCoordDim; ++i)
763  {
764  normal[i] = Array<OneD, NekDouble>(nq_face);
765  }
766 
767  size_t nqb = nq_face;
768  size_t nbnd = face;
771 
772  // Regular geometry case
773  if (type == SpatialDomains::eRegular ||
775  {
776  NekDouble fac;
777  // Set up normals
778  switch (face)
779  {
780  case 0:
781  {
782  for (i = 0; i < vCoordDim; ++i)
783  {
784  normal[i][0] = -df[3 * i + 2][0];
785  }
786  break;
787  }
788  case 1:
789  {
790  for (i = 0; i < vCoordDim; ++i)
791  {
792  normal[i][0] = -df[3 * i + 1][0];
793  }
794  break;
795  }
796  case 2:
797  {
798  for (i = 0; i < vCoordDim; ++i)
799  {
800  normal[i][0] = df[3 * i][0] + df[3 * i + 2][0];
801  }
802  break;
803  }
804  case 3:
805  {
806  for (i = 0; i < vCoordDim; ++i)
807  {
808  normal[i][0] = df[3 * i + 1][0] + df[3 * i + 2][0];
809  }
810  break;
811  }
812  case 4:
813  {
814  for (i = 0; i < vCoordDim; ++i)
815  {
816  normal[i][0] = -df[3 * i][0];
817  }
818  break;
819  }
820  default:
821  ASSERTL0(false, "face is out of range (face < 4)");
822  }
823 
824  // Normalise resulting vector.
825  fac = 0.0;
826  for (i = 0; i < vCoordDim; ++i)
827  {
828  fac += normal[i][0] * normal[i][0];
829  }
830  fac = 1.0 / sqrt(fac);
831 
832  Vmath::Fill(nqb, fac, length, 1);
833 
834  for (i = 0; i < vCoordDim; ++i)
835  {
836  Vmath::Fill(nq_face, fac * normal[i][0], normal[i], 1);
837  }
838  }
839  else
840  {
841  // Set up deformed normals.
842  int j, k;
843 
844  int nq0 = ptsKeys[0].GetNumPoints();
845  int nq1 = ptsKeys[1].GetNumPoints();
846  int nq2 = ptsKeys[2].GetNumPoints();
847  int nq01 = nq0 * nq1;
848  int nqtot;
849 
850  // Determine number of quadrature points on the face.
851  if (face == 0)
852  {
853  nqtot = nq0 * nq1;
854  }
855  else if (face == 1 || face == 3)
856  {
857  nqtot = nq0 * nq2;
858  }
859  else
860  {
861  nqtot = nq1 * nq2;
862  }
863 
864  LibUtilities::PointsKey points0;
865  LibUtilities::PointsKey points1;
866 
867  Array<OneD, NekDouble> faceJac(nqtot);
868  Array<OneD, NekDouble> normals(vCoordDim * nqtot, 0.0);
869 
870  // Extract Jacobian along face and recover local derivatives
871  // (dx/dr) for polynomial interpolation by multiplying m_gmat by
872  // jacobian
873  switch (face)
874  {
875  case 0:
876  {
877  for (j = 0; j < nq01; ++j)
878  {
879  normals[j] = -df[2][j] * jac[j];
880  normals[nqtot + j] = -df[5][j] * jac[j];
881  normals[2 * nqtot + j] = -df[8][j] * jac[j];
882  faceJac[j] = jac[j];
883  }
884 
885  points0 = ptsKeys[0];
886  points1 = ptsKeys[1];
887  break;
888  }
889 
890  case 1:
891  {
892  for (j = 0; j < nq0; ++j)
893  {
894  for (k = 0; k < nq2; ++k)
895  {
896  int tmp = j + nq01 * k;
897  normals[j + k * nq0] = -df[1][tmp] * jac[tmp];
898  normals[nqtot + j + k * nq0] = -df[4][tmp] * jac[tmp];
899  normals[2 * nqtot + j + k * nq0] =
900  -df[7][tmp] * jac[tmp];
901  faceJac[j + k * nq0] = jac[tmp];
902  }
903  }
904 
905  points0 = ptsKeys[0];
906  points1 = ptsKeys[2];
907  break;
908  }
909 
910  case 2:
911  {
912  for (j = 0; j < nq1; ++j)
913  {
914  for (k = 0; k < nq2; ++k)
915  {
916  int tmp = nq0 - 1 + nq0 * j + nq01 * k;
917  normals[j + k * nq1] =
918  (df[0][tmp] + df[2][tmp]) * jac[tmp];
919  normals[nqtot + j + k * nq1] =
920  (df[3][tmp] + df[5][tmp]) * jac[tmp];
921  normals[2 * nqtot + j + k * nq1] =
922  (df[6][tmp] + df[8][tmp]) * jac[tmp];
923  faceJac[j + k * nq1] = jac[tmp];
924  }
925  }
926 
927  points0 = ptsKeys[1];
928  points1 = ptsKeys[2];
929  break;
930  }
931 
932  case 3:
933  {
934  for (j = 0; j < nq0; ++j)
935  {
936  for (k = 0; k < nq2; ++k)
937  {
938  int tmp = nq0 * (nq1 - 1) + j + nq01 * k;
939  normals[j + k * nq0] =
940  (df[1][tmp] + df[2][tmp]) * jac[tmp];
941  normals[nqtot + j + k * nq0] =
942  (df[4][tmp] + df[5][tmp]) * jac[tmp];
943  normals[2 * nqtot + j + k * nq0] =
944  (df[7][tmp] + df[8][tmp]) * jac[tmp];
945  faceJac[j + k * nq0] = jac[tmp];
946  }
947  }
948 
949  points0 = ptsKeys[0];
950  points1 = ptsKeys[2];
951  break;
952  }
953 
954  case 4:
955  {
956  for (j = 0; j < nq1; ++j)
957  {
958  for (k = 0; k < nq2; ++k)
959  {
960  int tmp = j * nq0 + nq01 * k;
961  normals[j + k * nq1] = -df[0][tmp] * jac[tmp];
962  normals[nqtot + j + k * nq1] = -df[3][tmp] * jac[tmp];
963  normals[2 * nqtot + j + k * nq1] =
964  -df[6][tmp] * jac[tmp];
965  faceJac[j + k * nq1] = jac[tmp];
966  }
967  }
968 
969  points0 = ptsKeys[1];
970  points1 = ptsKeys[2];
971  break;
972  }
973 
974  default:
975  ASSERTL0(false, "face is out of range (face < 4)");
976  }
977 
978  Array<OneD, NekDouble> work(nq_face, 0.0);
979  // Interpolate Jacobian and invert
980  LibUtilities::Interp2D(points0, points1, faceJac,
981  tobasis0.GetPointsKey(), tobasis1.GetPointsKey(),
982  work);
983  Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
984 
985  // Interpolate normal and multiply by inverse Jacobian.
986  for (i = 0; i < vCoordDim; ++i)
987  {
988  LibUtilities::Interp2D(points0, points1, &normals[i * nqtot],
989  tobasis0.GetPointsKey(),
990  tobasis1.GetPointsKey(), &normal[i][0]);
991  Vmath::Vmul(nq_face, work, 1, normal[i], 1, normal[i], 1);
992  }
993 
994  // Normalise to obtain unit normals.
995  Vmath::Zero(nq_face, work, 1);
996  for (i = 0; i < GetCoordim(); ++i)
997  {
998  Vmath::Vvtvp(nq_face, normal[i], 1, normal[i], 1, work, 1, work, 1);
999  }
1000 
1001  Vmath::Vsqrt(nq_face, work, 1, work, 1);
1002  Vmath::Sdiv(nq_face, 1.0, work, 1, work, 1);
1003 
1004  Vmath::Vcopy(nqb, work, 1, length, 1);
1005 
1006  for (i = 0; i < GetCoordim(); ++i)
1007  {
1008  Vmath::Vmul(nq_face, normal[i], 1, work, 1, normal[i], 1);
1009  }
1010  }
1011 }
1012 
1014  const StdRegions::StdMatrixKey &mkey)
1015 {
1016  int nq = GetTotPoints();
1017 
1018  // Calculate sqrt of the Jacobian
1020  Array<OneD, NekDouble> sqrt_jac(nq);
1021  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1022  {
1023  Vmath::Vsqrt(nq, jac, 1, sqrt_jac, 1);
1024  }
1025  else
1026  {
1027  Vmath::Fill(nq, sqrt(jac[0]), sqrt_jac, 1);
1028  }
1029 
1030  // Multiply array by sqrt(Jac)
1031  Vmath::Vmul(nq, sqrt_jac, 1, array, 1, array, 1);
1032 
1033  // Apply std region filter
1034  StdPyrExp::v_SVVLaplacianFilter(array, mkey);
1035 
1036  // Divide by sqrt(Jac)
1037  Vmath::Vdiv(nq, array, 1, sqrt_jac, 1, array, 1);
1038 }
1039 
1040 //---------------------------------------
1041 // Matrix creation functions
1042 //---------------------------------------
1043 
1045 {
1046  DNekMatSharedPtr returnval;
1047 
1048  switch (mkey.GetMatrixType())
1049  {
1056  returnval = Expansion3D::v_GenMatrix(mkey);
1057  break;
1058  default:
1059  returnval = StdPyrExp::v_GenMatrix(mkey);
1060  }
1061 
1062  return returnval;
1063 }
1064 
1066 {
1067  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1068  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1069  LibUtilities::BasisKey bkey2 = m_base[2]->GetBasisKey();
1071  MemoryManager<StdPyrExp>::AllocateSharedPtr(bkey0, bkey1, bkey2);
1072 
1073  return tmp->GetStdMatrix(mkey);
1074 }
1075 
1077 {
1078  return m_matrixManager[mkey];
1079 }
1080 
1082 {
1083  return m_staticCondMatrixManager[mkey];
1084 }
1085 
1087 {
1088  m_staticCondMatrixManager.DeleteObject(mkey);
1089 }
1090 
1092 {
1093  if (m_metrics.count(eMetricQuadrature) == 0)
1094  {
1096  }
1097 
1098  int i, j;
1099  const unsigned int nqtot = GetTotPoints();
1100  const unsigned int dim = 3;
1101  const MetricType m[3][3] = {
1105 
1106  for (unsigned int i = 0; i < dim; ++i)
1107  {
1108  for (unsigned int j = i; j < dim; ++j)
1109  {
1110  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1111  }
1112  }
1113 
1114  // Define shorthand synonyms for m_metrics storage
1115  Array<OneD, NekDouble> g0(m_metrics[m[0][0]]);
1116  Array<OneD, NekDouble> g1(m_metrics[m[1][1]]);
1117  Array<OneD, NekDouble> g2(m_metrics[m[2][2]]);
1118  Array<OneD, NekDouble> g3(m_metrics[m[0][1]]);
1119  Array<OneD, NekDouble> g4(m_metrics[m[0][2]]);
1120  Array<OneD, NekDouble> g5(m_metrics[m[1][2]]);
1121 
1122  // Allocate temporary storage
1123  Array<OneD, NekDouble> alloc(9 * nqtot, 0.0);
1124  Array<OneD, NekDouble> h0(nqtot, alloc);
1125  Array<OneD, NekDouble> h1(nqtot, alloc + 1 * nqtot);
1126  Array<OneD, NekDouble> h2(nqtot, alloc + 2 * nqtot);
1127  Array<OneD, NekDouble> wsp1(nqtot, alloc + 3 * nqtot);
1128  Array<OneD, NekDouble> wsp2(nqtot, alloc + 4 * nqtot);
1129  Array<OneD, NekDouble> wsp3(nqtot, alloc + 5 * nqtot);
1130  Array<OneD, NekDouble> wsp4(nqtot, alloc + 6 * nqtot);
1131  Array<OneD, NekDouble> wsp5(nqtot, alloc + 7 * nqtot);
1132  Array<OneD, NekDouble> wsp6(nqtot, alloc + 8 * nqtot);
1133 
1134  const Array<TwoD, const NekDouble> &df =
1135  m_metricinfo->GetDerivFactors(GetPointsKeys());
1136  const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
1137  const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
1138  const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
1139  const unsigned int nquad0 = m_base[0]->GetNumPoints();
1140  const unsigned int nquad1 = m_base[1]->GetNumPoints();
1141  const unsigned int nquad2 = m_base[2]->GetNumPoints();
1142 
1143  // Populate collapsed coordinate arrays h0, h1 and h2.
1144  for (j = 0; j < nquad2; ++j)
1145  {
1146  for (i = 0; i < nquad1; ++i)
1147  {
1148  Vmath::Fill(nquad0, 2.0 / (1.0 - z2[j]),
1149  &h0[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1150  Vmath::Fill(nquad0, 1.0 / (1.0 - z2[j]),
1151  &h1[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1152  Vmath::Fill(nquad0, (1.0 + z1[i]) / (1.0 - z2[j]),
1153  &h2[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1154  }
1155  }
1156  for (i = 0; i < nquad0; i++)
1157  {
1158  Blas::Dscal(nquad1 * nquad2, 1 + z0[i], &h1[0] + i, nquad0);
1159  }
1160 
1161  // Step 3. Construct combined metric terms for physical space to
1162  // collapsed coordinate system.
1163  // Order of construction optimised to minimise temporary storage
1164  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1165  {
1166  // f_{1k}
1167  Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &df[2][0], 1, &h1[0], 1,
1168  &wsp1[0], 1);
1169  Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &df[5][0], 1, &h1[0], 1,
1170  &wsp2[0], 1);
1171  Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &df[8][0], 1, &h1[0], 1,
1172  &wsp3[0], 1);
1173 
1174  // g0
1175  Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp1[0], 1, &wsp2[0], 1, &wsp2[0],
1176  1, &g0[0], 1);
1177  Vmath::Vvtvp(nqtot, &wsp3[0], 1, &wsp3[0], 1, &g0[0], 1, &g0[0], 1);
1178 
1179  // g4
1180  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp1[0], 1, &df[5][0], 1, &wsp2[0],
1181  1, &g4[0], 1);
1182  Vmath::Vvtvp(nqtot, &df[8][0], 1, &wsp3[0], 1, &g4[0], 1, &g4[0], 1);
1183 
1184  // f_{2k}
1185  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &h0[0], 1, &df[2][0], 1, &h2[0], 1,
1186  &wsp4[0], 1);
1187  Vmath::Vvtvvtp(nqtot, &df[4][0], 1, &h0[0], 1, &df[5][0], 1, &h2[0], 1,
1188  &wsp5[0], 1);
1189  Vmath::Vvtvvtp(nqtot, &df[7][0], 1, &h0[0], 1, &df[8][0], 1, &h2[0], 1,
1190  &wsp6[0], 1);
1191 
1192  // g1
1193  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1194  1, &g1[0], 1);
1195  Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g1[0], 1, &g1[0], 1);
1196 
1197  // g3
1198  Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp4[0], 1, &wsp2[0], 1, &wsp5[0],
1199  1, &g3[0], 1);
1200  Vmath::Vvtvp(nqtot, &wsp3[0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1201 
1202  // g5
1203  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0],
1204  1, &g5[0], 1);
1205  Vmath::Vvtvp(nqtot, &df[8][0], 1, &wsp6[0], 1, &g5[0], 1, &g5[0], 1);
1206 
1207  // g2
1208  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &df[2][0], 1, &df[5][0], 1,
1209  &df[5][0], 1, &g2[0], 1);
1210  Vmath::Vvtvp(nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1211  }
1212  else
1213  {
1214  // f_{1k}
1215  Vmath::Svtsvtp(nqtot, df[0][0], &h0[0], 1, df[2][0], &h1[0], 1,
1216  &wsp1[0], 1);
1217  Vmath::Svtsvtp(nqtot, df[3][0], &h0[0], 1, df[5][0], &h1[0], 1,
1218  &wsp2[0], 1);
1219  Vmath::Svtsvtp(nqtot, df[6][0], &h0[0], 1, df[8][0], &h1[0], 1,
1220  &wsp3[0], 1);
1221 
1222  // g0
1223  Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp1[0], 1, &wsp2[0], 1, &wsp2[0],
1224  1, &g0[0], 1);
1225  Vmath::Vvtvp(nqtot, &wsp3[0], 1, &wsp3[0], 1, &g0[0], 1, &g0[0], 1);
1226 
1227  // g4
1228  Vmath::Svtsvtp(nqtot, df[2][0], &wsp1[0], 1, df[5][0], &wsp2[0], 1,
1229  &g4[0], 1);
1230  Vmath::Svtvp(nqtot, df[8][0], &wsp3[0], 1, &g4[0], 1, &g4[0], 1);
1231 
1232  // f_{2k}
1233  Vmath::Svtsvtp(nqtot, df[1][0], &h0[0], 1, df[2][0], &h2[0], 1,
1234  &wsp4[0], 1);
1235  Vmath::Svtsvtp(nqtot, df[4][0], &h0[0], 1, df[5][0], &h2[0], 1,
1236  &wsp5[0], 1);
1237  Vmath::Svtsvtp(nqtot, df[7][0], &h0[0], 1, df[8][0], &h2[0], 1,
1238  &wsp6[0], 1);
1239 
1240  // g1
1241  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1242  1, &g1[0], 1);
1243  Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g1[0], 1, &g1[0], 1);
1244 
1245  // g3
1246  Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp4[0], 1, &wsp2[0], 1, &wsp5[0],
1247  1, &g3[0], 1);
1248  Vmath::Vvtvp(nqtot, &wsp3[0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1249 
1250  // g5
1251  Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1,
1252  &g5[0], 1);
1253  Vmath::Svtvp(nqtot, df[8][0], &wsp6[0], 1, &g5[0], 1, &g5[0], 1);
1254 
1255  // g2
1256  Vmath::Fill(nqtot,
1257  df[2][0] * df[2][0] + df[5][0] * df[5][0] +
1258  df[8][0] * df[8][0],
1259  &g2[0], 1);
1260  }
1261 
1262  for (unsigned int i = 0; i < dim; ++i)
1263  {
1264  for (unsigned int j = i; j < dim; ++j)
1265  {
1266  MultiplyByQuadratureMetric(m_metrics[m[i][j]], m_metrics[m[i][j]]);
1267  }
1268  }
1269 }
1270 
1272  const Array<OneD, const NekDouble> &inarray,
1274 {
1275  // This implementation is only valid when there are no coefficients
1276  // associated to the Laplacian operator
1277  if (m_metrics.count(eMetricLaplacian00) == 0)
1278  {
1280  }
1281 
1282  int nquad0 = m_base[0]->GetNumPoints();
1283  int nquad1 = m_base[1]->GetNumPoints();
1284  int nq2 = m_base[2]->GetNumPoints();
1285  int nqtot = nquad0 * nquad1 * nq2;
1286 
1287  ASSERTL1(wsp.size() >= 6 * nqtot, "Insufficient workspace size.");
1288  ASSERTL1(m_ncoeffs <= nqtot, "Workspace not set up for ncoeffs > nqtot");
1289 
1290  const Array<OneD, const NekDouble> &base0 = m_base[0]->GetBdata();
1291  const Array<OneD, const NekDouble> &base1 = m_base[1]->GetBdata();
1292  const Array<OneD, const NekDouble> &base2 = m_base[2]->GetBdata();
1293  const Array<OneD, const NekDouble> &dbase0 = m_base[0]->GetDbdata();
1294  const Array<OneD, const NekDouble> &dbase1 = m_base[1]->GetDbdata();
1295  const Array<OneD, const NekDouble> &dbase2 = m_base[2]->GetDbdata();
1296  const Array<OneD, const NekDouble> &metric00 =
1297  m_metrics[eMetricLaplacian00];
1298  const Array<OneD, const NekDouble> &metric01 =
1299  m_metrics[eMetricLaplacian01];
1300  const Array<OneD, const NekDouble> &metric02 =
1301  m_metrics[eMetricLaplacian02];
1302  const Array<OneD, const NekDouble> &metric11 =
1303  m_metrics[eMetricLaplacian11];
1304  const Array<OneD, const NekDouble> &metric12 =
1305  m_metrics[eMetricLaplacian12];
1306  const Array<OneD, const NekDouble> &metric22 =
1307  m_metrics[eMetricLaplacian22];
1308 
1309  // Allocate temporary storage
1310  Array<OneD, NekDouble> wsp0(2 * nqtot, wsp);
1311  Array<OneD, NekDouble> wsp1(nqtot, wsp + 1 * nqtot);
1312  Array<OneD, NekDouble> wsp2(nqtot, wsp + 2 * nqtot);
1313  Array<OneD, NekDouble> wsp3(nqtot, wsp + 3 * nqtot);
1314  Array<OneD, NekDouble> wsp4(nqtot, wsp + 4 * nqtot);
1315  Array<OneD, NekDouble> wsp5(nqtot, wsp + 5 * nqtot);
1316 
1317  // LAPLACIAN MATRIX OPERATION
1318  // wsp1 = du_dxi1 = D_xi1 * inarray = D_xi1 * u
1319  // wsp2 = du_dxi2 = D_xi2 * inarray = D_xi2 * u
1320  // wsp2 = du_dxi3 = D_xi3 * inarray = D_xi3 * u
1321  StdExpansion3D::PhysTensorDeriv(inarray, wsp0, wsp1, wsp2);
1322 
1323  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1324  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1325  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1326  // especially for this purpose
1327  Vmath::Vvtvvtp(nqtot, &metric00[0], 1, &wsp0[0], 1, &metric01[0], 1,
1328  &wsp1[0], 1, &wsp3[0], 1);
1329  Vmath::Vvtvp(nqtot, &metric02[0], 1, &wsp2[0], 1, &wsp3[0], 1, &wsp3[0], 1);
1330  Vmath::Vvtvvtp(nqtot, &metric01[0], 1, &wsp0[0], 1, &metric11[0], 1,
1331  &wsp1[0], 1, &wsp4[0], 1);
1332  Vmath::Vvtvp(nqtot, &metric12[0], 1, &wsp2[0], 1, &wsp4[0], 1, &wsp4[0], 1);
1333  Vmath::Vvtvvtp(nqtot, &metric02[0], 1, &wsp0[0], 1, &metric12[0], 1,
1334  &wsp1[0], 1, &wsp5[0], 1);
1335  Vmath::Vvtvp(nqtot, &metric22[0], 1, &wsp2[0], 1, &wsp5[0], 1, &wsp5[0], 1);
1336 
1337  // outarray = m = (D_xi1 * B)^T * k
1338  // wsp1 = n = (D_xi2 * B)^T * l
1339  IProductWRTBase_SumFacKernel(dbase0, base1, base2, wsp3, outarray, wsp0,
1340  false, true, true);
1341  IProductWRTBase_SumFacKernel(base0, dbase1, base2, wsp4, wsp2, wsp0, true,
1342  false, true);
1343  Vmath::Vadd(m_ncoeffs, wsp2.get(), 1, outarray.get(), 1, outarray.get(), 1);
1344  IProductWRTBase_SumFacKernel(base0, base1, dbase2, wsp5, wsp2, wsp0, true,
1345  true, false);
1346  Vmath::Vadd(m_ncoeffs, wsp2.get(), 1, outarray.get(), 1, outarray.get(), 1);
1347 }
1348 
1349 /** @brief: This method gets all of the factors which are
1350  required as part of the Gradient Jump Penalty
1351  stabilisation and involves the product of the normal and
1352  geometric factors along the element trace.
1353 */
1354 void PyrExp::v_NormalTraceDerivFactors(
1355  Array<OneD, Array<OneD, NekDouble>> &d0factors,
1356  Array<OneD, Array<OneD, NekDouble>> &d1factors,
1357  Array<OneD, Array<OneD, NekDouble>> &d2factors)
1358 {
1359  int nquad0 = GetNumPoints(0);
1360  int nquad1 = GetNumPoints(1);
1361  int nquad2 = GetNumPoints(2);
1362 
1363  const Array<TwoD, const NekDouble> &df =
1364  m_metricinfo->GetDerivFactors(GetPointsKeys());
1365 
1366  if (d0factors.size() != 5)
1367  {
1368  d0factors = Array<OneD, Array<OneD, NekDouble>>(5);
1369  d1factors = Array<OneD, Array<OneD, NekDouble>>(5);
1370  d2factors = Array<OneD, Array<OneD, NekDouble>>(5);
1371  }
1372 
1373  if (d0factors[0].size() != nquad0 * nquad1)
1374  {
1375  d0factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1376  d1factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1377  d2factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1378  }
1379 
1380  if (d0factors[1].size() != nquad0 * nquad2)
1381  {
1382  d0factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1383  d0factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1384  d1factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1385  d1factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1386  d2factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1387  d2factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1388  }
1389 
1390  if (d0factors[2].size() != nquad1 * nquad2)
1391  {
1392  d0factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1393  d0factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1394  d1factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1395  d1factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1396  d2factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1397  d2factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1398  }
1399 
1400  // Outwards normals
1401  const Array<OneD, const Array<OneD, NekDouble>> &normal_0 =
1402  GetTraceNormal(0);
1403  const Array<OneD, const Array<OneD, NekDouble>> &normal_1 =
1404  GetTraceNormal(1);
1405  const Array<OneD, const Array<OneD, NekDouble>> &normal_2 =
1406  GetTraceNormal(2);
1407  const Array<OneD, const Array<OneD, NekDouble>> &normal_3 =
1408  GetTraceNormal(3);
1409  const Array<OneD, const Array<OneD, NekDouble>> &normal_4 =
1410  GetTraceNormal(4);
1411 
1412  int ncoords = normal_0.size();
1413 
1414  // first gather together standard cartesian inner products
1415  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1416  {
1417  // face 0
1418  for (int i = 0; i < nquad0 * nquad1; ++i)
1419  {
1420  d0factors[0][i] = df[0][i] * normal_0[0][i];
1421  d1factors[0][i] = df[1][i] * normal_0[0][i];
1422  d2factors[0][i] = df[2][i] * normal_0[0][i];
1423  }
1424 
1425  for (int n = 1; n < ncoords; ++n)
1426  {
1427  for (int i = 0; i < nquad0 * nquad1; ++i)
1428  {
1429  d0factors[0][i] += df[3 * n][i] * normal_0[n][i];
1430  d1factors[0][i] += df[3 * n + 1][i] * normal_0[n][i];
1431  d2factors[0][i] += df[3 * n + 2][i] * normal_0[n][i];
1432  }
1433  }
1434 
1435  // faces 1 and 3
1436  for (int j = 0; j < nquad2; ++j)
1437  {
1438  for (int i = 0; i < nquad0; ++i)
1439  {
1440  d0factors[1][i] = df[0][j * nquad0 * nquad1 + i] *
1441  normal_1[0][j * nquad0 + i];
1442  d0factors[3][i] =
1443  df[0][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1444  normal_3[0][j * nquad0 + i];
1445  d1factors[1][i] = df[1][j * nquad0 * nquad1 + i] *
1446  normal_1[0][j * nquad0 + i];
1447  d1factors[3][i] =
1448  df[1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1449  normal_3[0][j * nquad0 + i];
1450  d2factors[1][i] = df[2][j * nquad0 * nquad1 + i] *
1451  normal_1[0][j * nquad0 + i];
1452  d2factors[3][i] =
1453  df[2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1454  normal_3[0][j * nquad0 + i];
1455  }
1456  }
1457 
1458  for (int n = 1; n < ncoords; ++n)
1459  {
1460  for (int j = 0; j < nquad2; ++j)
1461  {
1462  for (int i = 0; i < nquad0; ++i)
1463  {
1464  d0factors[1][i] = df[3 * n][j * nquad0 * nquad1 + i] *
1465  normal_1[0][j * nquad0 + i];
1466  d0factors[3][i] =
1467  df[3 * n][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1468  normal_3[0][j * nquad0 + i];
1469  d1factors[1][i] = df[3 * n + 1][j * nquad0 * nquad1 + i] *
1470  normal_1[0][j * nquad0 + i];
1471  d1factors[3][i] =
1472  df[3 * n + 1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1473  normal_3[0][j * nquad0 + i];
1474  d2factors[1][i] = df[3 * n + 2][j * nquad0 * nquad1 + i] *
1475  normal_1[0][j * nquad0 + i];
1476  d2factors[3][i] =
1477  df[3 * n + 2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1478  normal_3[0][j * nquad0 + i];
1479  }
1480  }
1481  }
1482 
1483  // faces 2 and 4
1484  for (int j = 0; j < nquad2; ++j)
1485  {
1486  for (int i = 0; i < nquad1; ++i)
1487  {
1488  d0factors[2][j * nquad1 + i] =
1489  df[0][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1490  normal_2[0][j * nquad1 + i];
1491  d0factors[4][j * nquad1 + i] =
1492  df[0][j * nquad0 * nquad1 + i * nquad0] *
1493  normal_4[0][j * nquad1 + i];
1494  d1factors[2][j * nquad1 + i] =
1495  df[1][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1496  normal_2[0][j * nquad1 + i];
1497  d1factors[4][j * nquad1 + i] =
1498  df[1][j * nquad0 * nquad1 + i * nquad0] *
1499  normal_4[0][j * nquad1 + i];
1500  d2factors[2][j * nquad1 + i] =
1501  df[2][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1502  normal_2[0][j * nquad1 + i];
1503  d2factors[4][j * nquad1 + i] =
1504  df[2][j * nquad0 * nquad1 + i * nquad0] *
1505  normal_4[0][j * nquad1 + i];
1506  }
1507  }
1508 
1509  for (int n = 1; n < ncoords; ++n)
1510  {
1511  for (int j = 0; j < nquad2; ++j)
1512  {
1513  for (int i = 0; i < nquad1; ++i)
1514  {
1515  d0factors[2][j * nquad1 + i] +=
1516  df[3 * n][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1517  normal_2[n][j * nquad0 + i];
1518  d0factors[4][j * nquad0 + i] +=
1519  df[3 * n][i * nquad0 + j * nquad0 * nquad1] *
1520  normal_4[n][j * nquad0 + i];
1521  d1factors[2][j * nquad1 + i] +=
1522  df[3 * n + 1]
1523  [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1524  normal_2[n][j * nquad0 + i];
1525  d1factors[4][j * nquad0 + i] +=
1526  df[3 * n + 1][i * nquad0 + j * nquad0 * nquad1] *
1527  normal_4[n][j * nquad0 + i];
1528  d2factors[2][j * nquad1 + i] +=
1529  df[3 * n + 2]
1530  [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1531  normal_2[n][j * nquad0 + i];
1532  d2factors[4][j * nquad0 + i] +=
1533  df[3 * n + 2][i * nquad0 + j * nquad0 * nquad1] *
1534  normal_4[n][j * nquad0 + i];
1535  }
1536  }
1537  }
1538  }
1539  else
1540  {
1541  // Face 0
1542  for (int i = 0; i < nquad0 * nquad1; ++i)
1543  {
1544  d0factors[0][i] = df[0][0] * normal_0[0][i];
1545  d1factors[0][i] = df[1][0] * normal_0[0][i];
1546  d2factors[0][i] = df[2][0] * normal_0[0][i];
1547  }
1548 
1549  for (int n = 1; n < ncoords; ++n)
1550  {
1551  for (int i = 0; i < nquad0 * nquad1; ++i)
1552  {
1553  d0factors[0][i] += df[3 * n][0] * normal_0[n][i];
1554  d1factors[0][i] += df[3 * n + 1][0] * normal_0[n][i];
1555  d2factors[0][i] += df[3 * n + 2][0] * normal_0[n][i];
1556  }
1557  }
1558 
1559  // faces 1 and 3
1560  for (int i = 0; i < nquad0 * nquad2; ++i)
1561  {
1562  d0factors[1][i] = df[0][0] * normal_1[0][i];
1563  d0factors[3][i] = df[0][0] * normal_3[0][i];
1564 
1565  d1factors[1][i] = df[1][0] * normal_1[0][i];
1566  d1factors[3][i] = df[1][0] * normal_3[0][i];
1567 
1568  d2factors[1][i] = df[2][0] * normal_1[0][i];
1569  d2factors[3][i] = df[2][0] * normal_3[0][i];
1570  }
1571 
1572  for (int n = 1; n < ncoords; ++n)
1573  {
1574  for (int i = 0; i < nquad0 * nquad2; ++i)
1575  {
1576  d0factors[1][i] += df[3 * n][0] * normal_1[n][i];
1577  d0factors[3][i] += df[3 * n][0] * normal_3[n][i];
1578 
1579  d1factors[1][i] += df[3 * n + 1][0] * normal_1[n][i];
1580  d1factors[3][i] += df[3 * n + 1][0] * normal_3[n][i];
1581 
1582  d2factors[1][i] += df[3 * n + 2][0] * normal_1[n][i];
1583  d2factors[3][i] += df[3 * n + 2][0] * normal_3[n][i];
1584  }
1585  }
1586 
1587  // faces 2 and 4
1588  for (int i = 0; i < nquad1 * nquad2; ++i)
1589  {
1590  d0factors[2][i] = df[0][0] * normal_2[0][i];
1591  d0factors[4][i] = df[0][0] * normal_4[0][i];
1592 
1593  d1factors[2][i] = df[1][0] * normal_2[0][i];
1594  d1factors[4][i] = df[1][0] * normal_4[0][i];
1595 
1596  d2factors[2][i] = df[2][0] * normal_2[0][i];
1597  d2factors[4][i] = df[2][0] * normal_4[0][i];
1598  }
1599 
1600  for (int n = 1; n < ncoords; ++n)
1601  {
1602  for (int i = 0; i < nquad1 * nquad2; ++i)
1603  {
1604  d0factors[2][i] += df[3 * n][0] * normal_2[n][i];
1605  d0factors[4][i] += df[3 * n][0] * normal_4[n][i];
1606 
1607  d1factors[2][i] += df[3 * n + 1][0] * normal_2[n][i];
1608  d1factors[4][i] += df[3 * n + 1][0] * normal_4[n][i];
1609 
1610  d2factors[2][i] += df[3 * n + 2][0] * normal_2[n][i];
1611  d2factors[4][i] += df[3 * n + 2][0] * normal_4[n][i];
1612  }
1613  }
1614  }
1615 }
1616 
1617 } // namespace LocalRegions
1618 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
Describes the specification for a Basis.
Definition: Basis.h:50
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:130
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:147
Defines a specification for a set of points.
Definition: Points.h:59
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
std::map< int, NormalVector > m_traceNormals
Definition: Expansion.h:275
std::map< int, Array< OneD, NekDouble > > m_elmtBndNormDirElmtLen
the element length in each element boundary(Vertex, edge or face) normal direction calculated based o...
Definition: Expansion.h:285
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:166
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:272
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:273
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: Expansion.cpp:524
virtual int v_GetCoordim()
Definition: PyrExp.cpp:627
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
Definition: PyrExp.cpp:506
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: PyrExp.cpp:1044
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: PyrExp.h:169
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
Definition: PyrExp.cpp:1013
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: PyrExp.cpp:1076
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: PyrExp.cpp:1086
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:128
virtual void v_GetTracePhysMap(const int face, Array< OneD, int > &outarray)
Definition: PyrExp.cpp:632
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: PyrExp.cpp:1065
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:282
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: PyrExp.cpp:355
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:610
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: PyrExp.h:171
NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
Definition: PyrExp.cpp:602
void v_ComputeTraceNormal(const int face)
Definition: PyrExp.cpp:730
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: PyrExp.cpp:540
virtual void v_ComputeLaplacianMetric()
Definition: PyrExp.cpp:1091
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
Definition: PyrExp.cpp:499
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:45
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: PyrExp.cpp:1081
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
Definition: PyrExp.cpp:523
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:100
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:227
virtual void v_ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int mode_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType)
Definition: PyrExp.cpp:547
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Calculates the inner product .
Definition: PyrExp.cpp:348
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: PyrExp.cpp:288
virtual void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
Definition: PyrExp.cpp:1271
virtual void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: PyrExp.cpp:400
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
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)
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:130
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:140
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:163
const LibUtilities::PointsKeyVector GetPointsKeys() const
const LibUtilities::BasisKey GetTraceBasisKey(const int i, int k=-1) const
This function returns the basis key belonging to the i-th trace.
Definition: StdExpansion.h:307
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:375
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space.
Definition: StdExpansion.h:432
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:213
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space.
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:226
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:731
Array< OneD, LibUtilities::BasisSharedPtr > m_base
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:85
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:168
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
Definition: Blas.hpp:154
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:106
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:250
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
GeomType
Indicates the type of element geometry.
@ eRegular
Geometry is straight-sided with constant geometric factors.
@ eMovingRegular
Currently unused.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< PyrGeom > PyrGeomSharedPtr
Definition: PyrGeom.h:77
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::shared_ptr< StdPyrExp > StdPyrExpSharedPtr
Definition: StdPyrExp.h:248
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:534
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:751
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:209
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:622
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:574
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:359
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
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:324
void Vdiv(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:284
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
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:692
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291