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