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