Nektar++
PrismExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File PrismExp.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: PrismExp routines
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
39 #include <LocalRegions/PrismExp.h>
40 #include <SpatialDomains/SegGeom.h>
41 
42 using namespace std;
43 
44 namespace Nektar
45 {
46 namespace LocalRegions
47 {
48 
49 PrismExp::PrismExp(const LibUtilities::BasisKey &Ba,
50  const LibUtilities::BasisKey &Bb,
51  const LibUtilities::BasisKey &Bc,
53  : StdExpansion(LibUtilities::StdPrismData::getNumberOfCoefficients(
54  Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
55  3, Ba, Bb, Bc),
56  StdExpansion3D(LibUtilities::StdPrismData::getNumberOfCoefficients(
57  Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
58  Ba, Bb, Bc),
59  StdPrismExp(Ba, Bb, Bc), Expansion(geom), Expansion3D(geom),
60  m_matrixManager(
61  std::bind(&Expansion3D::CreateMatrix, this, std::placeholders::_1),
62  std::string("PrismExpMatrix")),
63  m_staticCondMatrixManager(std::bind(&Expansion::CreateStaticCondMatrix,
64  this, std::placeholders::_1),
65  std::string("PrismExpStaticCondMatrix"))
66 {
67 }
68 
70  : StdExpansion(T), StdExpansion3D(T), StdPrismExp(T), Expansion(T),
71  Expansion3D(T), m_matrixManager(T.m_matrixManager),
72  m_staticCondMatrixManager(T.m_staticCondMatrixManager)
73 {
74 }
75 
77 {
78 }
79 
80 //-------------------------------
81 // Integration Methods
82 //-------------------------------
83 
84 /**
85  * \brief Integrate the physical point list \a inarray over prismatic
86  * region and return the value.
87  *
88  * Inputs:\n
89  *
90  * - \a inarray: definition of function to be returned at quadrature
91  * point of expansion.
92  *
93  * Outputs:\n
94  *
95  * - returns \f$\int^1_{-1}\int^1_{-1}\int^1_{-1} u(\bar \eta_1,
96  * \xi_2, \xi_3) J[i,j,k] d \bar \eta_1 d \xi_2 d \xi_3 \f$ \n \f$ =
97  * \sum_{i=0}^{Q_1 - 1} \sum_{j=0}^{Q_2 - 1} \sum_{k=0}^{Q_3 - 1}
98  * u(\bar \eta_{1i}^{0,0}, \xi_{2j}^{0,0},\xi_{3k}^{1,0})w_{i}^{0,0}
99  * w_{j}^{0,0} \hat w_{k}^{1,0} \f$ \n where \f$ inarray[i,j, k] =
100  * u(\bar \eta_{1i}^{0,0}, \xi_{2j}^{0,0},\xi_{3k}^{1,0}) \f$, \n
101  * \f$\hat w_{i}^{1,0} = \frac {w_{j}^{1,0}} {2} \f$ \n and \f$
102  * J[i,j,k] \f$ is the Jacobian evaluated at the quadrature point.
103  */
105 {
106  int nquad0 = m_base[0]->GetNumPoints();
107  int nquad1 = m_base[1]->GetNumPoints();
108  int nquad2 = m_base[2]->GetNumPoints();
110  Array<OneD, NekDouble> tmp(nquad0 * nquad1 * nquad2);
111 
112  // Multiply inarray with Jacobian
113  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
114  {
115  Vmath::Vmul(nquad0 * nquad1 * nquad2, &jac[0], 1,
116  (NekDouble *)&inarray[0], 1, &tmp[0], 1);
117  }
118  else
119  {
120  Vmath::Smul(nquad0 * nquad1 * nquad2, (NekDouble)jac[0],
121  (NekDouble *)&inarray[0], 1, &tmp[0], 1);
122  }
123 
124  // Call StdPrismExp version.
125  return StdPrismExp::v_Integral(tmp);
126 }
127 
128 //----------------------------
129 // Differentiation Methods
130 //----------------------------
132  Array<OneD, NekDouble> &out_d0,
133  Array<OneD, NekDouble> &out_d1,
134  Array<OneD, NekDouble> &out_d2)
135 {
136  int nqtot = GetTotPoints();
137 
139  m_metricinfo->GetDerivFactors(GetPointsKeys());
140  Array<OneD, NekDouble> diff0(nqtot);
141  Array<OneD, NekDouble> diff1(nqtot);
142  Array<OneD, NekDouble> diff2(nqtot);
143 
144  StdPrismExp::v_PhysDeriv(inarray, diff0, diff1, diff2);
145 
146  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
147  {
148  if (out_d0.size())
149  {
150  Vmath::Vmul(nqtot, &df[0][0], 1, &diff0[0], 1, &out_d0[0], 1);
151  Vmath::Vvtvp(nqtot, &df[1][0], 1, &diff1[0], 1, &out_d0[0], 1,
152  &out_d0[0], 1);
153  Vmath::Vvtvp(nqtot, &df[2][0], 1, &diff2[0], 1, &out_d0[0], 1,
154  &out_d0[0], 1);
155  }
156 
157  if (out_d1.size())
158  {
159  Vmath::Vmul(nqtot, &df[3][0], 1, &diff0[0], 1, &out_d1[0], 1);
160  Vmath::Vvtvp(nqtot, &df[4][0], 1, &diff1[0], 1, &out_d1[0], 1,
161  &out_d1[0], 1);
162  Vmath::Vvtvp(nqtot, &df[5][0], 1, &diff2[0], 1, &out_d1[0], 1,
163  &out_d1[0], 1);
164  }
165 
166  if (out_d2.size())
167  {
168  Vmath::Vmul(nqtot, &df[6][0], 1, &diff0[0], 1, &out_d2[0], 1);
169  Vmath::Vvtvp(nqtot, &df[7][0], 1, &diff1[0], 1, &out_d2[0], 1,
170  &out_d2[0], 1);
171  Vmath::Vvtvp(nqtot, &df[8][0], 1, &diff2[0], 1, &out_d2[0], 1,
172  &out_d2[0], 1);
173  }
174  }
175  else // regular geometry
176  {
177  if (out_d0.size())
178  {
179  Vmath::Smul(nqtot, df[0][0], &diff0[0], 1, &out_d0[0], 1);
180  Blas::Daxpy(nqtot, df[1][0], &diff1[0], 1, &out_d0[0], 1);
181  Blas::Daxpy(nqtot, df[2][0], &diff2[0], 1, &out_d0[0], 1);
182  }
183 
184  if (out_d1.size())
185  {
186  Vmath::Smul(nqtot, df[3][0], &diff0[0], 1, &out_d1[0], 1);
187  Blas::Daxpy(nqtot, df[4][0], &diff1[0], 1, &out_d1[0], 1);
188  Blas::Daxpy(nqtot, df[5][0], &diff2[0], 1, &out_d1[0], 1);
189  }
190 
191  if (out_d2.size())
192  {
193  Vmath::Smul(nqtot, df[6][0], &diff0[0], 1, &out_d2[0], 1);
194  Blas::Daxpy(nqtot, df[7][0], &diff1[0], 1, &out_d2[0], 1);
195  Blas::Daxpy(nqtot, df[8][0], &diff2[0], 1, &out_d2[0], 1);
196  }
197  }
198 }
199 
200 //---------------------------------------
201 // Transforms
202 //---------------------------------------
203 
204 /**
205  * \brief Forward transform from physical quadrature space stored in
206  * \a inarray and evaluate the expansion coefficients and store in \a
207  * (this)->m_coeffs
208  *
209  * Inputs:\n
210  *
211  * - \a inarray: array of physical quadrature points to be transformed
212  *
213  * Outputs:\n
214  *
215  * - (this)->_coeffs: updated array of expansion coefficients.
216  */
218  Array<OneD, NekDouble> &outarray)
219 {
220  if (m_base[0]->Collocation() && m_base[1]->Collocation() &&
221  m_base[2]->Collocation())
222  {
223  Vmath::Vcopy(GetNcoeffs(), &inarray[0], 1, &outarray[0], 1);
224  }
225  else
226  {
227  v_IProductWRTBase(inarray, outarray);
228 
229  // get Mass matrix inverse
230  MatrixKey masskey(StdRegions::eInvMass, DetShapeType(), *this);
231  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
232 
233  // copy inarray in case inarray == outarray
234  DNekVec in(m_ncoeffs, outarray);
235  DNekVec out(m_ncoeffs, outarray, eWrapper);
236 
237  out = (*matsys) * in;
238  }
239 }
240 
241 //---------------------------------------
242 // Inner product functions
243 //---------------------------------------
244 
245 /**
246  * \brief Calculate the inner product of inarray with respect to the
247  * basis B=base0*base1*base2 and put into outarray:
248  *
249  * \f$ \begin{array}{rcl} I_{pqr} = (\phi_{pqr}, u)_{\delta} & = &
250  * \sum_{i=0}^{nq_0} \sum_{j=0}^{nq_1} \sum_{k=0}^{nq_2} \psi_{p}^{a}
251  * (\bar \eta_{1i}) \psi_{q}^{a} (\xi_{2j}) \psi_{pr}^{b} (\xi_{3k})
252  * w_i w_j w_k u(\bar \eta_{1,i} \xi_{2,j} \xi_{3,k}) J_{i,j,k}\\ & =
253  * & \sum_{i=0}^{nq_0} \psi_p^a(\bar \eta_{1,i}) \sum_{j=0}^{nq_1}
254  * \psi_{q}^a(\xi_{2,j}) \sum_{k=0}^{nq_2} \psi_{pr}^b u(\bar
255  * \eta_{1i},\xi_{2j},\xi_{3k}) J_{i,j,k} \end{array} \f$ \n
256  *
257  * where
258  *
259  * \f$ \phi_{pqr} (\xi_1 , \xi_2 , \xi_3) = \psi_p^a (\bar \eta_1)
260  * \psi_{q}^a (\xi_2) \psi_{pr}^b (\xi_3) \f$ \n
261  *
262  * which can be implemented as \n \f$f_{pr} (\xi_{3k}) =
263  * \sum_{k=0}^{nq_3} \psi_{pr}^b u(\bar \eta_{1i},\xi_{2j},\xi_{3k})
264  * J_{i,j,k} = {\bf B_3 U} \f$ \n \f$ g_{q} (\xi_{3k}) =
265  * \sum_{j=0}^{nq_1} \psi_{q}^a (\xi_{2j}) f_{pr} (\xi_{3k}) = {\bf
266  * B_2 F} \f$ \n \f$ (\phi_{pqr}, u)_{\delta} = \sum_{k=0}^{nq_0}
267  * \psi_{p}^a (\xi_{3k}) g_{q} (\xi_{3k}) = {\bf B_1 G} \f$
268  */
270  Array<OneD, NekDouble> &outarray)
271 {
272  v_IProductWRTBase_SumFac(inarray, outarray);
273 }
274 
276  const Array<OneD, const NekDouble> &inarray,
277  Array<OneD, NekDouble> &outarray, bool multiplybyweights)
278 {
279  const int nquad0 = m_base[0]->GetNumPoints();
280  const int nquad1 = m_base[1]->GetNumPoints();
281  const int nquad2 = m_base[2]->GetNumPoints();
282  const int order0 = m_base[0]->GetNumModes();
283  const int order1 = m_base[1]->GetNumModes();
284 
285  Array<OneD, NekDouble> wsp(order0 * nquad2 * (nquad1 + order1));
286 
287  if (multiplybyweights)
288  {
289  Array<OneD, NekDouble> tmp(nquad0 * nquad1 * nquad2);
290 
291  MultiplyByQuadratureMetric(inarray, tmp);
292 
294  m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
295  tmp, outarray, wsp, true, true, true);
296  }
297  else
298  {
300  m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
301  inarray, outarray, wsp, true, true, true);
302  }
303 }
304 
305 /**
306  * @brief Calculates the inner product \f$ I_{pqr} = (u,
307  * \partial_{x_i} \phi_{pqr}) \f$.
308  *
309  * The derivative of the basis functions is performed using the chain
310  * rule in order to incorporate the geometric factors. Assuming that
311  * the basis functions are a tensor product
312  * \f$\phi_{pqr}(\eta_1,\eta_2,\eta_3) =
313  * \phi_1(\eta_1)\phi_2(\eta_2)\phi_3(\eta_3)\f$, this yields the
314  * result
315  *
316  * \f[
317  * I_{pqr} = \sum_{j=1}^3 \left(u, \frac{\partial u}{\partial \eta_j}
318  * \frac{\partial \eta_j}{\partial x_i}\right)
319  * \f]
320  *
321  * In the tetrahedral element, we must also incorporate a second set
322  * of geometric factors which incorporate the collapsed co-ordinate
323  * system, so that
324  *
325  * \f[ \frac{\partial\eta_j}{\partial x_i} = \sum_{k=1}^3
326  * \frac{\partial\eta_j}{\partial\xi_k}\frac{\partial\xi_k}{\partial
327  * x_i} \f]
328  *
329  * These derivatives can be found on p152 of Sherwin & Karniadakis.
330  *
331  * @param dir Direction in which to take the derivative.
332  * @param inarray The function \f$ u \f$.
333  * @param outarray Value of the inner product.
334  */
336  const int dir, const Array<OneD, const NekDouble> &inarray,
337  Array<OneD, NekDouble> &outarray)
338 {
339  v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
340 }
341 
343  const int dir, const Array<OneD, const NekDouble> &inarray,
344  Array<OneD, NekDouble> &outarray)
345 {
346  const int nquad0 = m_base[0]->GetNumPoints();
347  const int nquad1 = m_base[1]->GetNumPoints();
348  const int nquad2 = m_base[2]->GetNumPoints();
349  const int order0 = m_base[0]->GetNumModes();
350  const int order1 = m_base[1]->GetNumModes();
351  const int nqtot = nquad0 * nquad1 * nquad2;
352 
353  Array<OneD, NekDouble> tmp1(nqtot);
354  Array<OneD, NekDouble> tmp2(nqtot);
355  Array<OneD, NekDouble> tmp3(nqtot);
356  Array<OneD, NekDouble> tmp4(nqtot);
358  Array<OneD, NekDouble> wsp(order0 * nquad2 * (nquad1 + order1));
359 
360  MultiplyByQuadratureMetric(inarray, tmp1);
361 
363  tmp2D[0] = tmp2;
364  tmp2D[1] = tmp3;
365  tmp2D[2] = tmp4;
366 
367  PrismExp::v_AlignVectorToCollapsedDir(dir, tmp1, tmp2D);
368 
369  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
370  m_base[2]->GetBdata(), tmp2, outarray, wsp,
371  true, true, true);
372 
373  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
374  m_base[2]->GetBdata(), tmp3, tmp6, wsp, true,
375  true, true);
376 
377  Vmath::Vadd(m_ncoeffs, tmp6, 1, outarray, 1, outarray, 1);
378 
379  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(),
380  m_base[2]->GetDbdata(), tmp4, tmp6, wsp, true,
381  true, true);
382 
383  Vmath::Vadd(m_ncoeffs, tmp6, 1, outarray, 1, outarray, 1);
384 }
385 
387  const int dir, const Array<OneD, const NekDouble> &inarray,
388  Array<OneD, Array<OneD, NekDouble>> &outarray)
389 {
390  const int nquad0 = m_base[0]->GetNumPoints();
391  const int nquad1 = m_base[1]->GetNumPoints();
392  const int nquad2 = m_base[2]->GetNumPoints();
393  const int order0 = m_base[0]->GetNumModes();
394  const int order1 = m_base[1]->GetNumModes();
395  const int nqtot = nquad0 * nquad1 * nquad2;
396 
397  const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
398  const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
399 
400  Array<OneD, NekDouble> gfac0(nquad0);
401  Array<OneD, NekDouble> gfac2(nquad2);
402  Array<OneD, NekDouble> tmp1(nqtot);
403  Array<OneD, NekDouble> tmp5(nqtot);
405  Array<OneD, NekDouble> wsp(order0 * nquad2 * (nquad1 + order1));
406 
407  Array<OneD, NekDouble> tmp2 = outarray[0];
408  Array<OneD, NekDouble> tmp3 = outarray[1];
409  Array<OneD, NekDouble> tmp4 = outarray[2];
410 
411  const Array<TwoD, const NekDouble> &df =
412  m_metricinfo->GetDerivFactors(GetPointsKeys());
413 
414  Vmath::Vcopy(nqtot, inarray, 1, tmp1, 1); // Dir3 metric
415 
416  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
417  {
418  Vmath::Vmul(nqtot, &df[3 * dir][0], 1, tmp1.get(), 1, tmp2.get(), 1);
419  Vmath::Vmul(nqtot, &df[3 * dir + 1][0], 1, tmp1.get(), 1, tmp3.get(),
420  1);
421  Vmath::Vmul(nqtot, &df[3 * dir + 2][0], 1, tmp1.get(), 1, tmp4.get(),
422  1);
423  }
424  else
425  {
426  Vmath::Smul(nqtot, df[3 * dir][0], tmp1.get(), 1, tmp2.get(), 1);
427  Vmath::Smul(nqtot, df[3 * dir + 1][0], tmp1.get(), 1, tmp3.get(), 1);
428  Vmath::Smul(nqtot, df[3 * dir + 2][0], tmp1.get(), 1, tmp4.get(), 1);
429  }
430 
431  // set up geometric factor: (1+z0)/2
432  for (int i = 0; i < nquad0; ++i)
433  {
434  gfac0[i] = 0.5 * (1 + z0[i]);
435  }
436 
437  // Set up geometric factor: 2/(1-z2)
438  for (int i = 0; i < nquad2; ++i)
439  {
440  gfac2[i] = 2.0 / (1 - z2[i]);
441  }
442 
443  const int nq01 = nquad0 * nquad1;
444 
445  for (int i = 0; i < nquad2; ++i)
446  {
447  Vmath::Smul(nq01, gfac2[i], &tmp2[0] + i * nq01, 1, &tmp2[0] + i * nq01,
448  1);
449  Vmath::Smul(nq01, gfac2[i], &tmp4[0] + i * nq01, 1, &tmp5[0] + i * nq01,
450  1);
451  }
452 
453  for (int i = 0; i < nquad1 * nquad2; ++i)
454  {
455  Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp5[0] + i * nquad0, 1,
456  &tmp5[0] + i * nquad0, 1);
457  }
458 
459  Vmath::Vadd(nqtot, &tmp2[0], 1, &tmp5[0], 1, &tmp2[0], 1);
460 }
461 
462 //---------------------------------------
463 // Evaluation functions
464 //---------------------------------------
465 
467 {
469  m_base[0]->GetBasisKey(), m_base[1]->GetBasisKey(),
470  m_base[2]->GetBasisKey());
471 }
472 
474 {
476  m_base[0]->GetPointsKey());
478  m_base[1]->GetPointsKey());
480  m_base[2]->GetPointsKey());
481 
483  bkey0, bkey1, bkey2);
484 }
485 
486 /**
487  * @brief Get the coordinates #coords at the local coordinates
488  * #Lcoords.
489  */
491  Array<OneD, NekDouble> &coords)
492 {
493  int i;
494 
495  ASSERTL1(Lcoords[0] <= -1.0 && Lcoords[0] >= 1.0 && Lcoords[1] <= -1.0 &&
496  Lcoords[1] >= 1.0 && Lcoords[2] <= -1.0 && Lcoords[2] >= 1.0,
497  "Local coordinates are not in region [-1,1]");
498 
499  m_geom->FillGeom();
500 
501  for (i = 0; i < m_geom->GetCoordim(); ++i)
502  {
503  coords[i] = m_geom->GetCoord(i, Lcoords);
504  }
505 }
506 
508  Array<OneD, NekDouble> &coords_1,
509  Array<OneD, NekDouble> &coords_2)
510 {
511  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
512 }
513 
514 /**
515  * Given the local cartesian coordinate \a Lcoord evaluate the
516  * value of physvals at this point by calling through to the
517  * StdExpansion method
518  */
520  const Array<OneD, const NekDouble> &Lcoord,
521  const Array<OneD, const NekDouble> &physvals)
522 {
523  // Evaluate point in local (eta) coordinates.
524  return StdPrismExp::v_PhysEvaluate(Lcoord, physvals);
525 }
526 
528  const Array<OneD, const NekDouble> &physvals)
529 {
530  Array<OneD, NekDouble> Lcoord(3);
531 
532  ASSERTL0(m_geom, "m_geom not defined");
533 
534  m_geom->GetLocCoords(coord, Lcoord);
535 
536  return StdPrismExp::v_PhysEvaluate(Lcoord, physvals);
537 }
538 
539 //---------------------------------------
540 // Helper functions
541 //---------------------------------------
542 
544 {
545  return m_geom->GetCoordim();
546 }
547 
549  const NekDouble *data, const std::vector<unsigned int> &nummodes,
550  const int mode_offset, NekDouble *coeffs,
551  std::vector<LibUtilities::BasisType> &fromType)
552 {
553  boost::ignore_unused(fromType);
554 
555  int data_order0 = nummodes[mode_offset];
556  int fillorder0 = min(m_base[0]->GetNumModes(), data_order0);
557  int data_order1 = nummodes[mode_offset + 1];
558  int order1 = m_base[1]->GetNumModes();
559  int fillorder1 = min(order1, data_order1);
560  int data_order2 = nummodes[mode_offset + 2];
561  int order2 = m_base[2]->GetNumModes();
562  int fillorder2 = min(order2, data_order2);
563 
564  switch (m_base[0]->GetBasisType())
565  {
567  {
568  int i, j;
569  int cnt = 0;
570  int cnt1 = 0;
571 
573  "Extraction routine not set up for this basis");
575  "Extraction routine not set up for this basis");
576 
577  Vmath::Zero(m_ncoeffs, coeffs, 1);
578  for (j = 0; j < fillorder0; ++j)
579  {
580  for (i = 0; i < fillorder1; ++i)
581  {
582  Vmath::Vcopy(fillorder2 - j, &data[cnt], 1, &coeffs[cnt1],
583  1);
584  cnt += data_order2 - j;
585  cnt1 += order2 - j;
586  }
587 
588  // count out data for j iteration
589  for (i = fillorder1; i < data_order1; ++i)
590  {
591  cnt += data_order2 - j;
592  }
593 
594  for (i = fillorder1; i < order1; ++i)
595  {
596  cnt1 += order2 - j;
597  }
598  }
599  }
600  break;
601  default:
602  ASSERTL0(false, "basis is either not set up or not "
603  "hierarchicial");
604  }
605 }
606 
607 void PrismExp::v_GetTracePhysMap(const int face, Array<OneD, int> &outarray)
608 {
609  int nquad0 = m_base[0]->GetNumPoints();
610  int nquad1 = m_base[1]->GetNumPoints();
611  int nquad2 = m_base[2]->GetNumPoints();
612  int nq0 = 0;
613  int nq1 = 0;
614 
615  switch (face)
616  {
617  case 0:
618  nq0 = nquad0;
619  nq1 = nquad1;
620  if (outarray.size() != nq0 * nq1)
621  {
622  outarray = Array<OneD, int>(nq0 * nq1);
623  }
624 
625  // Directions A and B positive
626  for (int i = 0; i < nquad0 * nquad1; ++i)
627  {
628  outarray[i] = i;
629  }
630  break;
631  case 1:
632 
633  nq0 = nquad0;
634  nq1 = nquad2;
635  if (outarray.size() != nq0 * nq1)
636  {
637  outarray = Array<OneD, int>(nq0 * nq1);
638  }
639 
640  // Direction A and B positive
641  for (int k = 0; k < nquad2; k++)
642  {
643  for (int i = 0; i < nquad0; ++i)
644  {
645  outarray[k * nquad0 + i] = (nquad0 * nquad1 * k) + i;
646  }
647  }
648 
649  break;
650  case 2:
651 
652  nq0 = nquad1;
653  nq1 = nquad2;
654  if (outarray.size() != nq0 * nq1)
655  {
656  outarray = Array<OneD, int>(nq0 * nq1);
657  }
658 
659  // Directions A and B positive
660  for (int j = 0; j < nquad1 * nquad2; ++j)
661  {
662  outarray[j] = nquad0 - 1 + j * nquad0;
663  }
664  break;
665  case 3:
666  nq0 = nquad0;
667  nq1 = nquad2;
668  if (outarray.size() != nq0 * nq1)
669  {
670  outarray = Array<OneD, int>(nq0 * nq1);
671  }
672 
673  // Direction A and B positive
674  for (int k = 0; k < nquad2; k++)
675  {
676  for (int i = 0; i < nquad0; ++i)
677  {
678  outarray[k * nquad0 + i] =
679  nquad0 * (nquad1 - 1) + (nquad0 * nquad1 * k) + i;
680  }
681  }
682  break;
683  case 4:
684 
685  nq0 = nquad1;
686  nq1 = nquad2;
687  if (outarray.size() != nq0 * nq1)
688  {
689  outarray = Array<OneD, int>(nq0 * nq1);
690  }
691 
692  // Directions A and B positive
693  for (int j = 0; j < nquad1 * nquad2; ++j)
694  {
695  outarray[j] = j * nquad0;
696  }
697  break;
698  default:
699  ASSERTL0(false, "face value (> 4) is out of range");
700  break;
701  }
702 }
703 
704 /** \brief Get the normals along specficied face
705  * Get the face normals interplated to a points0 x points 0
706  * type distribution
707  **/
708 void PrismExp::v_ComputeTraceNormal(const int face)
709 {
710  const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
711  GetGeom()->GetMetricInfo();
712 
714  for (int i = 0; i < ptsKeys.size(); ++i)
715  {
716  // Need at least 2 points for computing normals
717  if (ptsKeys[i].GetNumPoints() == 1)
718  {
719  LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
720  ptsKeys[i] = pKey;
721  }
722  }
723 
724  SpatialDomains::GeomType type = geomFactors->GetGtype();
725  const Array<TwoD, const NekDouble> &df =
726  geomFactors->GetDerivFactors(ptsKeys);
727  const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
728 
729  int nq0 = ptsKeys[0].GetNumPoints();
730  int nq1 = ptsKeys[1].GetNumPoints();
731  int nq2 = ptsKeys[2].GetNumPoints();
732  int nq01 = nq0 * nq1;
733  int nqtot;
734 
735  LibUtilities::BasisKey tobasis0 = GetTraceBasisKey(face, 0);
736  LibUtilities::BasisKey tobasis1 = GetTraceBasisKey(face, 1);
737 
738  // Number of quadrature points in face expansion.
739  int nq_face = tobasis0.GetNumPoints() * tobasis1.GetNumPoints();
740 
741  int vCoordDim = GetCoordim();
742  int i;
743 
746  for (i = 0; i < vCoordDim; ++i)
747  {
748  normal[i] = Array<OneD, NekDouble>(nq_face);
749  }
750 
751  size_t nqb = nq_face;
752  size_t nbnd = face;
755 
756  // Regular geometry case
757  if (type == SpatialDomains::eRegular ||
759  {
760  NekDouble fac;
761  // Set up normals
762  switch (face)
763  {
764  case 0:
765  {
766  for (i = 0; i < vCoordDim; ++i)
767  {
768  normal[i][0] = -df[3 * i + 2][0];
769  ;
770  }
771  break;
772  }
773  case 1:
774  {
775  for (i = 0; i < vCoordDim; ++i)
776  {
777  normal[i][0] = -df[3 * i + 1][0];
778  }
779  break;
780  }
781  case 2:
782  {
783  for (i = 0; i < vCoordDim; ++i)
784  {
785  normal[i][0] = df[3 * i][0] + df[3 * i + 2][0];
786  }
787  break;
788  }
789  case 3:
790  {
791  for (i = 0; i < vCoordDim; ++i)
792  {
793  normal[i][0] = df[3 * i + 1][0];
794  }
795  break;
796  }
797  case 4:
798  {
799  for (i = 0; i < vCoordDim; ++i)
800  {
801  normal[i][0] = -df[3 * i][0];
802  }
803  break;
804  }
805  default:
806  ASSERTL0(false, "face is out of range (face < 4)");
807  }
808 
809  // Normalise resulting vector.
810  fac = 0.0;
811  for (i = 0; i < vCoordDim; ++i)
812  {
813  fac += normal[i][0] * normal[i][0];
814  }
815  fac = 1.0 / sqrt(fac);
816 
817  Vmath::Fill(nqb, fac, length, 1);
818 
819  for (i = 0; i < vCoordDim; ++i)
820  {
821  Vmath::Fill(nq_face, fac * normal[i][0], normal[i], 1);
822  }
823  }
824  else
825  {
826  // Set up deformed normals.
827  int j, k;
828 
829  // Determine number of quadrature points on the face of 3D elmt
830  if (face == 0)
831  {
832  nqtot = nq0 * nq1;
833  }
834  else if (face == 1 || face == 3)
835  {
836  nqtot = nq0 * nq2;
837  }
838  else
839  {
840  nqtot = nq1 * nq2;
841  }
842 
843  LibUtilities::PointsKey points0;
844  LibUtilities::PointsKey points1;
845 
846  Array<OneD, NekDouble> faceJac(nqtot);
847  Array<OneD, NekDouble> normals(vCoordDim * nqtot, 0.0);
848 
849  // Extract Jacobian along face and recover local derivatives
850  // (dx/dr) for polynomial interpolation by multiplying m_gmat by
851  // jacobian
852  switch (face)
853  {
854  case 0:
855  {
856  for (j = 0; j < nq01; ++j)
857  {
858  normals[j] = -df[2][j] * jac[j];
859  normals[nqtot + j] = -df[5][j] * jac[j];
860  normals[2 * nqtot + j] = -df[8][j] * jac[j];
861  faceJac[j] = jac[j];
862  }
863 
864  points0 = ptsKeys[0];
865  points1 = ptsKeys[1];
866  break;
867  }
868 
869  case 1:
870  {
871  for (j = 0; j < nq0; ++j)
872  {
873  for (k = 0; k < nq2; ++k)
874  {
875  int tmp = j + nq01 * k;
876  normals[j + k * nq0] = -df[1][tmp] * jac[tmp];
877  normals[nqtot + j + k * nq0] = -df[4][tmp] * jac[tmp];
878  normals[2 * nqtot + j + k * nq0] =
879  -df[7][tmp] * jac[tmp];
880  faceJac[j + k * nq0] = jac[tmp];
881  }
882  }
883 
884  points0 = ptsKeys[0];
885  points1 = ptsKeys[2];
886  break;
887  }
888 
889  case 2:
890  {
891  for (j = 0; j < nq1; ++j)
892  {
893  for (k = 0; k < nq2; ++k)
894  {
895  int tmp = nq0 - 1 + nq0 * j + nq01 * k;
896  normals[j + k * nq1] =
897  (df[0][tmp] + df[2][tmp]) * jac[tmp];
898  normals[nqtot + j + k * nq1] =
899  (df[3][tmp] + df[5][tmp]) * jac[tmp];
900  normals[2 * nqtot + j + k * nq1] =
901  (df[6][tmp] + df[8][tmp]) * jac[tmp];
902  faceJac[j + k * nq1] = jac[tmp];
903  }
904  }
905 
906  points0 = ptsKeys[1];
907  points1 = ptsKeys[2];
908  break;
909  }
910 
911  case 3:
912  {
913  for (j = 0; j < nq0; ++j)
914  {
915  for (k = 0; k < nq2; ++k)
916  {
917  int tmp = nq0 * (nq1 - 1) + j + nq01 * k;
918  normals[j + k * nq0] = df[1][tmp] * jac[tmp];
919  normals[nqtot + j + k * nq0] = df[4][tmp] * jac[tmp];
920  normals[2 * nqtot + j + k * nq0] =
921  df[7][tmp] * jac[tmp];
922  faceJac[j + k * nq0] = jac[tmp];
923  }
924  }
925 
926  points0 = ptsKeys[0];
927  points1 = ptsKeys[2];
928  break;
929  }
930 
931  case 4:
932  {
933  for (j = 0; j < nq1; ++j)
934  {
935  for (k = 0; k < nq2; ++k)
936  {
937  int tmp = j * nq0 + nq01 * k;
938  normals[j + k * nq1] = -df[0][tmp] * jac[tmp];
939  normals[nqtot + j + k * nq1] = -df[3][tmp] * jac[tmp];
940  normals[2 * nqtot + j + k * nq1] =
941  -df[6][tmp] * jac[tmp];
942  faceJac[j + k * nq1] = jac[tmp];
943  }
944  }
945 
946  points0 = ptsKeys[1];
947  points1 = ptsKeys[2];
948  break;
949  }
950 
951  default:
952  ASSERTL0(false, "face is out of range (face < 4)");
953  }
954 
955  Array<OneD, NekDouble> work(nq_face, 0.0);
956  // Interpolate Jacobian and invert
957  LibUtilities::Interp2D(points0, points1, faceJac,
958  tobasis0.GetPointsKey(), tobasis1.GetPointsKey(),
959  work);
960  Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
961 
962  // Interpolate normal and multiply by inverse Jacobian.
963  for (i = 0; i < vCoordDim; ++i)
964  {
965  LibUtilities::Interp2D(points0, points1, &normals[i * nqtot],
966  tobasis0.GetPointsKey(),
967  tobasis1.GetPointsKey(), &normal[i][0]);
968  Vmath::Vmul(nq_face, work, 1, normal[i], 1, normal[i], 1);
969  }
970 
971  // Normalise to obtain unit normals.
972  Vmath::Zero(nq_face, work, 1);
973  for (i = 0; i < GetCoordim(); ++i)
974  {
975  Vmath::Vvtvp(nq_face, normal[i], 1, normal[i], 1, work, 1, work, 1);
976  }
977 
978  Vmath::Vsqrt(nq_face, work, 1, work, 1);
979  Vmath::Sdiv(nq_face, 1.0, work, 1, work, 1);
980 
981  Vmath::Vcopy(nqb, work, 1, length, 1);
982 
983  for (i = 0; i < GetCoordim(); ++i)
984  {
985  Vmath::Vmul(nq_face, normal[i], 1, work, 1, normal[i], 1);
986  }
987  }
988 }
989 
991  Array<OneD, NekDouble> &outarray,
992  const StdRegions::StdMatrixKey &mkey)
993 {
994  StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
995 }
996 
998  Array<OneD, NekDouble> &outarray,
999  const StdRegions::StdMatrixKey &mkey)
1000 {
1001  PrismExp::LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
1002 }
1003 
1004 void PrismExp::v_LaplacianMatrixOp(const int k1, const int k2,
1005  const Array<OneD, const NekDouble> &inarray,
1006  Array<OneD, NekDouble> &outarray,
1007  const StdRegions::StdMatrixKey &mkey)
1008 {
1009  StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1010 }
1011 
1013  Array<OneD, NekDouble> &outarray,
1014  const StdRegions::StdMatrixKey &mkey)
1015 {
1016  PrismExp::v_HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
1017 }
1018 
1020  const Array<OneD, const NekDouble> &inarray,
1021  Array<OneD, NekDouble> &outarray, const StdRegions::StdMatrixKey &mkey)
1022 {
1023  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1024 
1025  if (inarray.get() == outarray.get())
1026  {
1028  Vmath::Vcopy(m_ncoeffs, inarray.get(), 1, tmp.get(), 1);
1029 
1030  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, mat->Scale(),
1031  (mat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
1032  tmp.get(), 1, 0.0, outarray.get(), 1);
1033  }
1034  else
1035  {
1036  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, mat->Scale(),
1037  (mat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
1038  inarray.get(), 1, 0.0, outarray.get(), 1);
1039  }
1040 }
1041 
1043  const StdRegions::StdMatrixKey &mkey)
1044 {
1045  int nq = GetTotPoints();
1046 
1047  // Calculate sqrt of the Jacobian
1049  Array<OneD, NekDouble> sqrt_jac(nq);
1050  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1051  {
1052  Vmath::Vsqrt(nq, jac, 1, sqrt_jac, 1);
1053  }
1054  else
1055  {
1056  Vmath::Fill(nq, sqrt(jac[0]), sqrt_jac, 1);
1057  }
1058 
1059  // Multiply array by sqrt(Jac)
1060  Vmath::Vmul(nq, sqrt_jac, 1, array, 1, array, 1);
1061 
1062  // Apply std region filter
1063  StdPrismExp::v_SVVLaplacianFilter(array, mkey);
1064 
1065  // Divide by sqrt(Jac)
1066  Vmath::Vdiv(nq, array, 1, sqrt_jac, 1, array, 1);
1067 }
1068 
1069 //---------------------------------------
1070 // Matrix creation functions
1071 //---------------------------------------
1072 
1074 {
1075  DNekMatSharedPtr returnval;
1076 
1077  switch (mkey.GetMatrixType())
1078  {
1086  returnval = Expansion3D::v_GenMatrix(mkey);
1087  break;
1088  default:
1089  returnval = StdPrismExp::v_GenMatrix(mkey);
1090  break;
1091  }
1092 
1093  return returnval;
1094 }
1095 
1097  const StdRegions::StdMatrixKey &mkey)
1098 {
1099  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1100  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1101  LibUtilities::BasisKey bkey2 = m_base[2]->GetBasisKey();
1103  MemoryManager<StdPrismExp>::AllocateSharedPtr(bkey0, bkey1, bkey2);
1104 
1105  return tmp->GetStdMatrix(mkey);
1106 }
1107 
1109 {
1110  return m_matrixManager[mkey];
1111 }
1112 
1114  const MatrixKey &mkey)
1115 {
1116  return m_staticCondMatrixManager[mkey];
1117 }
1118 
1120 {
1121  m_staticCondMatrixManager.DeleteObject(mkey);
1122 }
1123 
1124 /**
1125  * @brief Calculate the Laplacian multiplication in a matrix-free
1126  * manner.
1127  *
1128  * This function is the kernel of the Laplacian matrix-free operator,
1129  * and is used in #v_HelmholtzMatrixOp_MatFree to determine the effect
1130  * of the Helmholtz operator in a similar fashion.
1131  *
1132  * The majority of the calculation is precisely the same as in the
1133  * hexahedral expansion; however the collapsed co-ordinate system must
1134  * be taken into account when constructing the geometric factors. How
1135  * this is done is detailed more exactly in the tetrahedral expansion.
1136  * On entry to this function, the input #inarray must be in its
1137  * backwards-transformed state (i.e. \f$\mathbf{u} =
1138  * \mathbf{B}\hat{\mathbf{u}}\f$). The output is in coefficient space.
1139  *
1140  * @see %TetExp::v_HelmholtzMatrixOp_MatFree
1141  */
1143  const Array<OneD, const NekDouble> &inarray,
1145 {
1146  int nquad0 = m_base[0]->GetNumPoints();
1147  int nquad1 = m_base[1]->GetNumPoints();
1148  int nquad2 = m_base[2]->GetNumPoints();
1149  int nqtot = nquad0 * nquad1 * nquad2;
1150  int i;
1151 
1152  // Set up temporary storage.
1153  Array<OneD, NekDouble> alloc(11 * nqtot, 0.0);
1154  Array<OneD, NekDouble> wsp1(alloc); // TensorDeriv 1
1155  Array<OneD, NekDouble> wsp2(alloc + 1 * nqtot); // TensorDeriv 2
1156  Array<OneD, NekDouble> wsp3(alloc + 2 * nqtot); // TensorDeriv 3
1157  Array<OneD, NekDouble> g0(alloc + 3 * nqtot); // g0
1158  Array<OneD, NekDouble> g1(alloc + 4 * nqtot); // g1
1159  Array<OneD, NekDouble> g2(alloc + 5 * nqtot); // g2
1160  Array<OneD, NekDouble> g3(alloc + 6 * nqtot); // g3
1161  Array<OneD, NekDouble> g4(alloc + 7 * nqtot); // g4
1162  Array<OneD, NekDouble> g5(alloc + 8 * nqtot); // g5
1163  Array<OneD, NekDouble> h0(alloc + 3 * nqtot); // h0 == g0
1164  Array<OneD, NekDouble> h1(alloc + 6 * nqtot); // h1 == g3
1165  Array<OneD, NekDouble> wsp4(alloc + 4 * nqtot); // wsp4 == g1
1166  Array<OneD, NekDouble> wsp5(alloc + 5 * nqtot); // wsp5 == g2
1167  Array<OneD, NekDouble> wsp6(alloc + 8 * nqtot); // wsp6 == g5
1168  Array<OneD, NekDouble> wsp7(alloc + 3 * nqtot); // wsp7 == g0
1169  Array<OneD, NekDouble> wsp8(alloc + 9 * nqtot); // wsp8
1170  Array<OneD, NekDouble> wsp9(alloc + 10 * nqtot); // wsp9
1171 
1172  const Array<OneD, const NekDouble> &base0 = m_base[0]->GetBdata();
1173  const Array<OneD, const NekDouble> &base1 = m_base[1]->GetBdata();
1174  const Array<OneD, const NekDouble> &base2 = m_base[2]->GetBdata();
1175  const Array<OneD, const NekDouble> &dbase0 = m_base[0]->GetDbdata();
1176  const Array<OneD, const NekDouble> &dbase1 = m_base[1]->GetDbdata();
1177  const Array<OneD, const NekDouble> &dbase2 = m_base[2]->GetDbdata();
1178 
1179  // Step 1. LAPLACIAN MATRIX OPERATION
1180  // wsp1 = du_dxi1 = D_xi1 * wsp0 = D_xi1 * u
1181  // wsp2 = du_dxi2 = D_xi2 * wsp0 = D_xi2 * u
1182  // wsp3 = du_dxi3 = D_xi3 * wsp0 = D_xi3 * u
1183  StdExpansion3D::PhysTensorDeriv(inarray, wsp1, wsp2, wsp3);
1184 
1185  const Array<TwoD, const NekDouble> &df =
1186  m_metricinfo->GetDerivFactors(GetPointsKeys());
1187  const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
1188  const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
1189 
1190  // Step 2. Calculate the metric terms of the collapsed
1191  // coordinate transformation (Spencer's book P152)
1192  for (i = 0; i < nquad2; ++i)
1193  {
1194  Vmath::Fill(nquad0 * nquad1, 2.0 / (1.0 - z2[i]),
1195  &h0[0] + i * nquad0 * nquad1, 1);
1196  Vmath::Fill(nquad0 * nquad1, 2.0 / (1.0 - z2[i]),
1197  &h1[0] + i * nquad0 * nquad1, 1);
1198  }
1199  for (i = 0; i < nquad0; i++)
1200  {
1201  Blas::Dscal(nquad1 * nquad2, 0.5 * (1 + z0[i]), &h1[0] + i, nquad0);
1202  }
1203 
1204  // Step 3. Construct combined metric terms for physical space to
1205  // collapsed coordinate system. Order of construction optimised
1206  // to minimise temporary storage
1207  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1208  {
1209  // wsp4 = d eta_1/d x_1
1210  Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &df[2][0], 1, &h1[0], 1,
1211  &wsp4[0], 1);
1212  // wsp5 = d eta_2/d x_1
1213  Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &df[5][0], 1, &h1[0], 1,
1214  &wsp5[0], 1);
1215  // wsp6 = d eta_3/d x_1d
1216  Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &df[8][0], 1, &h1[0], 1,
1217  &wsp6[0], 1);
1218 
1219  // g0 (overwrites h0)
1220  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1221  1, &g0[0], 1);
1222  Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1223 
1224  // g3 (overwrites h1)
1225  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &wsp4[0], 1, &df[4][0], 1, &wsp5[0],
1226  1, &g3[0], 1);
1227  Vmath::Vvtvp(nqtot, &df[7][0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1228 
1229  // g4
1230  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0],
1231  1, &g4[0], 1);
1232  Vmath::Vvtvp(nqtot, &df[8][0], 1, &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1233 
1234  // Overwrite wsp4/5/6 with g1/2/5
1235  // g1
1236  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &df[1][0], 1, &df[4][0], 1,
1237  &df[4][0], 1, &g1[0], 1);
1238  Vmath::Vvtvp(nqtot, &df[7][0], 1, &df[7][0], 1, &g1[0], 1, &g1[0], 1);
1239 
1240  // g2
1241  Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &df[2][0], 1, &df[5][0], 1,
1242  &df[5][0], 1, &g2[0], 1);
1243  Vmath::Vvtvp(nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1244 
1245  // g5
1246  Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &df[2][0], 1, &df[4][0], 1,
1247  &df[5][0], 1, &g5[0], 1);
1248  Vmath::Vvtvp(nqtot, &df[7][0], 1, &df[8][0], 1, &g5[0], 1, &g5[0], 1);
1249  }
1250  else
1251  {
1252  // wsp4 = d eta_1/d x_1
1253  Vmath::Svtsvtp(nqtot, df[0][0], &h0[0], 1, df[2][0], &h1[0], 1,
1254  &wsp4[0], 1);
1255  // wsp5 = d eta_2/d x_1
1256  Vmath::Svtsvtp(nqtot, df[3][0], &h0[0], 1, df[5][0], &h1[0], 1,
1257  &wsp5[0], 1);
1258  // wsp6 = d eta_3/d x_1
1259  Vmath::Svtsvtp(nqtot, df[6][0], &h0[0], 1, df[8][0], &h1[0], 1,
1260  &wsp6[0], 1);
1261 
1262  // g0 (overwrites h0)
1263  Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1264  1, &g0[0], 1);
1265  Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1266 
1267  // g3 (overwrites h1)
1268  Vmath::Svtsvtp(nqtot, df[1][0], &wsp4[0], 1, df[4][0], &wsp5[0], 1,
1269  &g3[0], 1);
1270  Vmath::Svtvp(nqtot, df[7][0], &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
1271 
1272  // g4
1273  Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1,
1274  &g4[0], 1);
1275  Vmath::Svtvp(nqtot, df[8][0], &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1276 
1277  // Overwrite wsp4/5/6 with g1/2/5
1278  // g1
1279  Vmath::Fill(nqtot,
1280  df[1][0] * df[1][0] + df[4][0] * df[4][0] +
1281  df[7][0] * df[7][0],
1282  &g1[0], 1);
1283 
1284  // g2
1285  Vmath::Fill(nqtot,
1286  df[2][0] * df[2][0] + df[5][0] * df[5][0] +
1287  df[8][0] * df[8][0],
1288  &g2[0], 1);
1289 
1290  // g5
1291  Vmath::Fill(nqtot,
1292  df[1][0] * df[2][0] + df[4][0] * df[5][0] +
1293  df[7][0] * df[8][0],
1294  &g5[0], 1);
1295  }
1296  // Compute component derivatives into wsp7, 8, 9 (wsp7 overwrites
1297  // g0).
1298  Vmath::Vvtvvtp(nqtot, &g0[0], 1, &wsp1[0], 1, &g3[0], 1, &wsp2[0], 1,
1299  &wsp7[0], 1);
1300  Vmath::Vvtvp(nqtot, &g4[0], 1, &wsp3[0], 1, &wsp7[0], 1, &wsp7[0], 1);
1301  Vmath::Vvtvvtp(nqtot, &g1[0], 1, &wsp2[0], 1, &g3[0], 1, &wsp1[0], 1,
1302  &wsp8[0], 1);
1303  Vmath::Vvtvp(nqtot, &g5[0], 1, &wsp3[0], 1, &wsp8[0], 1, &wsp8[0], 1);
1304  Vmath::Vvtvvtp(nqtot, &g2[0], 1, &wsp3[0], 1, &g4[0], 1, &wsp1[0], 1,
1305  &wsp9[0], 1);
1306  Vmath::Vvtvp(nqtot, &g5[0], 1, &wsp2[0], 1, &wsp9[0], 1, &wsp9[0], 1);
1307 
1308  // Step 4.
1309  // Multiply by quadrature metric
1310  MultiplyByQuadratureMetric(wsp7, wsp7);
1311  MultiplyByQuadratureMetric(wsp8, wsp8);
1312  MultiplyByQuadratureMetric(wsp9, wsp9);
1313 
1314  // Perform inner product w.r.t derivative bases.
1315  IProductWRTBase_SumFacKernel(dbase0, base1, base2, wsp7, wsp1, wsp, false,
1316  true, true);
1317  IProductWRTBase_SumFacKernel(base0, dbase1, base2, wsp8, wsp2, wsp, true,
1318  false, true);
1319  IProductWRTBase_SumFacKernel(base0, base1, dbase2, wsp9, outarray, wsp,
1320  true, true, false);
1321 
1322  // Step 5.
1323  // Sum contributions from wsp1, wsp2 and outarray.
1324  Vmath::Vadd(m_ncoeffs, wsp1.get(), 1, outarray.get(), 1, outarray.get(), 1);
1325  Vmath::Vadd(m_ncoeffs, wsp2.get(), 1, outarray.get(), 1, outarray.get(), 1);
1326 }
1327 
1329  bool oldstandard)
1330 {
1331  boost::ignore_unused(oldstandard);
1332 
1333  int np0 = m_base[0]->GetNumPoints();
1334  int np1 = m_base[1]->GetNumPoints();
1335  int np2 = m_base[2]->GetNumPoints();
1336  int np = max(np0, max(np1, np2));
1337  Array<OneD, int> prismpt(6);
1338  bool standard = true;
1339 
1340  int vid0 = m_geom->GetVid(0);
1341  int vid1 = m_geom->GetVid(1);
1342  int vid2 = m_geom->GetVid(4);
1343  int rotate = 0;
1344 
1345  // sort out prism rotation according to
1346  if ((vid2 < vid1) && (vid2 < vid0)) // top triangle vertex is lowest id
1347  {
1348  rotate = 0;
1349  if (vid0 > vid1)
1350  {
1351  standard = false; // reverse base direction
1352  }
1353  }
1354  else if ((vid1 < vid2) && (vid1 < vid0))
1355  {
1356  rotate = 1;
1357  if (vid2 > vid0)
1358  {
1359  standard = false; // reverse base direction
1360  }
1361  }
1362  else if ((vid0 < vid2) && (vid0 < vid1))
1363  {
1364  rotate = 2;
1365  if (vid1 > vid2)
1366  {
1367  standard = false; // reverse base direction
1368  }
1369  }
1370 
1371  conn = Array<OneD, int>(12 * (np - 1) * (np - 1) * (np - 1));
1372 
1373  int row = 0;
1374  int rowp1 = 0;
1375  int plane = 0;
1376  int row1 = 0;
1377  int row1p1 = 0;
1378  int planep1 = 0;
1379  int cnt = 0;
1380 
1381  Array<OneD, int> rot(3);
1382 
1383  rot[0] = (0 + rotate) % 3;
1384  rot[1] = (1 + rotate) % 3;
1385  rot[2] = (2 + rotate) % 3;
1386 
1387  // lower diagonal along 1-3 on base
1388  for (int i = 0; i < np - 1; ++i)
1389  {
1390  planep1 += (np - i) * np;
1391  row = 0; // current plane row offset
1392  rowp1 = 0; // current plane row plus one offset
1393  row1 = 0; // next plane row offset
1394  row1p1 = 0; // nex plane row plus one offset
1395  if (standard == false)
1396  {
1397  for (int j = 0; j < np - 1; ++j)
1398  {
1399  rowp1 += np - i;
1400  row1p1 += np - i - 1;
1401  for (int k = 0; k < np - i - 2; ++k)
1402  {
1403  // bottom prism block
1404  prismpt[rot[0]] = plane + row + k;
1405  prismpt[rot[1]] = plane + row + k + 1;
1406  prismpt[rot[2]] = planep1 + row1 + k;
1407 
1408  prismpt[3 + rot[0]] = plane + rowp1 + k;
1409  prismpt[3 + rot[1]] = plane + rowp1 + k + 1;
1410  prismpt[3 + rot[2]] = planep1 + row1p1 + k;
1411 
1412  conn[cnt++] = prismpt[0];
1413  conn[cnt++] = prismpt[1];
1414  conn[cnt++] = prismpt[3];
1415  conn[cnt++] = prismpt[2];
1416 
1417  conn[cnt++] = prismpt[5];
1418  conn[cnt++] = prismpt[2];
1419  conn[cnt++] = prismpt[3];
1420  conn[cnt++] = prismpt[4];
1421 
1422  conn[cnt++] = prismpt[3];
1423  conn[cnt++] = prismpt[1];
1424  conn[cnt++] = prismpt[4];
1425  conn[cnt++] = prismpt[2];
1426 
1427  // upper prism block.
1428  prismpt[rot[0]] = planep1 + row1 + k + 1;
1429  prismpt[rot[1]] = planep1 + row1 + k;
1430  prismpt[rot[2]] = plane + row + k + 1;
1431 
1432  prismpt[3 + rot[0]] = planep1 + row1p1 + k + 1;
1433  prismpt[3 + rot[1]] = planep1 + row1p1 + k;
1434  prismpt[3 + rot[2]] = plane + rowp1 + k + 1;
1435 
1436  conn[cnt++] = prismpt[0];
1437  conn[cnt++] = prismpt[1];
1438  conn[cnt++] = prismpt[2];
1439  conn[cnt++] = prismpt[5];
1440 
1441  conn[cnt++] = prismpt[5];
1442  conn[cnt++] = prismpt[0];
1443  conn[cnt++] = prismpt[4];
1444  conn[cnt++] = prismpt[1];
1445 
1446  conn[cnt++] = prismpt[3];
1447  conn[cnt++] = prismpt[4];
1448  conn[cnt++] = prismpt[0];
1449  conn[cnt++] = prismpt[5];
1450  }
1451 
1452  // bottom prism block
1453  prismpt[rot[0]] = plane + row + np - i - 2;
1454  prismpt[rot[1]] = plane + row + np - i - 1;
1455  prismpt[rot[2]] = planep1 + row1 + np - i - 2;
1456 
1457  prismpt[3 + rot[0]] = plane + rowp1 + np - i - 2;
1458  prismpt[3 + rot[1]] = plane + rowp1 + np - i - 1;
1459  prismpt[3 + rot[2]] = planep1 + row1p1 + np - i - 2;
1460 
1461  conn[cnt++] = prismpt[0];
1462  conn[cnt++] = prismpt[1];
1463  conn[cnt++] = prismpt[3];
1464  conn[cnt++] = prismpt[2];
1465 
1466  conn[cnt++] = prismpt[5];
1467  conn[cnt++] = prismpt[2];
1468  conn[cnt++] = prismpt[3];
1469  conn[cnt++] = prismpt[4];
1470 
1471  conn[cnt++] = prismpt[3];
1472  conn[cnt++] = prismpt[1];
1473  conn[cnt++] = prismpt[4];
1474  conn[cnt++] = prismpt[2];
1475 
1476  row += np - i;
1477  row1 += np - i - 1;
1478  }
1479  }
1480  else
1481  { // lower diagonal along 0-4 on base
1482  for (int j = 0; j < np - 1; ++j)
1483  {
1484  rowp1 += np - i;
1485  row1p1 += np - i - 1;
1486  for (int k = 0; k < np - i - 2; ++k)
1487  {
1488  // bottom prism block
1489  prismpt[rot[0]] = plane + row + k;
1490  prismpt[rot[1]] = plane + row + k + 1;
1491  prismpt[rot[2]] = planep1 + row1 + k;
1492 
1493  prismpt[3 + rot[0]] = plane + rowp1 + k;
1494  prismpt[3 + rot[1]] = plane + rowp1 + k + 1;
1495  prismpt[3 + rot[2]] = planep1 + row1p1 + k;
1496 
1497  conn[cnt++] = prismpt[0];
1498  conn[cnt++] = prismpt[1];
1499  conn[cnt++] = prismpt[4];
1500  conn[cnt++] = prismpt[2];
1501 
1502  conn[cnt++] = prismpt[4];
1503  conn[cnt++] = prismpt[3];
1504  conn[cnt++] = prismpt[0];
1505  conn[cnt++] = prismpt[2];
1506 
1507  conn[cnt++] = prismpt[3];
1508  conn[cnt++] = prismpt[4];
1509  conn[cnt++] = prismpt[5];
1510  conn[cnt++] = prismpt[2];
1511 
1512  // upper prism block.
1513  prismpt[rot[0]] = planep1 + row1 + k + 1;
1514  prismpt[rot[1]] = planep1 + row1 + k;
1515  prismpt[rot[2]] = plane + row + k + 1;
1516 
1517  prismpt[3 + rot[0]] = planep1 + row1p1 + k + 1;
1518  prismpt[3 + rot[1]] = planep1 + row1p1 + k;
1519  prismpt[3 + rot[2]] = plane + rowp1 + k + 1;
1520 
1521  conn[cnt++] = prismpt[0];
1522  conn[cnt++] = prismpt[2];
1523  conn[cnt++] = prismpt[1];
1524  conn[cnt++] = prismpt[5];
1525 
1526  conn[cnt++] = prismpt[3];
1527  conn[cnt++] = prismpt[5];
1528  conn[cnt++] = prismpt[0];
1529  conn[cnt++] = prismpt[1];
1530 
1531  conn[cnt++] = prismpt[5];
1532  conn[cnt++] = prismpt[3];
1533  conn[cnt++] = prismpt[4];
1534  conn[cnt++] = prismpt[1];
1535  }
1536 
1537  // bottom prism block
1538  prismpt[rot[0]] = plane + row + np - i - 2;
1539  prismpt[rot[1]] = plane + row + np - i - 1;
1540  prismpt[rot[2]] = planep1 + row1 + np - i - 2;
1541 
1542  prismpt[3 + rot[0]] = plane + rowp1 + np - i - 2;
1543  prismpt[3 + rot[1]] = plane + rowp1 + np - i - 1;
1544  prismpt[3 + rot[2]] = planep1 + row1p1 + np - i - 2;
1545 
1546  conn[cnt++] = prismpt[0];
1547  conn[cnt++] = prismpt[1];
1548  conn[cnt++] = prismpt[4];
1549  conn[cnt++] = prismpt[2];
1550 
1551  conn[cnt++] = prismpt[4];
1552  conn[cnt++] = prismpt[3];
1553  conn[cnt++] = prismpt[0];
1554  conn[cnt++] = prismpt[2];
1555 
1556  conn[cnt++] = prismpt[3];
1557  conn[cnt++] = prismpt[4];
1558  conn[cnt++] = prismpt[5];
1559  conn[cnt++] = prismpt[2];
1560 
1561  row += np - i;
1562  row1 += np - i - 1;
1563  }
1564  }
1565  plane += (np - i) * np;
1566  }
1567 }
1568 
1569 /** @brief: This method gets all of the factors which are
1570  required as part of the Gradient Jump Penalty
1571  stabilisation and involves the product of the normal and
1572  geometric factors along the element trace.
1573 */
1575  Array<OneD, Array<OneD, NekDouble>> &d0factors,
1576  Array<OneD, Array<OneD, NekDouble>> &d1factors,
1577  Array<OneD, Array<OneD, NekDouble>> &d2factors)
1578 {
1579  int nquad0 = GetNumPoints(0);
1580  int nquad1 = GetNumPoints(1);
1581  int nquad2 = GetNumPoints(2);
1582 
1583  const Array<TwoD, const NekDouble> &df =
1584  m_metricinfo->GetDerivFactors(GetPointsKeys());
1585 
1586  if (d0factors.size() != 5)
1587  {
1588  d0factors = Array<OneD, Array<OneD, NekDouble>>(5);
1589  d1factors = Array<OneD, Array<OneD, NekDouble>>(5);
1590  d2factors = Array<OneD, Array<OneD, NekDouble>>(5);
1591  }
1592 
1593  if (d0factors[0].size() != nquad0 * nquad1)
1594  {
1595  d0factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1596  d1factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1597  d2factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1598  }
1599 
1600  if (d0factors[1].size() != nquad0 * nquad2)
1601  {
1602  d0factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1603  d0factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1604  d1factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1605  d1factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1606  d2factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1607  d2factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1608  }
1609 
1610  if (d0factors[2].size() != nquad1 * nquad2)
1611  {
1612  d0factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1613  d0factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1614  d1factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1615  d1factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1616  d2factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1617  d2factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1618  }
1619 
1620  // Outwards normals
1621  const Array<OneD, const Array<OneD, NekDouble>> &normal_0 =
1622  GetTraceNormal(0);
1623  const Array<OneD, const Array<OneD, NekDouble>> &normal_1 =
1624  GetTraceNormal(1);
1625  const Array<OneD, const Array<OneD, NekDouble>> &normal_2 =
1626  GetTraceNormal(2);
1627  const Array<OneD, const Array<OneD, NekDouble>> &normal_3 =
1628  GetTraceNormal(3);
1629  const Array<OneD, const Array<OneD, NekDouble>> &normal_4 =
1630  GetTraceNormal(4);
1631 
1632  int ncoords = normal_0.size();
1633 
1634  // first gather together standard cartesian inner products
1635  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1636  {
1637  // face 0
1638  for (int i = 0; i < nquad0 * nquad1; ++i)
1639  {
1640  d0factors[0][i] = df[0][i] * normal_0[0][i];
1641  d1factors[0][i] = df[1][i] * normal_0[0][i];
1642  d2factors[0][i] = df[2][i] * normal_0[0][i];
1643  }
1644 
1645  for (int n = 1; n < ncoords; ++n)
1646  {
1647  for (int i = 0; i < nquad0 * nquad1; ++i)
1648  {
1649  d0factors[0][i] += df[3 * n][i] * normal_0[n][i];
1650  d1factors[0][i] += df[3 * n + 1][i] * normal_0[n][i];
1651  d2factors[0][i] += df[3 * n + 2][i] * normal_0[n][i];
1652  }
1653  }
1654 
1655  // faces 1 and 3
1656  for (int j = 0; j < nquad2; ++j)
1657  {
1658  for (int i = 0; i < nquad0; ++i)
1659  {
1660  d0factors[1][i] = df[0][j * nquad0 * nquad1 + i] *
1661  normal_1[0][j * nquad0 + i];
1662  d0factors[3][i] =
1663  df[0][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1664  normal_3[0][j * nquad0 + i];
1665  d1factors[1][i] = df[1][j * nquad0 * nquad1 + i] *
1666  normal_1[0][j * nquad0 + i];
1667  d1factors[3][i] =
1668  df[1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1669  normal_3[0][j * nquad0 + i];
1670  d2factors[1][i] = df[2][j * nquad0 * nquad1 + i] *
1671  normal_1[0][j * nquad0 + i];
1672  d2factors[3][i] =
1673  df[2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1674  normal_3[0][j * nquad0 + i];
1675  }
1676  }
1677 
1678  for (int n = 1; n < ncoords; ++n)
1679  {
1680  for (int j = 0; j < nquad2; ++j)
1681  {
1682  for (int i = 0; i < nquad0; ++i)
1683  {
1684  d0factors[1][i] = df[3 * n][j * nquad0 * nquad1 + i] *
1685  normal_1[0][j * nquad0 + i];
1686  d0factors[3][i] =
1687  df[3 * n][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1688  normal_3[0][j * nquad0 + i];
1689  d1factors[1][i] = df[3 * n + 1][j * nquad0 * nquad1 + i] *
1690  normal_1[0][j * nquad0 + i];
1691  d1factors[3][i] =
1692  df[3 * n + 1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1693  normal_3[0][j * nquad0 + i];
1694  d2factors[1][i] = df[3 * n + 2][j * nquad0 * nquad1 + i] *
1695  normal_1[0][j * nquad0 + i];
1696  d2factors[3][i] =
1697  df[3 * n + 2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1698  normal_3[0][j * nquad0 + i];
1699  }
1700  }
1701  }
1702 
1703  // faces 2 and 4
1704  for (int j = 0; j < nquad2; ++j)
1705  {
1706  for (int i = 0; i < nquad1; ++i)
1707  {
1708  d0factors[2][j * nquad1 + i] =
1709  df[0][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1710  normal_2[0][j * nquad1 + i];
1711  d0factors[4][j * nquad1 + i] =
1712  df[0][j * nquad0 * nquad1 + i * nquad0] *
1713  normal_4[0][j * nquad1 + i];
1714  d1factors[2][j * nquad1 + i] =
1715  df[1][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1716  normal_2[0][j * nquad1 + i];
1717  d1factors[4][j * nquad1 + i] =
1718  df[1][j * nquad0 * nquad1 + i * nquad0] *
1719  normal_4[0][j * nquad1 + i];
1720  d2factors[2][j * nquad1 + i] =
1721  df[2][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1722  normal_2[0][j * nquad1 + i];
1723  d2factors[4][j * nquad1 + i] =
1724  df[2][j * nquad0 * nquad1 + i * nquad0] *
1725  normal_4[0][j * nquad1 + i];
1726  }
1727  }
1728 
1729  for (int n = 1; n < ncoords; ++n)
1730  {
1731  for (int j = 0; j < nquad2; ++j)
1732  {
1733  for (int i = 0; i < nquad1; ++i)
1734  {
1735  d0factors[2][j * nquad1 + i] +=
1736  df[3 * n][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1737  normal_2[n][j * nquad0 + i];
1738  d0factors[4][j * nquad0 + i] +=
1739  df[3 * n][i * nquad0 + j * nquad0 * nquad1] *
1740  normal_4[n][j * nquad0 + i];
1741  d1factors[2][j * nquad1 + i] +=
1742  df[3 * n + 1]
1743  [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1744  normal_2[n][j * nquad0 + i];
1745  d1factors[4][j * nquad0 + i] +=
1746  df[3 * n + 1][i * nquad0 + j * nquad0 * nquad1] *
1747  normal_4[n][j * nquad0 + i];
1748  d2factors[2][j * nquad1 + i] +=
1749  df[3 * n + 2]
1750  [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1751  normal_2[n][j * nquad0 + i];
1752  d2factors[4][j * nquad0 + i] +=
1753  df[3 * n + 2][i * nquad0 + j * nquad0 * nquad1] *
1754  normal_4[n][j * nquad0 + i];
1755  }
1756  }
1757  }
1758  }
1759  else
1760  {
1761  // Face 0
1762  for (int i = 0; i < nquad0 * nquad1; ++i)
1763  {
1764  d0factors[0][i] = df[0][0] * normal_0[0][i];
1765  d1factors[0][i] = df[1][0] * normal_0[0][i];
1766  d2factors[0][i] = df[2][0] * normal_0[0][i];
1767  }
1768 
1769  for (int n = 1; n < ncoords; ++n)
1770  {
1771  for (int i = 0; i < nquad0 * nquad1; ++i)
1772  {
1773  d0factors[0][i] += df[3 * n][0] * normal_0[n][i];
1774  d1factors[0][i] += df[3 * n + 1][0] * normal_0[n][i];
1775  d2factors[0][i] += df[3 * n + 2][0] * normal_0[n][i];
1776  }
1777  }
1778 
1779  // faces 1 and 3
1780  for (int i = 0; i < nquad0 * nquad2; ++i)
1781  {
1782  d0factors[1][i] = df[0][0] * normal_1[0][i];
1783  d0factors[3][i] = df[0][0] * normal_3[0][i];
1784 
1785  d1factors[1][i] = df[1][0] * normal_1[0][i];
1786  d1factors[3][i] = df[1][0] * normal_3[0][i];
1787 
1788  d2factors[1][i] = df[2][0] * normal_1[0][i];
1789  d2factors[3][i] = df[2][0] * normal_3[0][i];
1790  }
1791 
1792  for (int n = 1; n < ncoords; ++n)
1793  {
1794  for (int i = 0; i < nquad0 * nquad2; ++i)
1795  {
1796  d0factors[1][i] += df[3 * n][0] * normal_1[n][i];
1797  d0factors[3][i] += df[3 * n][0] * normal_3[n][i];
1798 
1799  d1factors[1][i] += df[3 * n + 1][0] * normal_1[n][i];
1800  d1factors[3][i] += df[3 * n + 1][0] * normal_3[n][i];
1801 
1802  d2factors[1][i] += df[3 * n + 2][0] * normal_1[n][i];
1803  d2factors[3][i] += df[3 * n + 2][0] * normal_3[n][i];
1804  }
1805  }
1806 
1807  // faces 2 and 4
1808  for (int i = 0; i < nquad1 * nquad2; ++i)
1809  {
1810  d0factors[2][i] = df[0][0] * normal_2[0][i];
1811  d0factors[4][i] = df[0][0] * normal_4[0][i];
1812 
1813  d1factors[2][i] = df[1][0] * normal_2[0][i];
1814  d1factors[4][i] = df[1][0] * normal_4[0][i];
1815 
1816  d2factors[2][i] = df[2][0] * normal_2[0][i];
1817  d2factors[4][i] = df[2][0] * normal_4[0][i];
1818  }
1819 
1820  for (int n = 1; n < ncoords; ++n)
1821  {
1822  for (int i = 0; i < nquad1 * nquad2; ++i)
1823  {
1824  d0factors[2][i] += df[3 * n][0] * normal_2[n][i];
1825  d0factors[4][i] += df[3 * n][0] * normal_4[n][i];
1826 
1827  d1factors[2][i] += df[3 * n + 1][0] * normal_2[n][i];
1828  d1factors[4][i] += df[3 * n + 1][0] * normal_4[n][i];
1829 
1830  d2factors[2][i] += df[3 * n + 2][0] * normal_2[n][i];
1831  d2factors[4][i] += df[3 * n + 2][0] * normal_4[n][i];
1832  }
1833  }
1834  }
1835 }
1836 } // namespace LocalRegions
1837 } // 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
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:88
const NormalVector & GetTraceNormal(const int id)
Definition: Expansion.cpp:250
void v_ComputeTraceNormal(const int face)
Get the normals along specficied face Get the face normals interplated to a points0 x points 0 type d...
Definition: PrismExp.cpp:708
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: PrismExp.cpp:1119
virtual void v_NormalTraceDerivFactors(Array< OneD, Array< OneD, NekDouble >> &d0factors, Array< OneD, Array< OneD, NekDouble >> &d1factors, Array< OneD, Array< OneD, NekDouble >> &d2factors)
: This method gets all of the factors which are required as part of the Gradient Jump Penalty stabili...
Definition: PrismExp.cpp:1574
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: PrismExp.cpp:997
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
Definition: PrismExp.cpp:1042
virtual void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: PrismExp.cpp:1019
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: PrismExp.cpp:548
virtual void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true)
Definition: PrismExp.cpp:1328
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: PrismExp.cpp:527
virtual NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
Definition: PrismExp.cpp:519
virtual void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
Calculate the Laplacian multiplication in a matrix-free manner.
Definition: PrismExp.cpp:1142
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: PrismExp.cpp:1012
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: PrismExp.cpp:269
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: PrismExp.h:193
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: PrismExp.cpp:1113
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: PrismExp.cpp:1073
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
Definition: PrismExp.cpp:473
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: PrismExp.cpp:1108
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: PrismExp.cpp:507
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Calculates the inner product .
Definition: PrismExp.cpp:335
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: PrismExp.cpp:131
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: PrismExp.cpp:990
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrate the physical point list inarray over prismatic region and return the value.
Definition: PrismExp.cpp:104
PrismExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc, const SpatialDomains::PrismGeomSharedPtr &geom)
Constructor using BasisKey class for quadrature points and order definition.
Definition: PrismExp.cpp:49
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: PrismExp.cpp:275
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
Definition: PrismExp.cpp:466
virtual void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: PrismExp.cpp:386
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: PrismExp.cpp:217
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: PrismExp.cpp:342
virtual void v_GetTracePhysMap(const int face, Array< OneD, int > &outarray)
Definition: PrismExp.cpp:607
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: PrismExp.h:191
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
Get the coordinates #coords at the local coordinates #Lcoords.
Definition: PrismExp.cpp:490
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: PrismExp.cpp:1096
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)
virtual void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
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
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:375
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:213
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 Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = A x where A[m x n].
Definition: Blas.hpp:246
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
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:51
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:50
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
std::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:84
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< StdPrismExp > StdPrismExpSharedPtr
Definition: StdPrismExp.h:231
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
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