Nektar++
HexExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File HexExp.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: Methods for Hex expansion in local regoins
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
37 #include <LocalRegions/HexExp.h>
38 #include <SpatialDomains/HexGeom.h>
39 
40 using namespace std;
41 
42 namespace Nektar
43 {
44 namespace LocalRegions
45 {
46 /**
47  * @class HexExp
48  * Defines a hexahedral local expansion.
49  */
50 
51 /**
52  * \brief Constructor using BasisKey class for quadrature points and
53  * order definition
54  *
55  * @param Ba Basis key for first coordinate.
56  * @param Bb Basis key for second coordinate.
57  * @param Bc Basis key for third coordinate.
58  */
59 HexExp::HexExp(const LibUtilities::BasisKey &Ba,
60  const LibUtilities::BasisKey &Bb,
61  const LibUtilities::BasisKey &Bc,
63  : StdExpansion(Ba.GetNumModes() * Bb.GetNumModes() * Bc.GetNumModes(), 3,
64  Ba, Bb, Bc),
65  StdExpansion3D(Ba.GetNumModes() * Bb.GetNumModes() * Bc.GetNumModes(), Ba,
66  Bb, Bc),
67  StdHexExp(Ba, Bb, Bc), Expansion(geom), Expansion3D(geom),
68  m_matrixManager(
69  std::bind(&Expansion3D::CreateMatrix, this, std::placeholders::_1),
70  std::string("HexExpMatrix")),
71  m_staticCondMatrixManager(std::bind(&Expansion::CreateStaticCondMatrix,
72  this, std::placeholders::_1),
73  std::string("HexExpStaticCondMatrix"))
74 {
75 }
76 
77 /**
78  * \brief Copy Constructor
79  *
80  * @param T HexExp to copy.
81  */
83  : StdExpansion(T), StdExpansion3D(T), StdHexExp(T), Expansion(T),
84  Expansion3D(T), m_matrixManager(T.m_matrixManager),
85  m_staticCondMatrixManager(T.m_staticCondMatrixManager)
86 {
87 }
88 
89 /**
90  * \brief Destructor
91  */
93 {
94 }
95 
96 //-----------------------------
97 // Integration Methods
98 //-----------------------------
99 /**
100  * \brief Integrate the physical point list \a inarray over region
101  *
102  * @param inarray definition of function to be returned at
103  * quadrature points of expansion.
104  * @returns \f$\int^1_{-1}\int^1_{-1} \int^1_{-1}
105  * u(\eta_1, \eta_2, \eta_3) J[i,j,k] d \eta_1 d \eta_2 d \eta_3 \f$
106  * where \f$inarray[i,j,k] = u(\eta_{1i},\eta_{2j},\eta_{3k}) \f$
107  * and \f$ J[i,j,k] \f$ is the Jacobian evaluated at the quadrature
108  * point.
109  */
111 {
112  int nquad0 = m_base[0]->GetNumPoints();
113  int nquad1 = m_base[1]->GetNumPoints();
114  int nquad2 = m_base[2]->GetNumPoints();
116  NekDouble returnVal;
117  Array<OneD, NekDouble> tmp(nquad0 * nquad1 * nquad2);
118 
119  // multiply inarray with Jacobian
120 
121  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
122  {
123  Vmath::Vmul(nquad0 * nquad1 * nquad2, &jac[0], 1,
124  (NekDouble *)&inarray[0], 1, &tmp[0], 1);
125  }
126  else
127  {
128  Vmath::Smul(nquad0 * nquad1 * nquad2, (NekDouble)jac[0],
129  (NekDouble *)&inarray[0], 1, &tmp[0], 1);
130  }
131 
132  // call StdHexExp version;
133  returnVal = StdHexExp::v_Integral(tmp);
134 
135  return returnVal;
136 }
137 
138 //-----------------------------
139 // Differentiation Methods
140 //-----------------------------
141 /**
142  * \brief Calculate the derivative of the physical points
143  *
144  * For Hexahedral region can use the Tensor_Deriv function defined
145  * under StdExpansion.
146  * @param inarray Input array
147  * @param out_d0 Derivative of \a inarray in first direction.
148  * @param out_d1 Derivative of \a inarray in second direction.
149  * @param out_d2 Derivative of \a inarray in third direction.
150  */
152  Array<OneD, NekDouble> &out_d0,
153  Array<OneD, NekDouble> &out_d1,
154  Array<OneD, NekDouble> &out_d2)
155 {
156  int nquad0 = m_base[0]->GetNumPoints();
157  int nquad1 = m_base[1]->GetNumPoints();
158  int nquad2 = m_base[2]->GetNumPoints();
159  int ntot = nquad0 * nquad1 * nquad2;
160 
162  m_metricinfo->GetDerivFactors(GetPointsKeys());
166 
167  StdHexExp::v_PhysDeriv(inarray, Diff0, Diff1, Diff2);
168 
169  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
170  {
171  if (out_d0.size())
172  {
173  Vmath::Vmul(ntot, &df[0][0], 1, &Diff0[0], 1, &out_d0[0], 1);
174  Vmath::Vvtvp(ntot, &df[1][0], 1, &Diff1[0], 1, &out_d0[0], 1,
175  &out_d0[0], 1);
176  Vmath::Vvtvp(ntot, &df[2][0], 1, &Diff2[0], 1, &out_d0[0], 1,
177  &out_d0[0], 1);
178  }
179 
180  if (out_d1.size())
181  {
182  Vmath::Vmul(ntot, &df[3][0], 1, &Diff0[0], 1, &out_d1[0], 1);
183  Vmath::Vvtvp(ntot, &df[4][0], 1, &Diff1[0], 1, &out_d1[0], 1,
184  &out_d1[0], 1);
185  Vmath::Vvtvp(ntot, &df[5][0], 1, &Diff2[0], 1, &out_d1[0], 1,
186  &out_d1[0], 1);
187  }
188 
189  if (out_d2.size())
190  {
191  Vmath::Vmul(ntot, &df[6][0], 1, &Diff0[0], 1, &out_d2[0], 1);
192  Vmath::Vvtvp(ntot, &df[7][0], 1, &Diff1[0], 1, &out_d2[0], 1,
193  &out_d2[0], 1);
194  Vmath::Vvtvp(ntot, &df[8][0], 1, &Diff2[0], 1, &out_d2[0], 1,
195  &out_d2[0], 1);
196  }
197  }
198  else // regular geometry
199  {
200  if (out_d0.size())
201  {
202  Vmath::Smul(ntot, df[0][0], &Diff0[0], 1, &out_d0[0], 1);
203  Blas::Daxpy(ntot, df[1][0], &Diff1[0], 1, &out_d0[0], 1);
204  Blas::Daxpy(ntot, df[2][0], &Diff2[0], 1, &out_d0[0], 1);
205  }
206 
207  if (out_d1.size())
208  {
209  Vmath::Smul(ntot, df[3][0], &Diff0[0], 1, &out_d1[0], 1);
210  Blas::Daxpy(ntot, df[4][0], &Diff1[0], 1, &out_d1[0], 1);
211  Blas::Daxpy(ntot, df[5][0], &Diff2[0], 1, &out_d1[0], 1);
212  }
213 
214  if (out_d2.size())
215  {
216  Vmath::Smul(ntot, df[6][0], &Diff0[0], 1, &out_d2[0], 1);
217  Blas::Daxpy(ntot, df[7][0], &Diff1[0], 1, &out_d2[0], 1);
218  Blas::Daxpy(ntot, df[8][0], &Diff2[0], 1, &out_d2[0], 1);
219  }
220  }
221 }
222 
223 /**
224  * \brief Calculate the derivative of the physical points in a single
225  * direction.
226  *
227  * @param dir Direction in which to compute derivative.
228  * Valid values are 0, 1, 2.
229  * @param inarray Input array.
230  * @param outarray Output array.
231  */
232 void HexExp::v_PhysDeriv(const int dir,
233  const Array<OneD, const NekDouble> &inarray,
234  Array<OneD, NekDouble> &outarray)
235 {
236  switch (dir)
237  {
238  case 0:
239  {
240  PhysDeriv(inarray, outarray, NullNekDouble1DArray,
242  }
243  break;
244  case 1:
245  {
246  PhysDeriv(inarray, NullNekDouble1DArray, outarray,
248  }
249  break;
250  case 2:
251  {
253  outarray);
254  }
255  break;
256  default:
257  {
258  ASSERTL1(false, "input dir is out of range");
259  }
260  break;
261  }
262 }
263 
265  const Array<OneD, const NekDouble> &inarray,
266  const Array<OneD, const NekDouble> &direction,
267  Array<OneD, NekDouble> &outarray)
268 {
269 
270  int shapedim = 3;
271  int nquad0 = m_base[0]->GetNumPoints();
272  int nquad1 = m_base[1]->GetNumPoints();
273  int nquad2 = m_base[2]->GetNumPoints();
274  int ntot = nquad0 * nquad1 * nquad2;
275 
277  m_metricinfo->GetDerivFactors(GetPointsKeys());
281 
282  StdHexExp::v_PhysDeriv(inarray, Diff0, Diff1, Diff2);
283 
284  Array<OneD, Array<OneD, NekDouble>> dfdir(shapedim);
285  Expansion::ComputeGmatcdotMF(df, direction, dfdir);
286 
287  Vmath::Vmul(ntot, &dfdir[0][0], 1, &Diff0[0], 1, &outarray[0], 1);
288  Vmath::Vvtvp(ntot, &dfdir[1][0], 1, &Diff1[0], 1, &outarray[0], 1,
289  &outarray[0], 1);
290  Vmath::Vvtvp(ntot, &dfdir[2][0], 1, &Diff2[0], 1, &outarray[0], 1,
291  &outarray[0], 1);
292 }
293 
294 //-----------------------------
295 // Transforms
296 //-----------------------------
297 
298 /**
299  * \brief Forward transform from physical quadrature space stored in \a
300  * inarray and evaluate the expansion coefficients and store in
301  * \a (this)->_coeffs
302  *
303  * @param inarray Input array
304  * @param outarray Output array
305  */
307  Array<OneD, NekDouble> &outarray)
308 {
309  if (m_base[0]->Collocation() && m_base[1]->Collocation() &&
310  m_base[2]->Collocation())
311  {
312  Vmath::Vcopy(GetNcoeffs(), &inarray[0], 1, &outarray[0], 1);
313  }
314  else
315  {
316  IProductWRTBase(inarray, outarray);
317 
318  // get Mass matrix inverse
319  MatrixKey masskey(StdRegions::eInvMass, DetShapeType(), *this);
320  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
321 
322  // copy inarray in case inarray == outarray
323  DNekVec in(m_ncoeffs, outarray);
324  DNekVec out(m_ncoeffs, outarray, eWrapper);
325 
326  out = (*matsys) * in;
327  }
328 }
329 
330 //-----------------------------
331 // Inner product functions
332 //-----------------------------
333 
334 /**
335  * \brief Calculate the inner product of inarray with respect to the
336  * elements basis.
337  *
338  * @param inarray Input array of physical space data.
339  * @param outarray Output array of data.
340  */
342  Array<OneD, NekDouble> &outarray)
343 {
344  HexExp::v_IProductWRTBase_SumFac(inarray, outarray);
345 }
346 
347 /**
348  * \brief Calculate the inner product of inarray with respect to the
349  * given basis B = base0 * base1 * base2.
350  *
351  * \f$ \begin{array}{rcl} I_{pqr} = (\phi_{pqr}, u)_{\delta}
352  * & = & \sum_{i=0}^{nq_0} \sum_{j=0}^{nq_1} \sum_{k=0}^{nq_2}
353  * \psi_{p}^{a} (\xi_{1i}) \psi_{q}^{a} (\xi_{2j}) \psi_{r}^{a}
354  * (\xi_{3k}) w_i w_j w_k u(\xi_{1,i} \xi_{2,j} \xi_{3,k})
355  * J_{i,j,k}\\ & = & \sum_{i=0}^{nq_0} \psi_p^a(\xi_{1,i})
356  * \sum_{j=0}^{nq_1} \psi_{q}^a(\xi_{2,j}) \sum_{k=0}^{nq_2}
357  * \psi_{r}^a u(\xi_{1i},\xi_{2j},\xi_{3k})
358  * J_{i,j,k} \end{array} \f$ \n
359  * where
360  * \f$ \phi_{pqr} (\xi_1 , \xi_2 , \xi_3)
361  * = \psi_p^a ( \xi_1) \psi_{q}^a (\xi_2) \psi_{r}^a (\xi_3) \f$ \n
362  * which can be implemented as \n
363  * \f$f_{r} (\xi_{3k})
364  * = \sum_{k=0}^{nq_3} \psi_{r}^a u(\xi_{1i},\xi_{2j},\xi_{3k})
365  * J_{i,j,k} = {\bf B_3 U} \f$ \n
366  * \f$ g_{q} (\xi_{3k}) = \sum_{j=0}^{nq_1} \psi_{q}^a (\xi_{2j})
367  * f_{r} (\xi_{3k}) = {\bf B_2 F} \f$ \n
368  * \f$ (\phi_{pqr}, u)_{\delta}
369  * = \sum_{k=0}^{nq_0} \psi_{p}^a (\xi_{3k}) g_{q} (\xi_{3k})
370  * = {\bf B_1 G} \f$
371  *
372  * @param base0 Basis to integrate wrt in first dimension.
373  * @param base1 Basis to integrate wrt in second dimension.
374  * @param base2 Basis to integrate wrt in third dimension.
375  * @param inarray Input array.
376  * @param outarray Output array.
377  * @param coll_check (not used)
378  */
380  const Array<OneD, const NekDouble> &inarray,
381  Array<OneD, NekDouble> &outarray, bool multiplybyweights)
382 {
383  int nquad0 = m_base[0]->GetNumPoints();
384  int nquad1 = m_base[1]->GetNumPoints();
385  int nquad2 = m_base[2]->GetNumPoints();
386  int order0 = m_base[0]->GetNumModes();
387  int order1 = m_base[1]->GetNumModes();
388 
389  Array<OneD, NekDouble> wsp(nquad0 * nquad1 * (nquad2 + order0) +
390  order0 * order1 * nquad2);
391 
392  if (multiplybyweights)
393  {
394  Array<OneD, NekDouble> tmp(inarray.size());
395 
396  MultiplyByQuadratureMetric(inarray, tmp);
398  m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
399  tmp, outarray, wsp, true, true, true);
400  }
401  else
402  {
404  m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
405  inarray, outarray, wsp, true, true, true);
406  }
407 }
408 
410  const Array<OneD, const NekDouble> &inarray,
411  Array<OneD, NekDouble> &outarray)
412 {
413  HexExp::IProductWRTDerivBase_SumFac(dir, inarray, outarray);
414 }
415 
416 /**
417  * @brief Calculates the inner product \f$ I_{pqr} = (u,
418  * \partial_{x_i} \phi_{pqr}) \f$.
419  *
420  * The derivative of the basis functions is performed using the chain
421  * rule in order to incorporate the geometric factors. Assuming that
422  * the basis functions are a tensor product
423  * \f$\phi_{pqr}(\xi_1,\xi_2,\xi_3) =
424  * \phi_1(\xi_1)\phi_2(\xi_2)\phi_3(\xi_3)\f$, in the hexahedral
425  * element, this is straightforward and yields the result
426  *
427  * \f[
428  * I_{pqr} = \sum_{k=1}^3 \left(u, \frac{\partial u}{\partial \xi_k}
429  * \frac{\partial \xi_k}{\partial x_i}\right)
430  * \f]
431  *
432  * @param dir Direction in which to take the derivative.
433  * @param inarray The function \f$ u \f$.
434  * @param outarray Value of the inner product.
435  */
437  const int dir, const Array<OneD, const NekDouble> &inarray,
438  Array<OneD, NekDouble> &outarray)
439 {
440  ASSERTL1((dir == 0) || (dir == 1) || (dir == 2), "Invalid direction.");
441 
442  const int nq0 = m_base[0]->GetNumPoints();
443  const int nq1 = m_base[1]->GetNumPoints();
444  const int nq2 = m_base[2]->GetNumPoints();
445  const int nq = nq0 * nq1 * nq2;
446  const int nm0 = m_base[0]->GetNumModes();
447  const int nm1 = m_base[1]->GetNumModes();
448 
449  Array<OneD, NekDouble> alloc(4 * nq + m_ncoeffs + nm0 * nq2 * (nq1 + nm1));
450  Array<OneD, NekDouble> tmp1(alloc); // Quad metric
451  Array<OneD, NekDouble> tmp2(alloc + nq); // Dir1 metric
452  Array<OneD, NekDouble> tmp3(alloc + 2 * nq); // Dir2 metric
453  Array<OneD, NekDouble> tmp4(alloc + 3 * nq); // Dir3 metric
454  Array<OneD, NekDouble> tmp5(alloc + 4 * nq); // iprod tmp
455  Array<OneD, NekDouble> wsp(tmp5 + m_ncoeffs); // Wsp
456 
457  MultiplyByQuadratureMetric(inarray, tmp1);
458 
460  tmp2D[0] = tmp2;
461  tmp2D[1] = tmp3;
462  tmp2D[2] = tmp4;
463 
464  HexExp::v_AlignVectorToCollapsedDir(dir, tmp1, tmp2D);
465 
466  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
467  m_base[2]->GetBdata(), tmp2, outarray, wsp,
468  false, true, true);
469 
470  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
471  m_base[2]->GetBdata(), tmp3, tmp5, wsp, true,
472  false, true);
473  Vmath::Vadd(m_ncoeffs, tmp5, 1, outarray, 1, outarray, 1);
474 
475  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(),
476  m_base[2]->GetDbdata(), tmp4, tmp5, wsp, true,
477  true, false);
478  Vmath::Vadd(m_ncoeffs, tmp5, 1, outarray, 1, outarray, 1);
479 }
480 
482  const int dir, const Array<OneD, const NekDouble> &inarray,
483  Array<OneD, Array<OneD, NekDouble>> &outarray)
484 {
485  ASSERTL1((dir == 0) || (dir == 1) || (dir == 2), "Invalid direction.");
486 
487  const int nq0 = m_base[0]->GetNumPoints();
488  const int nq1 = m_base[1]->GetNumPoints();
489  const int nq2 = m_base[2]->GetNumPoints();
490  const int nq = nq0 * nq1 * nq2;
491 
492  const Array<TwoD, const NekDouble> &df =
493  m_metricinfo->GetDerivFactors(GetPointsKeys());
494 
495  Array<OneD, NekDouble> tmp1(nq); // Quad metric
496 
497  Array<OneD, NekDouble> tmp2 = outarray[0]; // Dir1 metric
498  Array<OneD, NekDouble> tmp3 = outarray[1]; // Dir2 metric
499  Array<OneD, NekDouble> tmp4 = outarray[2];
500 
501  Vmath::Vcopy(nq, inarray, 1, tmp1, 1); // Dir3 metric
502 
503  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
504  {
505  Vmath::Vmul(nq, &df[3 * dir][0], 1, tmp1.get(), 1, tmp2.get(), 1);
506  Vmath::Vmul(nq, &df[3 * dir + 1][0], 1, tmp1.get(), 1, tmp3.get(), 1);
507  Vmath::Vmul(nq, &df[3 * dir + 2][0], 1, tmp1.get(), 1, tmp4.get(), 1);
508  }
509  else
510  {
511  Vmath::Smul(nq, df[3 * dir][0], tmp1.get(), 1, tmp2.get(), 1);
512  Vmath::Smul(nq, df[3 * dir + 1][0], tmp1.get(), 1, tmp3.get(), 1);
513  Vmath::Smul(nq, df[3 * dir + 2][0], tmp1.get(), 1, tmp4.get(), 1);
514  }
515 }
516 
518  const int dir, const Array<OneD, const NekDouble> &inarray,
519  Array<OneD, NekDouble> &outarray)
520 {
521  int nq = GetTotPoints();
523 
524  switch (dir)
525  {
526  case 0:
527  {
529  }
530  break;
531  case 1:
532  {
534  }
535  break;
536  case 2:
537  {
539  }
540  break;
541  default:
542  {
543  ASSERTL1(false, "input dir is out of range");
544  }
545  break;
546  }
547 
548  MatrixKey iprodmatkey(mtype, DetShapeType(), *this);
549  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
550 
551  Blas::Dgemv('N', m_ncoeffs, nq, iprodmat->Scale(),
552  (iprodmat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
553  inarray.get(), 1, 0.0, outarray.get(), 1);
554 }
555 
556 /**
557  *
558  * @param dir Vector direction in which to take the derivative.
559  * @param inarray The function \f$ u \f$.
560  * @param outarray Value of the inner product.
561  */
563  const Array<OneD, const NekDouble> &direction,
564  const Array<OneD, const NekDouble> &inarray,
565  Array<OneD, NekDouble> &outarray)
566 {
567  int shapedim = 3;
568  const int nq0 = m_base[0]->GetNumPoints();
569  const int nq1 = m_base[1]->GetNumPoints();
570  const int nq2 = m_base[2]->GetNumPoints();
571  const int nq = nq0 * nq1 * nq2;
572  const int nm0 = m_base[0]->GetNumModes();
573  const int nm1 = m_base[1]->GetNumModes();
574 
575  const Array<TwoD, const NekDouble> &df =
576  m_metricinfo->GetDerivFactors(GetPointsKeys());
577 
578  Array<OneD, NekDouble> alloc(4 * nq + m_ncoeffs + nm0 * nq2 * (nq1 + nm1));
579  Array<OneD, NekDouble> tmp1(alloc); // Quad metric
580  Array<OneD, NekDouble> tmp2(alloc + nq); // Dir1 metric
581  Array<OneD, NekDouble> tmp3(alloc + 2 * nq); // Dir2 metric
582  Array<OneD, NekDouble> tmp4(alloc + 3 * nq); // Dir3 metric
583  Array<OneD, NekDouble> tmp5(alloc + 4 * nq); // iprod tmp
584  Array<OneD, NekDouble> wsp(tmp5 + m_ncoeffs); // Wsp
585 
586  MultiplyByQuadratureMetric(inarray, tmp1);
587 
588  Array<OneD, Array<OneD, NekDouble>> dfdir(shapedim);
589  Expansion::ComputeGmatcdotMF(df, direction, dfdir);
590 
591  Vmath::Vmul(nq, &dfdir[0][0], 1, tmp1.get(), 1, tmp2.get(), 1);
592  Vmath::Vmul(nq, &dfdir[1][0], 1, tmp1.get(), 1, tmp3.get(), 1);
593  Vmath::Vmul(nq, &dfdir[2][0], 1, tmp1.get(), 1, tmp4.get(), 1);
594 
595  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
596  m_base[2]->GetBdata(), tmp2, outarray, wsp,
597  false, true, true);
598 
599  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
600  m_base[2]->GetBdata(), tmp3, tmp5, wsp, true,
601  false, true);
602 
603  Vmath::Vadd(m_ncoeffs, tmp5, 1, outarray, 1, outarray, 1);
604 
605  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(),
606  m_base[2]->GetDbdata(), tmp4, tmp5, wsp, true,
607  true, false);
608 
609  Vmath::Vadd(m_ncoeffs, tmp5, 1, outarray, 1, outarray, 1);
610 }
611 
612 //-----------------------------
613 // Evaluation functions
614 //-----------------------------
615 
616 /**
617  * Given the local cartesian coordinate \a Lcoord evaluate the
618  * value of physvals at this point by calling through to the
619  * StdExpansion method
620  */
622  const Array<OneD, const NekDouble> &Lcoord,
623  const Array<OneD, const NekDouble> &physvals)
624 {
625  // Evaluate point in local coordinates.
626  return StdHexExp::v_PhysEvaluate(Lcoord, physvals);
627 }
628 
630  const Array<OneD, const NekDouble> &physvals)
631 {
633 
634  ASSERTL0(m_geom, "m_geom not defined");
635  m_geom->GetLocCoords(coord, Lcoord);
636  return StdHexExp::v_PhysEvaluate(Lcoord, physvals);
637 }
638 
640 {
642  m_base[0]->GetBasisKey(), m_base[1]->GetBasisKey(),
643  m_base[2]->GetBasisKey());
644 }
645 
647 {
649  m_base[0]->GetPointsKey());
651  m_base[1]->GetPointsKey());
653  m_base[2]->GetPointsKey());
654 
656  bkey2);
657 }
658 
659 /**
660  * \brief Retrieves the physical coordinates of a given set of
661  * reference coordinates.
662  *
663  * @param Lcoords Local coordinates in reference space.
664  * @param coords Corresponding coordinates in physical space.
665  */
667  Array<OneD, NekDouble> &coords)
668 {
669  int i;
670 
671  ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[0] <= 1.0 && Lcoords[1] >= -1.0 &&
672  Lcoords[1] <= 1.0 && Lcoords[2] >= -1.0 && Lcoords[2] <= 1.0,
673  "Local coordinates are not in region [-1,1]");
674 
675  m_geom->FillGeom();
676 
677  for (i = 0; i < m_geom->GetCoordim(); ++i)
678  {
679  coords[i] = m_geom->GetCoord(i, Lcoords);
680  }
681 }
682 
684  Array<OneD, NekDouble> &coords_1,
685  Array<OneD, NekDouble> &coords_2)
686 {
687  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
688 }
689 
690 //-----------------------------
691 // Helper functions
692 //-----------------------------
693 
694 /// Return the region shape using the enum-list of ShapeType
696 {
698 }
699 
701  const NekDouble *data, const std::vector<unsigned int> &nummodes,
702  const int mode_offset, NekDouble *coeffs,
703  std::vector<LibUtilities::BasisType> &fromType)
704 {
705  int data_order0 = nummodes[mode_offset];
706  int fillorder0 = min(m_base[0]->GetNumModes(), data_order0);
707  int data_order1 = nummodes[mode_offset + 1];
708  int order1 = m_base[1]->GetNumModes();
709  int fillorder1 = min(order1, data_order1);
710  int data_order2 = nummodes[mode_offset + 2];
711  int order2 = m_base[2]->GetNumModes();
712  int fillorder2 = min(order2, data_order2);
713 
714  // Check if same basis
715  if (fromType[0] != m_base[0]->GetBasisType() ||
716  fromType[1] != m_base[1]->GetBasisType() ||
717  fromType[2] != m_base[2]->GetBasisType())
718  {
719  // Construct a hex with the appropriate basis type at our
720  // quadrature points, and one more to do a forwards
721  // transform. We can then copy the output to coeffs.
722  StdRegions::StdHexExp tmpHex(
723  LibUtilities::BasisKey(fromType[0], data_order0,
724  m_base[0]->GetPointsKey()),
725  LibUtilities::BasisKey(fromType[1], data_order1,
726  m_base[1]->GetPointsKey()),
727  LibUtilities::BasisKey(fromType[2], data_order2,
728  m_base[2]->GetPointsKey()));
729  StdRegions::StdHexExp tmpHex2(m_base[0]->GetBasisKey(),
730  m_base[1]->GetBasisKey(),
731  m_base[2]->GetBasisKey());
732 
733  Array<OneD, const NekDouble> tmpData(tmpHex.GetNcoeffs(), data);
734  Array<OneD, NekDouble> tmpBwd(tmpHex2.GetTotPoints());
735  Array<OneD, NekDouble> tmpOut(tmpHex2.GetNcoeffs());
736 
737  tmpHex.BwdTrans(tmpData, tmpBwd);
738  tmpHex2.FwdTrans(tmpBwd, tmpOut);
739  Vmath::Vcopy(tmpOut.size(), &tmpOut[0], 1, coeffs, 1);
740 
741  return;
742  }
743 
744  switch (m_base[0]->GetBasisType())
745  {
747  {
748  int i, j;
749  int cnt = 0;
750  int cnt1 = 0;
751 
753  "Extraction routine not set up for this basis");
755  "Extraction routine not set up for this basis");
756 
757  Vmath::Zero(m_ncoeffs, coeffs, 1);
758  for (j = 0; j < fillorder0; ++j)
759  {
760  for (i = 0; i < fillorder1; ++i)
761  {
762  Vmath::Vcopy(fillorder2, &data[cnt], 1, &coeffs[cnt1], 1);
763  cnt += data_order2;
764  cnt1 += order2;
765  }
766 
767  // count out data for j iteration
768  for (i = fillorder1; i < data_order1; ++i)
769  {
770  cnt += data_order2;
771  }
772 
773  for (i = fillorder1; i < order1; ++i)
774  {
775  cnt1 += order2;
776  }
777  }
778  break;
779  }
781  {
782  LibUtilities::PointsKey p0(nummodes[0],
784  LibUtilities::PointsKey p1(nummodes[1],
786  LibUtilities::PointsKey p2(nummodes[2],
788  LibUtilities::PointsKey t0(m_base[0]->GetNumModes(),
790  LibUtilities::PointsKey t1(m_base[1]->GetNumModes(),
792  LibUtilities::PointsKey t2(m_base[2]->GetNumModes(),
794  LibUtilities::Interp3D(p0, p1, p2, data, t0, t1, t2, coeffs);
795  }
796  break;
797  default:
798  ASSERTL0(false, "basis is either not set up or not "
799  "hierarchicial");
800  }
801 }
802 
803 bool HexExp::v_GetFaceDGForwards(const int i) const
804 {
805  StdRegions::Orientation fo = GetGeom3D()->GetForient(i);
806 
811 }
812 
813 void HexExp::v_GetTracePhysMap(const int face, Array<OneD, int> &outarray)
814 {
815  int nquad0 = m_base[0]->GetNumPoints();
816  int nquad1 = m_base[1]->GetNumPoints();
817  int nquad2 = m_base[2]->GetNumPoints();
818 
819  int nq0 = 0;
820  int nq1 = 0;
821 
822  switch (face)
823  {
824  case 0:
825  nq0 = nquad0;
826  nq1 = nquad1;
827 
828  // Directions A and B positive
829  if (outarray.size() != nq0 * nq1)
830  {
831  outarray = Array<OneD, int>(nq0 * nq1);
832  }
833 
834  for (int i = 0; i < nquad0 * nquad1; ++i)
835  {
836  outarray[i] = i;
837  }
838 
839  break;
840  case 1:
841  nq0 = nquad0;
842  nq1 = nquad2;
843  // Direction A and B positive
844  if (outarray.size() != nq0 * nq1)
845  {
846  outarray = Array<OneD, int>(nq0 * nq1);
847  }
848 
849  // Direction A and B positive
850  for (int k = 0; k < nquad2; k++)
851  {
852  for (int i = 0; i < nquad0; ++i)
853  {
854  outarray[k * nquad0 + i] = nquad0 * nquad1 * k + i;
855  }
856  }
857  break;
858  case 2:
859  nq0 = nquad1;
860  nq1 = nquad2;
861 
862  // Direction A and B positive
863  if (outarray.size() != nq0 * nq1)
864  {
865  outarray = Array<OneD, int>(nq0 * nq1);
866  }
867 
868  for (int i = 0; i < nquad1 * nquad2; i++)
869  {
870  outarray[i] = nquad0 - 1 + i * nquad0;
871  }
872  break;
873  case 3:
874  nq0 = nquad0;
875  nq1 = nquad2;
876 
877  // Direction A and B positive
878  if (outarray.size() != nq0 * nq1)
879  {
880  outarray = Array<OneD, int>(nq0 * nq1);
881  }
882 
883  for (int k = 0; k < nquad2; k++)
884  {
885  for (int i = 0; i < nquad0; i++)
886  {
887  outarray[k * nquad0 + i] =
888  (nquad0 * (nquad1 - 1)) + (k * nquad0 * nquad1) + i;
889  }
890  }
891  break;
892  case 4:
893  nq0 = nquad1;
894  nq1 = nquad2;
895 
896  // Direction A and B positive
897  if (outarray.size() != nq0 * nq1)
898  {
899  outarray = Array<OneD, int>(nq0 * nq1);
900  }
901 
902  for (int i = 0; i < nquad1 * nquad2; i++)
903  {
904  outarray[i] = i * nquad0;
905  }
906  break;
907  case 5:
908  nq0 = nquad0;
909  nq1 = nquad1;
910  // Directions A and B positive
911  if (outarray.size() != nq0 * nq1)
912  {
913  outarray = Array<OneD, int>(nq0 * nq1);
914  }
915 
916  for (int i = 0; i < nquad0 * nquad1; i++)
917  {
918  outarray[i] = nquad0 * nquad1 * (nquad2 - 1) + i;
919  }
920 
921  break;
922  default:
923  ASSERTL0(false, "face value (> 5) is out of range");
924  break;
925  }
926 }
927 
928 void HexExp::v_ComputeTraceNormal(const int face)
929 {
930  int i;
931  const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
932  GetGeom()->GetMetricInfo();
933  SpatialDomains::GeomType type = geomFactors->GetGtype();
934 
936  for (i = 0; i < ptsKeys.size(); ++i)
937  {
938  // Need at least 2 points for computing normals
939  if (ptsKeys[i].GetNumPoints() == 1)
940  {
941  LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
942  ptsKeys[i] = pKey;
943  }
944  }
945 
946  const Array<TwoD, const NekDouble> &df =
947  geomFactors->GetDerivFactors(ptsKeys);
948  const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
949 
950  LibUtilities::BasisKey tobasis0 = GetTraceBasisKey(face, 0);
951  LibUtilities::BasisKey tobasis1 = GetTraceBasisKey(face, 1);
952 
953  // Number of quadrature points in face expansion.
954  int nq_face = tobasis0.GetNumPoints() * tobasis1.GetNumPoints();
955 
956  int vCoordDim = GetCoordim();
957 
960  for (i = 0; i < vCoordDim; ++i)
961  {
962  normal[i] = Array<OneD, NekDouble>(nq_face);
963  }
964 
965  size_t nqb = nq_face;
966  size_t nbnd = face;
969 
970  // Regular geometry case
971  if ((type == SpatialDomains::eRegular) ||
973  {
974  NekDouble fac;
975  // Set up normals
976  switch (face)
977  {
978  case 0:
979  for (i = 0; i < vCoordDim; ++i)
980  {
981  normal[i][0] = -df[3 * i + 2][0];
982  }
983  break;
984  case 1:
985  for (i = 0; i < vCoordDim; ++i)
986  {
987  normal[i][0] = -df[3 * i + 1][0];
988  }
989  break;
990  case 2:
991  for (i = 0; i < vCoordDim; ++i)
992  {
993  normal[i][0] = df[3 * i][0];
994  }
995  break;
996  case 3:
997  for (i = 0; i < vCoordDim; ++i)
998  {
999  normal[i][0] = df[3 * i + 1][0];
1000  }
1001  break;
1002  case 4:
1003  for (i = 0; i < vCoordDim; ++i)
1004  {
1005  normal[i][0] = -df[3 * i][0];
1006  }
1007  break;
1008  case 5:
1009  for (i = 0; i < vCoordDim; ++i)
1010  {
1011  normal[i][0] = df[3 * i + 2][0];
1012  }
1013  break;
1014  default:
1015  ASSERTL0(false, "face is out of range (edge < 5)");
1016  }
1017 
1018  // normalise
1019  fac = 0.0;
1020  for (i = 0; i < vCoordDim; ++i)
1021  {
1022  fac += normal[i][0] * normal[i][0];
1023  }
1024  fac = 1.0 / sqrt(fac);
1025 
1026  Vmath::Fill(nqb, fac, length, 1);
1027  for (i = 0; i < vCoordDim; ++i)
1028  {
1029  Vmath::Fill(nq_face, fac * normal[i][0], normal[i], 1);
1030  }
1031  }
1032  else // Set up deformed normals
1033  {
1034  int j, k;
1035 
1036  int nqe0 = ptsKeys[0].GetNumPoints();
1037  int nqe1 = ptsKeys[0].GetNumPoints();
1038  int nqe2 = ptsKeys[0].GetNumPoints();
1039  int nqe01 = nqe0 * nqe1;
1040  int nqe02 = nqe0 * nqe2;
1041  int nqe12 = nqe1 * nqe2;
1042 
1043  int nqe;
1044  if (face == 0 || face == 5)
1045  {
1046  nqe = nqe01;
1047  }
1048  else if (face == 1 || face == 3)
1049  {
1050  nqe = nqe02;
1051  }
1052  else
1053  {
1054  nqe = nqe12;
1055  }
1056 
1057  LibUtilities::PointsKey points0;
1058  LibUtilities::PointsKey points1;
1059 
1060  Array<OneD, NekDouble> faceJac(nqe);
1061  Array<OneD, NekDouble> normals(vCoordDim * nqe, 0.0);
1062 
1063  // Extract Jacobian along face and recover local
1064  // derivates (dx/dr) for polynomial interpolation by
1065  // multiplying m_gmat by jacobian
1066  switch (face)
1067  {
1068  case 0:
1069  for (j = 0; j < nqe; ++j)
1070  {
1071  normals[j] = -df[2][j] * jac[j];
1072  normals[nqe + j] = -df[5][j] * jac[j];
1073  normals[2 * nqe + j] = -df[8][j] * jac[j];
1074  faceJac[j] = jac[j];
1075  }
1076 
1077  points0 = ptsKeys[0];
1078  points1 = ptsKeys[1];
1079  break;
1080  case 1:
1081  for (j = 0; j < nqe0; ++j)
1082  {
1083  for (k = 0; k < nqe2; ++k)
1084  {
1085  int idx = j + nqe01 * k;
1086  normals[j + k * nqe0] = -df[1][idx] * jac[idx];
1087  normals[nqe + j + k * nqe0] = -df[4][idx] * jac[idx];
1088  normals[2 * nqe + j + k * nqe0] =
1089  -df[7][idx] * jac[idx];
1090  faceJac[j + k * nqe0] = jac[idx];
1091  }
1092  }
1093  points0 = ptsKeys[0];
1094  points1 = ptsKeys[2];
1095  break;
1096  case 2:
1097  for (j = 0; j < nqe1; ++j)
1098  {
1099  for (k = 0; k < nqe2; ++k)
1100  {
1101  int idx = nqe0 - 1 + nqe0 * j + nqe01 * k;
1102  normals[j + k * nqe1] = df[0][idx] * jac[idx];
1103  normals[nqe + j + k * nqe1] = df[3][idx] * jac[idx];
1104  normals[2 * nqe + j + k * nqe1] = df[6][idx] * jac[idx];
1105  faceJac[j + k * nqe1] = jac[idx];
1106  }
1107  }
1108  points0 = ptsKeys[1];
1109  points1 = ptsKeys[2];
1110  break;
1111  case 3:
1112  for (j = 0; j < nqe0; ++j)
1113  {
1114  for (k = 0; k < nqe2; ++k)
1115  {
1116  int idx = nqe0 * (nqe1 - 1) + j + nqe01 * k;
1117  normals[j + k * nqe0] = df[1][idx] * jac[idx];
1118  normals[nqe + j + k * nqe0] = df[4][idx] * jac[idx];
1119  normals[2 * nqe + j + k * nqe0] = df[7][idx] * jac[idx];
1120  faceJac[j + k * nqe0] = jac[idx];
1121  }
1122  }
1123  points0 = ptsKeys[0];
1124  points1 = ptsKeys[2];
1125  break;
1126  case 4:
1127  for (j = 0; j < nqe1; ++j)
1128  {
1129  for (k = 0; k < nqe2; ++k)
1130  {
1131  int idx = j * nqe0 + nqe01 * k;
1132  normals[j + k * nqe1] = -df[0][idx] * jac[idx];
1133  normals[nqe + j + k * nqe1] = -df[3][idx] * jac[idx];
1134  normals[2 * nqe + j + k * nqe1] =
1135  -df[6][idx] * jac[idx];
1136  faceJac[j + k * nqe1] = jac[idx];
1137  }
1138  }
1139  points0 = ptsKeys[1];
1140  points1 = ptsKeys[2];
1141  break;
1142  case 5:
1143  for (j = 0; j < nqe01; ++j)
1144  {
1145  int idx = j + nqe01 * (nqe2 - 1);
1146  normals[j] = df[2][idx] * jac[idx];
1147  normals[nqe + j] = df[5][idx] * jac[idx];
1148  normals[2 * nqe + j] = df[8][idx] * jac[idx];
1149  faceJac[j] = jac[idx];
1150  }
1151  points0 = ptsKeys[0];
1152  points1 = ptsKeys[1];
1153  break;
1154  default:
1155  ASSERTL0(false, "face is out of range (face < 5)");
1156  }
1157 
1158  Array<OneD, NekDouble> work(nq_face, 0.0);
1159  // Interpolate Jacobian and invert
1160  LibUtilities::Interp2D(points0, points1, faceJac,
1161  tobasis0.GetPointsKey(), tobasis1.GetPointsKey(),
1162  work);
1163 
1164  Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
1165 
1166  // interpolate
1167  for (i = 0; i < GetCoordim(); ++i)
1168  {
1169  LibUtilities::Interp2D(points0, points1, &normals[i * nqe],
1170  tobasis0.GetPointsKey(),
1171  tobasis1.GetPointsKey(), &normal[i][0]);
1172  Vmath::Vmul(nq_face, work, 1, normal[i], 1, normal[i], 1);
1173  }
1174 
1175  // normalise normal vectors
1176  Vmath::Zero(nq_face, work, 1);
1177  for (i = 0; i < GetCoordim(); ++i)
1178  {
1179  Vmath::Vvtvp(nq_face, normal[i], 1, normal[i], 1, work, 1, work, 1);
1180  }
1181 
1182  Vmath::Vsqrt(nq_face, work, 1, work, 1);
1183  Vmath::Sdiv(nq_face, 1.0, work, 1, work, 1);
1184 
1185  Vmath::Vcopy(nqb, work, 1, length, 1);
1186 
1187  for (i = 0; i < GetCoordim(); ++i)
1188  {
1189  Vmath::Vmul(nq_face, normal[i], 1, work, 1, normal[i], 1);
1190  }
1191  }
1192 }
1193 
1194 //-----------------------------
1195 // Operator creation functions
1196 //-----------------------------
1198  Array<OneD, NekDouble> &outarray,
1199  const StdRegions::StdMatrixKey &mkey)
1200 {
1201  StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
1202 }
1203 
1205  Array<OneD, NekDouble> &outarray,
1206  const StdRegions::StdMatrixKey &mkey)
1207 {
1208  HexExp::v_LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
1209 }
1210 
1211 void HexExp::v_LaplacianMatrixOp(const int k1, const int k2,
1212  const Array<OneD, const NekDouble> &inarray,
1213  Array<OneD, NekDouble> &outarray,
1214  const StdRegions::StdMatrixKey &mkey)
1215 {
1216  StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1217 }
1218 
1220  const Array<OneD, const NekDouble> &inarray,
1221  Array<OneD, NekDouble> &outarray,
1222  const StdRegions::StdMatrixKey &mkey)
1223 {
1224  StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
1225 }
1226 
1228  const Array<OneD, const NekDouble> &inarray,
1229  Array<OneD, NekDouble> &outarray, const StdRegions::StdMatrixKey &mkey)
1230 {
1231  StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray, outarray, mkey);
1232 }
1233 
1235  const Array<OneD, const NekDouble> &inarray,
1236  Array<OneD, NekDouble> &outarray, const StdRegions::StdMatrixKey &mkey)
1237 {
1238  StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray, outarray, mkey);
1239 }
1240 
1242  Array<OneD, NekDouble> &outarray,
1243  const StdRegions::StdMatrixKey &mkey)
1244 {
1245  HexExp::v_HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
1246 }
1247 
1249  const Array<OneD, const NekDouble> &inarray,
1250  Array<OneD, NekDouble> &outarray, const StdRegions::StdMatrixKey &mkey)
1251 {
1252  // int nConsts = mkey.GetNconstants();
1253  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1254 
1255  // switch(nConsts)
1256  // {
1257  // case 0:
1258  // {
1259  // mat = GetLocMatrix(mkey.GetMatrixType());
1260  // }
1261  // break;
1262  // case 1:
1263  // {
1264  // mat =
1265  // GetLocMatrix(mkey.GetMatrixType(),mkey.GetConstant(0));
1266  // }
1267  // break;
1268  // case 2:
1269  // {
1270  // mat =
1271  // GetLocMatrix(mkey.GetMatrixType(),mkey.GetConstant(0),mkey.GetConstant(1));
1272  // }
1273  // break;
1274  //
1275  // default:
1276  // {
1277  // NEKERROR(ErrorUtil::efatal, "Unknown number of
1278  // constants");
1279  // }
1280  // break;
1281  // }
1282 
1283  if (inarray.get() == outarray.get())
1284  {
1286  Vmath::Vcopy(m_ncoeffs, inarray.get(), 1, tmp.get(), 1);
1287 
1288  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, mat->Scale(),
1289  (mat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
1290  tmp.get(), 1, 0.0, outarray.get(), 1);
1291  }
1292  else
1293  {
1294  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, mat->Scale(),
1295  (mat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
1296  inarray.get(), 1, 0.0, outarray.get(), 1);
1297  }
1298 }
1299 
1300 /**
1301  * This function is used to compute exactly the advective numerical flux
1302  * on the interface of two elements with different expansions, hence an
1303  * appropriate number of Gauss points has to be used. The number of
1304  * Gauss points has to be equal to the number used by the highest
1305  * polynomial degree of the two adjacent elements
1306  *
1307  * @param numMin Is the reduced polynomial order
1308  * @param inarray Input array of coefficients
1309  * @param dumpVar Output array of reduced coefficients.
1310  */
1312  const Array<OneD, const NekDouble> &inarray,
1313  Array<OneD, NekDouble> &outarray)
1314 {
1315  int n_coeffs = inarray.size();
1316  int nmodes0 = m_base[0]->GetNumModes();
1317  int nmodes1 = m_base[1]->GetNumModes();
1318  int nmodes2 = m_base[2]->GetNumModes();
1319  int numMax = nmodes0;
1320 
1321  Array<OneD, NekDouble> coeff(n_coeffs);
1322  Array<OneD, NekDouble> coeff_tmp1(nmodes0 * nmodes1, 0.0);
1323  Array<OneD, NekDouble> coeff_tmp2(n_coeffs, 0.0);
1324  Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
1325 
1326  Vmath::Vcopy(n_coeffs, inarray, 1, coeff_tmp2, 1);
1327 
1328  const LibUtilities::PointsKey Pkey0(nmodes0,
1330  const LibUtilities::PointsKey Pkey1(nmodes1,
1332  const LibUtilities::PointsKey Pkey2(nmodes2,
1334 
1335  LibUtilities::BasisKey b0(m_base[0]->GetBasisType(), nmodes0, Pkey0);
1336  LibUtilities::BasisKey b1(m_base[1]->GetBasisType(), nmodes1, Pkey1);
1337  LibUtilities::BasisKey b2(m_base[2]->GetBasisType(), nmodes2, Pkey2);
1338  LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A, nmodes0, Pkey0);
1339  LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_A, nmodes1, Pkey1);
1340  LibUtilities::BasisKey bortho2(LibUtilities::eOrtho_A, nmodes2, Pkey2);
1341 
1342  LibUtilities::InterpCoeff3D(b0, b1, b2, coeff_tmp2, bortho0, bortho1,
1343  bortho2, coeff);
1344 
1345  Vmath::Zero(n_coeffs, coeff_tmp2, 1);
1346 
1347  int cnt = 0, cnt2 = 0;
1348 
1349  for (int u = 0; u < numMin + 1; ++u)
1350  {
1351  for (int i = 0; i < numMin; ++i)
1352  {
1353  Vmath::Vcopy(numMin, tmp = coeff + cnt + cnt2, 1,
1354  tmp2 = coeff_tmp1 + cnt, 1);
1355 
1356  cnt = i * numMax;
1357  }
1358 
1359  Vmath::Vcopy(nmodes0 * nmodes1, tmp3 = coeff_tmp1, 1,
1360  tmp4 = coeff_tmp2 + cnt2, 1);
1361 
1362  cnt2 = u * nmodes0 * nmodes1;
1363  }
1364 
1365  LibUtilities::InterpCoeff3D(bortho0, bortho1, bortho2, coeff_tmp2, b0, b1,
1366  b2, outarray);
1367 }
1368 
1370  const StdRegions::StdMatrixKey &mkey)
1371 {
1372  int nq = GetTotPoints();
1373 
1374  // Calculate sqrt of the Jacobian
1376  Array<OneD, NekDouble> sqrt_jac(nq);
1377  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1378  {
1379  Vmath::Vsqrt(nq, jac, 1, sqrt_jac, 1);
1380  }
1381  else
1382  {
1383  Vmath::Fill(nq, sqrt(jac[0]), sqrt_jac, 1);
1384  }
1385 
1386  // Multiply array by sqrt(Jac)
1387  Vmath::Vmul(nq, sqrt_jac, 1, array, 1, array, 1);
1388 
1389  // Apply std region filter
1390  StdHexExp::v_SVVLaplacianFilter(array, mkey);
1391 
1392  // Divide by sqrt(Jac)
1393  Vmath::Vdiv(nq, array, 1, sqrt_jac, 1, array, 1);
1394 }
1395 
1396 //-----------------------------
1397 // Matrix creation functions
1398 //-----------------------------
1400 {
1401  DNekMatSharedPtr returnval;
1402 
1403  switch (mkey.GetMatrixType())
1404  {
1412  returnval = Expansion3D::v_GenMatrix(mkey);
1413  break;
1414  default:
1415  returnval = StdHexExp::v_GenMatrix(mkey);
1416  }
1417 
1418  return returnval;
1419 }
1420 
1422 {
1423  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1424  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1425  LibUtilities::BasisKey bkey2 = m_base[2]->GetBasisKey();
1426 
1428  MemoryManager<StdHexExp>::AllocateSharedPtr(bkey0, bkey1, bkey2);
1429 
1430  return tmp->GetStdMatrix(mkey);
1431 }
1432 
1434 {
1435  return m_matrixManager[mkey];
1436 }
1437 
1439 {
1440  return m_staticCondMatrixManager[mkey];
1441 }
1442 
1444 {
1445  m_staticCondMatrixManager.DeleteObject(mkey);
1446 }
1447 
1449  const Array<OneD, const NekDouble> &inarray,
1451 {
1452  // This implementation is only valid when there are no
1453  // coefficients associated to the Laplacian operator
1454  if (m_metrics.count(eMetricLaplacian00) == 0)
1455  {
1457  }
1458 
1459  int nquad0 = m_base[0]->GetNumPoints();
1460  int nquad1 = m_base[1]->GetNumPoints();
1461  int nquad2 = m_base[2]->GetNumPoints();
1462  int nqtot = nquad0 * nquad1 * nquad2;
1463 
1464  ASSERTL1(wsp.size() >= 6 * nqtot, "Insufficient workspace size.");
1465 
1466  const Array<OneD, const NekDouble> &base0 = m_base[0]->GetBdata();
1467  const Array<OneD, const NekDouble> &base1 = m_base[1]->GetBdata();
1468  const Array<OneD, const NekDouble> &base2 = m_base[2]->GetBdata();
1469  const Array<OneD, const NekDouble> &dbase0 = m_base[0]->GetDbdata();
1470  const Array<OneD, const NekDouble> &dbase1 = m_base[1]->GetDbdata();
1471  const Array<OneD, const NekDouble> &dbase2 = m_base[2]->GetDbdata();
1472  const Array<OneD, const NekDouble> &metric00 =
1474  const Array<OneD, const NekDouble> &metric01 =
1476  const Array<OneD, const NekDouble> &metric02 =
1478  const Array<OneD, const NekDouble> &metric11 =
1480  const Array<OneD, const NekDouble> &metric12 =
1482  const Array<OneD, const NekDouble> &metric22 =
1484 
1485  // Allocate temporary storage
1486  Array<OneD, NekDouble> wsp0(wsp);
1487  Array<OneD, NekDouble> wsp1(wsp + 1 * nqtot);
1488  Array<OneD, NekDouble> wsp2(wsp + 2 * nqtot);
1489  Array<OneD, NekDouble> wsp3(wsp + 3 * nqtot);
1490  Array<OneD, NekDouble> wsp4(wsp + 4 * nqtot);
1491  Array<OneD, NekDouble> wsp5(wsp + 5 * nqtot);
1492 
1493  StdExpansion3D::PhysTensorDeriv(inarray, wsp0, wsp1, wsp2);
1494 
1495  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1496  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1497  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1498  // especially for this purpose
1499  Vmath::Vvtvvtp(nqtot, &metric00[0], 1, &wsp0[0], 1, &metric01[0], 1,
1500  &wsp1[0], 1, &wsp3[0], 1);
1501  Vmath::Vvtvp(nqtot, &metric02[0], 1, &wsp2[0], 1, &wsp3[0], 1, &wsp3[0], 1);
1502  Vmath::Vvtvvtp(nqtot, &metric01[0], 1, &wsp0[0], 1, &metric11[0], 1,
1503  &wsp1[0], 1, &wsp4[0], 1);
1504  Vmath::Vvtvp(nqtot, &metric12[0], 1, &wsp2[0], 1, &wsp4[0], 1, &wsp4[0], 1);
1505  Vmath::Vvtvvtp(nqtot, &metric02[0], 1, &wsp0[0], 1, &metric12[0], 1,
1506  &wsp1[0], 1, &wsp5[0], 1);
1507  Vmath::Vvtvp(nqtot, &metric22[0], 1, &wsp2[0], 1, &wsp5[0], 1, &wsp5[0], 1);
1508 
1509  // outarray = m = (D_xi1 * B)^T * k
1510  // wsp1 = n = (D_xi2 * B)^T * l
1511  IProductWRTBase_SumFacKernel(dbase0, base1, base2, wsp3, outarray, wsp0,
1512  false, true, true);
1513  IProductWRTBase_SumFacKernel(base0, dbase1, base2, wsp4, wsp2, wsp0, true,
1514  false, true);
1515  Vmath::Vadd(m_ncoeffs, wsp2.get(), 1, outarray.get(), 1, outarray.get(), 1);
1516  IProductWRTBase_SumFacKernel(base0, base1, dbase2, wsp5, wsp2, wsp0, true,
1517  true, false);
1518  Vmath::Vadd(m_ncoeffs, wsp2.get(), 1, outarray.get(), 1, outarray.get(), 1);
1519 }
1520 
1522 {
1523  if (m_metrics.count(eMetricQuadrature) == 0)
1524  {
1526  }
1527 
1528  const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
1529  const unsigned int nqtot = GetTotPoints();
1530  const unsigned int dim = 3;
1531  const MetricType m[3][3] = {
1535 
1536  for (unsigned int i = 0; i < dim; ++i)
1537  {
1538  for (unsigned int j = i; j < dim; ++j)
1539  {
1540  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1541  const Array<TwoD, const NekDouble> &gmat =
1542  m_metricinfo->GetGmat(GetPointsKeys());
1543  if (type == SpatialDomains::eDeformed)
1544  {
1545  Vmath::Vcopy(nqtot, &gmat[i * dim + j][0], 1,
1546  &m_metrics[m[i][j]][0], 1);
1547  }
1548  else
1549  {
1550  Vmath::Fill(nqtot, gmat[i * dim + j][0], &m_metrics[m[i][j]][0],
1551  1);
1552  }
1553  MultiplyByQuadratureMetric(m_metrics[m[i][j]], m_metrics[m[i][j]]);
1554  }
1555  }
1556 }
1557 
1558 /** @brief: This method gets all of the factors which are
1559  required as part of the Gradient Jump Penalty
1560  stabilisation and involves the product of the normal and
1561  geometric factors along the element trace.
1562 */
1564  Array<OneD, Array<OneD, NekDouble>> &d0factors,
1565  Array<OneD, Array<OneD, NekDouble>> &d1factors,
1566  Array<OneD, Array<OneD, NekDouble>> &d2factors)
1567 {
1568  int nquad0 = GetNumPoints(0);
1569  int nquad1 = GetNumPoints(1);
1570  int nquad2 = GetNumPoints(2);
1571 
1572  const Array<TwoD, const NekDouble> &df =
1573  m_metricinfo->GetDerivFactors(GetPointsKeys());
1574 
1575  if (d0factors.size() != 6)
1576  {
1577  d0factors = Array<OneD, Array<OneD, NekDouble>>(6);
1578  d1factors = Array<OneD, Array<OneD, NekDouble>>(6);
1579  d2factors = Array<OneD, Array<OneD, NekDouble>>(6);
1580  }
1581 
1582  if (d0factors[0].size() != nquad0 * nquad1)
1583  {
1584  d0factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1585  d0factors[5] = Array<OneD, NekDouble>(nquad0 * nquad1);
1586  d1factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1587  d1factors[5] = Array<OneD, NekDouble>(nquad0 * nquad1);
1588  d2factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1589  d2factors[5] = Array<OneD, NekDouble>(nquad0 * nquad1);
1590  }
1591 
1592  if (d0factors[1].size() != nquad0 * nquad2)
1593  {
1594  d0factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1595  d0factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1596  d1factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1597  d1factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1598  d2factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1599  d2factors[3] = Array<OneD, NekDouble>(nquad0 * nquad2);
1600  }
1601 
1602  if (d0factors[2].size() != nquad1 * nquad2)
1603  {
1604  d0factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1605  d0factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1606  d1factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1607  d1factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1608  d2factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1609  d2factors[4] = Array<OneD, NekDouble>(nquad1 * nquad2);
1610  }
1611 
1612  // Outwards normals
1613  const Array<OneD, const Array<OneD, NekDouble>> &normal_0 =
1614  GetTraceNormal(0);
1615  const Array<OneD, const Array<OneD, NekDouble>> &normal_1 =
1616  GetTraceNormal(1);
1617  const Array<OneD, const Array<OneD, NekDouble>> &normal_2 =
1618  GetTraceNormal(2);
1619  const Array<OneD, const Array<OneD, NekDouble>> &normal_3 =
1620  GetTraceNormal(3);
1621  const Array<OneD, const Array<OneD, NekDouble>> &normal_4 =
1622  GetTraceNormal(4);
1623  const Array<OneD, const Array<OneD, NekDouble>> &normal_5 =
1624  GetTraceNormal(5);
1625 
1626  int ncoords = normal_0.size();
1627 
1628  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1629  {
1630  // faces 0 and 5
1631  for (int i = 0; i < nquad0 * nquad1; ++i)
1632  {
1633  d0factors[0][i] = df[0][i] * normal_0[0][i];
1634  d0factors[5][i] =
1635  df[0][nquad0 * nquad1 * (nquad2 - 1) + i] * normal_5[0][i];
1636  d1factors[0][i] = df[1][i] * normal_0[0][i];
1637  d1factors[5][i] =
1638  df[1][nquad0 * nquad1 * (nquad2 - 1) + i] * normal_5[0][i];
1639  d2factors[0][i] = df[2][i] * normal_0[0][i];
1640  d2factors[5][i] =
1641  df[2][nquad0 * nquad1 * (nquad2 - 1) + i] * normal_5[0][i];
1642  }
1643 
1644  for (int n = 1; n < ncoords; ++n)
1645  {
1646  for (int i = 0; i < nquad0 * nquad1; ++i)
1647  {
1648  d0factors[0][i] += df[3 * n][i] * normal_0[n][i];
1649  d0factors[5][i] +=
1650  df[3 * n][nquad0 * nquad1 * (nquad2 - 1) + i] *
1651  normal_5[n][i];
1652  d1factors[0][i] += df[3 * n + 1][i] * normal_0[n][i];
1653  d1factors[5][i] +=
1654  df[3 * n + 1][nquad0 * nquad1 * (nquad2 - 1) + i] *
1655  normal_5[n][i];
1656  d2factors[0][i] += df[3 * n + 2][i] * normal_0[n][i];
1657  d2factors[5][i] +=
1658  df[3 * n + 2][nquad0 * nquad1 * (nquad2 - 1) + i] *
1659  normal_5[n][i];
1660  }
1661  }
1662 
1663  // faces 1 and 3
1664  for (int j = 0; j < nquad2; ++j)
1665  {
1666  for (int i = 0; i < nquad0; ++i)
1667  {
1668  d0factors[1][i] = df[0][j * nquad0 * nquad1 + i] *
1669  normal_1[0][j * nquad0 + i];
1670  d0factors[3][i] =
1671  df[0][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1672  normal_3[0][j * nquad0 + i];
1673  d1factors[1][i] = df[1][j * nquad0 * nquad1 + i] *
1674  normal_1[0][j * nquad0 + i];
1675  d1factors[3][i] =
1676  df[1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1677  normal_3[0][j * nquad0 + i];
1678  d2factors[1][i] = df[2][j * nquad0 * nquad1 + i] *
1679  normal_1[0][j * nquad0 + i];
1680  d2factors[3][i] =
1681  df[2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1682  normal_3[0][j * nquad0 + i];
1683  }
1684  }
1685 
1686  for (int n = 1; n < ncoords; ++n)
1687  {
1688  for (int j = 0; j < nquad2; ++j)
1689  {
1690  for (int i = 0; i < nquad0; ++i)
1691  {
1692  d0factors[1][i] = df[3 * n][j * nquad0 * nquad1 + i] *
1693  normal_1[0][j * nquad0 + i];
1694  d0factors[3][i] =
1695  df[3 * n][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1696  normal_3[0][j * nquad0 + i];
1697  d1factors[1][i] = df[3 * n + 1][j * nquad0 * nquad1 + i] *
1698  normal_1[0][j * nquad0 + i];
1699  d1factors[3][i] =
1700  df[3 * n + 1][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1701  normal_3[0][j * nquad0 + i];
1702  d2factors[1][i] = df[3 * n + 2][j * nquad0 * nquad1 + i] *
1703  normal_1[0][j * nquad0 + i];
1704  d2factors[3][i] =
1705  df[3 * n + 2][(j + 1) * nquad0 * nquad1 - nquad0 + i] *
1706  normal_3[0][j * nquad0 + i];
1707  }
1708  }
1709  }
1710 
1711  // faces 2 and 4
1712  for (int j = 0; j < nquad2; ++j)
1713  {
1714  for (int i = 0; i < nquad1; ++i)
1715  {
1716  d0factors[2][j * nquad1 + i] =
1717  df[0][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1718  normal_2[0][j * nquad1 + i];
1719  d0factors[4][j * nquad1 + i] =
1720  df[0][j * nquad0 * nquad1 + i * nquad0] *
1721  normal_4[0][j * nquad1 + i];
1722  d1factors[2][j * nquad1 + i] =
1723  df[1][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1724  normal_2[0][j * nquad1 + i];
1725  d1factors[4][j * nquad1 + i] =
1726  df[1][j * nquad0 * nquad1 + i * nquad0] *
1727  normal_4[0][j * nquad1 + i];
1728  d2factors[2][j * nquad1 + i] =
1729  df[2][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1730  normal_2[0][j * nquad1 + i];
1731  d2factors[4][j * nquad1 + i] =
1732  df[2][j * nquad0 * nquad1 + i * nquad0] *
1733  normal_4[0][j * nquad1 + i];
1734  }
1735  }
1736 
1737  for (int n = 1; n < ncoords; ++n)
1738  {
1739  for (int j = 0; j < nquad2; ++j)
1740  {
1741  for (int i = 0; i < nquad1; ++i)
1742  {
1743  d0factors[2][j * nquad1 + i] +=
1744  df[3 * n][j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1745  normal_2[n][j * nquad0 + i];
1746  d0factors[4][j * nquad0 + i] +=
1747  df[3 * n][i * nquad0 + j * nquad0 * nquad1] *
1748  normal_4[n][j * nquad0 + i];
1749  d1factors[2][j * nquad1 + i] +=
1750  df[3 * n + 1]
1751  [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1752  normal_2[n][j * nquad0 + i];
1753  d1factors[4][j * nquad0 + i] +=
1754  df[3 * n + 1][i * nquad0 + j * nquad0 * nquad1] *
1755  normal_4[n][j * nquad0 + i];
1756  d2factors[2][j * nquad1 + i] +=
1757  df[3 * n + 2]
1758  [j * nquad0 * nquad1 + (i + 1) * nquad0 - 1] *
1759  normal_2[n][j * nquad0 + i];
1760  d2factors[4][j * nquad0 + i] +=
1761  df[3 * n + 2][i * nquad0 + j * nquad0 * nquad1] *
1762  normal_4[n][j * nquad0 + i];
1763  }
1764  }
1765  }
1766  }
1767  else
1768  {
1769  // Faces 0 and 5
1770  for (int i = 0; i < nquad0 * nquad1; ++i)
1771  {
1772  d0factors[0][i] = df[0][0] * normal_0[0][i];
1773  d0factors[5][i] = df[0][0] * normal_5[0][i];
1774 
1775  d1factors[0][i] = df[1][0] * normal_0[0][i];
1776  d1factors[5][i] = df[1][0] * normal_5[0][i];
1777 
1778  d2factors[0][i] = df[2][0] * normal_0[0][i];
1779  d2factors[5][i] = df[2][0] * normal_5[0][i];
1780  }
1781 
1782  for (int n = 1; n < ncoords; ++n)
1783  {
1784  for (int i = 0; i < nquad0 * nquad1; ++i)
1785  {
1786  d0factors[0][i] += df[3 * n][0] * normal_0[n][i];
1787  d0factors[5][i] += df[3 * n][0] * normal_5[n][i];
1788 
1789  d1factors[0][i] += df[3 * n + 1][0] * normal_0[n][i];
1790  d1factors[5][i] += df[3 * n + 1][0] * normal_5[n][i];
1791 
1792  d2factors[0][i] += df[3 * n + 2][0] * normal_0[n][i];
1793  d2factors[5][i] += df[3 * n + 2][0] * normal_5[n][i];
1794  }
1795  }
1796 
1797  // faces 1 and 3
1798  for (int i = 0; i < nquad0 * nquad2; ++i)
1799  {
1800  d0factors[1][i] = df[0][0] * normal_1[0][i];
1801  d0factors[3][i] = df[0][0] * normal_3[0][i];
1802 
1803  d1factors[1][i] = df[1][0] * normal_1[0][i];
1804  d1factors[3][i] = df[1][0] * normal_3[0][i];
1805 
1806  d2factors[1][i] = df[2][0] * normal_1[0][i];
1807  d2factors[3][i] = df[2][0] * normal_3[0][i];
1808  }
1809 
1810  for (int n = 1; n < ncoords; ++n)
1811  {
1812  for (int i = 0; i < nquad0 * nquad2; ++i)
1813  {
1814  d0factors[1][i] += df[3 * n][0] * normal_1[n][i];
1815  d0factors[3][i] += df[3 * n][0] * normal_3[n][i];
1816 
1817  d1factors[1][i] += df[3 * n + 1][0] * normal_1[n][i];
1818  d1factors[3][i] += df[3 * n + 1][0] * normal_3[n][i];
1819 
1820  d2factors[1][i] += df[3 * n + 2][0] * normal_1[n][i];
1821  d2factors[3][i] += df[3 * n + 2][0] * normal_3[n][i];
1822  }
1823  }
1824 
1825  // faces 2 and 4
1826  for (int i = 0; i < nquad1 * nquad2; ++i)
1827  {
1828  d0factors[2][i] = df[0][0] * normal_2[0][i];
1829  d0factors[4][i] = df[0][0] * normal_4[0][i];
1830 
1831  d1factors[2][i] = df[1][0] * normal_2[0][i];
1832  d1factors[4][i] = df[1][0] * normal_4[0][i];
1833 
1834  d2factors[2][i] = df[2][0] * normal_2[0][i];
1835  d2factors[4][i] = df[2][0] * normal_4[0][i];
1836  }
1837 
1838  for (int n = 1; n < ncoords; ++n)
1839  {
1840  for (int i = 0; i < nquad1 * nquad2; ++i)
1841  {
1842  d0factors[2][i] += df[3 * n][0] * normal_2[n][i];
1843  d0factors[4][i] += df[3 * n][0] * normal_4[n][i];
1844 
1845  d1factors[2][i] += df[3 * n + 1][0] * normal_2[n][i];
1846  d1factors[4][i] += df[3 * n + 1][0] * normal_4[n][i];
1847 
1848  d2factors[2][i] += df[3 * n + 2][0] * normal_2[n][i];
1849  d2factors[4][i] += df[3 * n + 2][0] * normal_4[n][i];
1850  }
1851  }
1852  }
1853 }
1854 } // namespace LocalRegions
1855 } // 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)
SpatialDomains::Geometry3DSharedPtr GetGeom3D() const
Definition: Expansion3D.h:181
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
void ComputeGmatcdotMF(const Array< TwoD, const NekDouble > &df, const Array< OneD, const NekDouble > &direction, Array< OneD, Array< OneD, NekDouble >> &dfdir)
Definition: Expansion.cpp:597
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 IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: HexExp.cpp:517
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Calculates the inner product .
Definition: HexExp.cpp:436
virtual NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals)
This function evaluates the expansion at a single (arbitrary) point of the domain.
Definition: HexExp.cpp:629
virtual void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
Definition: HexExp.cpp:1448
virtual NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
Definition: HexExp.cpp:621
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: HexExp.h:249
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: HexExp.cpp:700
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1204
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1197
virtual void v_WeakDirectionalDerivMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1227
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: HexExp.h:251
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Calculate the inner product of inarray with respect to the given basis B = base0 * base1 * base2.
Definition: HexExp.cpp:379
virtual void v_GetTracePhysMap(const int face, Array< OneD, int > &outarray)
Definition: HexExp.cpp:813
virtual LibUtilities::ShapeType v_DetShapeType() const
Return the region shape using the enum-list of ShapeType.
Definition: HexExp.cpp:695
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: HexExp.cpp:1443
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: HexExp.cpp:409
void IProductWRTDirectionalDerivBase_SumFac(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: HexExp.cpp:562
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: HexExp.cpp:151
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1369
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: HexExp.cpp:306
virtual void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: HexExp.cpp:481
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: HexExp.cpp:1433
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1219
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1241
void v_ComputeTraceNormal(const int face)
Definition: HexExp.cpp:928
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: HexExp.cpp:683
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
Definition: HexExp.cpp:646
virtual void v_ComputeLaplacianMetric()
Definition: HexExp.cpp:1521
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
Definition: HexExp.cpp:639
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: HexExp.cpp:1438
void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1248
virtual void v_MassLevelCurvatureMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1234
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Calculate the inner product of inarray with respect to the elements basis.
Definition: HexExp.cpp:341
virtual void v_NormalTraceDerivFactors(Array< OneD, Array< OneD, NekDouble >> &factors, Array< OneD, Array< OneD, NekDouble >> &d0factors, Array< OneD, Array< OneD, NekDouble >> &d1factors)
: This method gets all of the factors which are required as part of the Gradient Jump Penalty stabili...
Definition: HexExp.cpp:1563
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1399
void v_PhysDirectionalDeriv(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &direction, Array< OneD, NekDouble > &out)
Physical derivative along a direction vector.
Definition: HexExp.cpp:264
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: HexExp.cpp:1421
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: HexExp.cpp:1311
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrate the physical point list inarray over region.
Definition: HexExp.cpp:110
virtual bool v_GetFaceDGForwards(const int i) const
Definition: HexExp.cpp:803
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
Retrieves the physical coordinates of a given set of reference coordinates.
Definition: HexExp.cpp:666
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
virtual void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
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 IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
this function calculates the inner product of a given function f with the different modes of the expa...
Definition: StdExpansion.h:536
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
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Definition: StdExpansion.h:850
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Class representing a hexehedral element in reference space.
Definition: StdHexExp.h:48
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 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 Interp3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
this function interpolates a 3D function evaluated at the quadrature points of the 3D basis,...
Definition: Interp.cpp:167
void InterpCoeff3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
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
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:53
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:44
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:58
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:50
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< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:86
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::shared_ptr< StdHexExp > StdHexExpSharedPtr
Definition: StdHexExp.h:254
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
static Array< OneD, NekDouble > NullNekDouble1DArray
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 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 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