Nektar++
StdTetExp.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: StdTetExp.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: Header field for tetrahedral routines built upon
32// StdExpansion3D
33//
34///////////////////////////////////////////////////////////////////////////////
35
37
38using namespace std;
39
40namespace Nektar::StdRegions
41{
43 const LibUtilities::BasisKey &Bb,
44 const LibUtilities::BasisKey &Bc)
45 : StdExpansion(LibUtilities::StdTetData::getNumberOfCoefficients(
46 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
47 3, Ba, Bb, Bc),
48 StdExpansion3D(LibUtilities::StdTetData::getNumberOfCoefficients(
49 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
50 Ba, Bb, Bc)
51{
52 ASSERTL0(Ba.GetNumModes() <= Bb.GetNumModes(),
53 "order in 'a' direction is higher than order "
54 "in 'b' direction");
55 ASSERTL0(Ba.GetNumModes() <= Bc.GetNumModes(),
56 "order in 'a' direction is higher than order "
57 "in 'c' direction");
58 ASSERTL0(Bb.GetNumModes() <= Bc.GetNumModes(),
59 "order in 'b' direction is higher than order "
60 "in 'c' direction");
61}
62
63//----------------------------
64// Differentiation Methods
65//----------------------------
66
67/**
68 * \brief Calculate the derivative of the physical points
69 *
70 * The derivative is evaluated at the nodal physical points.
71 * Derivatives with respect to the local Cartesian coordinates
72 *
73 * \f$\begin{Bmatrix} \frac {\partial} {\partial \xi_1} \\ \frac
74 * {\partial} {\partial \xi_2} \\ \frac {\partial} {\partial \xi_3}
75 * \end{Bmatrix} = \begin{Bmatrix} \frac 4 {(1-\eta_2)(1-\eta_3)}
76 * \frac \partial {\partial \eta_1} \ \ \frac {2(1+\eta_1)}
77 * {(1-\eta_2)(1-\eta_3)} \frac \partial {\partial \eta_1} + \frac 2
78 * {1-\eta_3} \frac \partial {\partial \eta_3} \\ \frac {2(1 +
79 * \eta_1)} {2(1 - \eta_2)(1-\eta_3)} \frac \partial {\partial \eta_1}
80 * + \frac {1 + \eta_2} {1 - \eta_3} \frac \partial {\partial \eta_2}
81 * + \frac \partial {\partial \eta_3} \end{Bmatrix}\f$
82 **/
84 Array<OneD, NekDouble> &out_dxi0,
85 Array<OneD, NekDouble> &out_dxi1,
86 Array<OneD, NekDouble> &out_dxi2)
87{
88 int Q0 = m_base[0]->GetNumPoints();
89 int Q1 = m_base[1]->GetNumPoints();
90 int Q2 = m_base[2]->GetNumPoints();
91 int Qtot = Q0 * Q1 * Q2;
92
93 // Compute the physical derivative
94 Array<OneD, NekDouble> out_dEta0(3 * Qtot, 0.0);
95 Array<OneD, NekDouble> out_dEta1 = out_dEta0 + Qtot;
96 Array<OneD, NekDouble> out_dEta2 = out_dEta1 + Qtot;
97
98 bool Do_2 = (out_dxi2.size() > 0) ? true : false;
99 bool Do_1 = (out_dxi1.size() > 0) ? true : false;
100
101 if (Do_2) // Need all local derivatives
102 {
103 PhysTensorDeriv(inarray, out_dEta0, out_dEta1, out_dEta2);
104 }
105 else if (Do_1) // Need 0 and 1 derivatives
106 {
107 PhysTensorDeriv(inarray, out_dEta0, out_dEta1, NullNekDouble1DArray);
108 }
109 else // Only need Eta0 derivaitve
110 {
111 PhysTensorDeriv(inarray, out_dEta0, NullNekDouble1DArray,
113 }
114
115 Array<OneD, const NekDouble> eta_0, eta_1, eta_2;
116 eta_0 = m_base[0]->GetZ();
117 eta_1 = m_base[1]->GetZ();
118 eta_2 = m_base[2]->GetZ();
119
120 // calculate 2.0/((1-eta_1)(1-eta_2)) Out_dEta0
121
122 NekDouble *dEta0 = &out_dEta0[0];
123 NekDouble fac;
124 for (int k = 0; k < Q2; ++k)
125 {
126 for (int j = 0; j < Q1; ++j, dEta0 += Q0)
127 {
128 Vmath::Smul(Q0, 2.0 / (1.0 - eta_1[j]), dEta0, 1, dEta0, 1);
129 }
130 fac = 1.0 / (1.0 - eta_2[k]);
131 Vmath::Smul(Q0 * Q1, fac, &out_dEta0[0] + k * Q0 * Q1, 1,
132 &out_dEta0[0] + k * Q0 * Q1, 1);
133 }
134
135 if (out_dxi0.size() > 0)
136 {
137 // out_dxi0 = 4.0/((1-eta_1)(1-eta_2)) Out_dEta0
138 Vmath::Smul(Qtot, 2.0, out_dEta0, 1, out_dxi0, 1);
139 }
140
141 if (Do_1 || Do_2)
142 {
143 Array<OneD, NekDouble> Fac0(Q0);
144 Vmath::Sadd(Q0, 1.0, eta_0, 1, Fac0, 1);
145
146 // calculate 2.0*(1+eta_0)/((1-eta_1)(1-eta_2)) Out_dEta0
147 for (int k = 0; k < Q1 * Q2; ++k)
148 {
149 Vmath::Vmul(Q0, &Fac0[0], 1, &out_dEta0[0] + k * Q0, 1,
150 &out_dEta0[0] + k * Q0, 1);
151 }
152 // calculate 2/(1.0-eta_2) out_dEta1
153 for (int k = 0; k < Q2; ++k)
154 {
155 Vmath::Smul(Q0 * Q1, 2.0 / (1.0 - eta_2[k]),
156 &out_dEta1[0] + k * Q0 * Q1, 1,
157 &out_dEta1[0] + k * Q0 * Q1, 1);
158 }
159
160 if (Do_1)
161 {
162 // calculate out_dxi1 = 2.0(1+eta_0)/((1-eta_1)(1-eta_2)) Out_dEta0
163 // + 2/(1.0-eta_2) out_dEta1
164 Vmath::Vadd(Qtot, out_dEta0, 1, out_dEta1, 1, out_dxi1, 1);
165 }
166
167 if (Do_2)
168 {
169 // calculate (1 + eta_1)/(1 -eta_2)*out_dEta1
170 NekDouble *dEta1 = &out_dEta1[0];
171 for (int k = 0; k < Q2; ++k)
172 {
173 for (int j = 0; j < Q1; ++j, dEta1 += Q0)
174 {
175 Vmath::Smul(Q0, (1.0 + eta_1[j]) / 2.0, dEta1, 1, dEta1, 1);
176 }
177 }
178
179 // calculate out_dxi2 =
180 // 2.0(1+eta_0)/((1-eta_1)(1-eta_2)) Out_dEta0 +
181 // (1 + eta_1)/(1 -eta_2)*out_dEta1 + out_dEta2
182 Vmath::Vadd(Qtot, out_dEta0, 1, out_dEta1, 1, out_dxi2, 1);
183 Vmath::Vadd(Qtot, out_dEta2, 1, out_dxi2, 1, out_dxi2, 1);
184 }
185 }
186}
187
188/**
189 * @param dir Direction in which to compute derivative.
190 * Valid values are 0, 1, 2.
191 * @param inarray Input array.
192 * @param outarray Output array.
193 */
194void StdTetExp::v_PhysDeriv(const int dir,
195 const Array<OneD, const NekDouble> &inarray,
196 Array<OneD, NekDouble> &outarray)
197{
198 switch (dir)
199 {
200 case 0:
201 {
202 v_PhysDeriv(inarray, outarray, NullNekDouble1DArray,
204 break;
205 }
206 case 1:
207 {
208 v_PhysDeriv(inarray, NullNekDouble1DArray, outarray,
210 break;
211 }
212 case 2:
213 {
215 outarray);
216 break;
217 }
218 default:
219 {
220 ASSERTL1(false, "input dir is out of range");
221 }
222 break;
223 }
224}
225
230{
231 StdTetExp::v_PhysDeriv(inarray, out_d0, out_d1, out_d2);
232}
233
234void StdTetExp::v_StdPhysDeriv(const int dir,
235 const Array<OneD, const NekDouble> &inarray,
236 Array<OneD, NekDouble> &outarray)
237{
238 StdTetExp::v_PhysDeriv(dir, inarray, outarray);
239}
240
241//---------------------------------------
242// Transforms
243//---------------------------------------
244
245/**
246 * @note 'r' (base[2]) runs fastest in this element
247 *
248 * \f$ u^{\delta} (\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{m(pqr)} \hat
249 * u_{pqr} \phi_{pqr} (\xi_{1i}, \xi_{2j}, \xi_{3k})\f$
250 *
251 * Backward transformation is three dimensional tensorial expansion
252 * \f$ u (\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{p=0}^{Q_x} \psi_p^a
253 * (\xi_{1i}) \lbrace { \sum_{q=0}^{Q_y} \psi_{pq}^b (\xi_{2j})
254 * \lbrace { \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{pqr}^c (\xi_{3k})
255 * \rbrace} \rbrace}. \f$ And sumfactorizing step of the form is as:\\
256 *
257 * \f$ f_{pq} (\xi_{3k}) = \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{pqr}^c
258 * (\xi_{3k}),\\
259 *
260 * g_{p} (\xi_{2j}, \xi_{3k}) = \sum_{r=0}^{Q_y} \psi_{pq}^b
261 * (\xi_{2j}) f_{pq} (\xi_{3k}),\\
262 *
263 * u(\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{p=0}^{Q_x} \psi_{p}^a
264 * (\xi_{1i}) g_{p} (\xi_{2j}, \xi_{3k}). \f$
265 */
267 Array<OneD, NekDouble> &outarray)
268{
271 "Basis[1] is not a general tensor type");
272
275 "Basis[2] is not a general tensor type");
276
277 if (m_base[0]->Collocation() && m_base[1]->Collocation() &&
278 m_base[2]->Collocation())
279 {
281 m_base[2]->GetNumPoints(),
282 inarray, 1, outarray, 1);
283 }
284 else
285 {
286 StdTetExp::v_BwdTrans_SumFac(inarray, outarray);
287 }
288}
289
290/**
291 * Sum-factorisation implementation of the BwdTrans operation.
292 */
294 Array<OneD, NekDouble> &outarray)
295{
296 int nquad1 = m_base[1]->GetNumPoints();
297 int nquad2 = m_base[2]->GetNumPoints();
298 int order0 = m_base[0]->GetNumModes();
299 int order1 = m_base[1]->GetNumModes();
300
301 Array<OneD, NekDouble> wsp(nquad2 * order0 * (2 * order1 - order0 + 1) / 2 +
302 nquad2 * nquad1 * order0);
303
304 BwdTrans_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(),
305 m_base[2]->GetBdata(), inarray, outarray, wsp, true,
306 true, true);
307}
308
309/**
310 * @param base0 x-dirn basis matrix
311 * @param base1 y-dirn basis matrix
312 * @param base2 z-dirn basis matrix
313 * @param inarray Input vector of modes.
314 * @param outarray Output vector of physical space data.
315 * @param wsp Workspace of size Q_x*P_z*(P_y+Q_y)
316 * @param doCheckCollDir0 Check for collocation of basis.
317 * @param doCheckCollDir1 Check for collocation of basis.
318 * @param doCheckCollDir2 Check for collocation of basis.
319 * @todo Account for some directions being collocated. See
320 * StdQuadExp as an example.
321 */
323 const Array<OneD, const NekDouble> &base0,
324 const Array<OneD, const NekDouble> &base1,
325 const Array<OneD, const NekDouble> &base2,
326 const Array<OneD, const NekDouble> &inarray,
328 [[maybe_unused]] bool doCheckCollDir0,
329 [[maybe_unused]] bool doCheckCollDir1,
330 [[maybe_unused]] bool doCheckCollDir2)
331{
332 int nquad0 = m_base[0]->GetNumPoints();
333 int nquad1 = m_base[1]->GetNumPoints();
334 int nquad2 = m_base[2]->GetNumPoints();
335
336 int order0 = m_base[0]->GetNumModes();
337 int order1 = m_base[1]->GetNumModes();
338 int order2 = m_base[2]->GetNumModes();
339
340 Array<OneD, NekDouble> tmp = wsp;
342 tmp + nquad2 * order0 * (2 * order1 - order0 + 1) / 2;
343
344 int i, j, mode, mode1, cnt;
345
346 // Perform summation over '2' direction
347 mode = mode1 = cnt = 0;
348 for (i = 0; i < order0; ++i)
349 {
350 for (j = 0; j < order1 - i; ++j, ++cnt)
351 {
352 Blas::Dgemv('N', nquad2, order2 - i - j, 1.0,
353 base2.data() + mode * nquad2, nquad2,
354 inarray.data() + mode1, 1, 0.0,
355 tmp.data() + cnt * nquad2, 1);
356 mode += order2 - i - j;
357 mode1 += order2 - i - j;
358 }
359 // increment mode in case order1!=order2
360 for (j = order1 - i; j < order2 - i; ++j)
361 {
362 mode += order2 - i - j;
363 }
364 }
365
366 // fix for modified basis by adding split of top singular
367 // vertex mode - currently (1+c)/2 x (1-b)/2 x (1-a)/2
368 // component is evaluated
370 {
371 // top singular vertex - (1+c)/2 x (1+b)/2 x (1-a)/2 component
372 Blas::Daxpy(nquad2, inarray[1], base2.data() + nquad2, 1,
373 &tmp[0] + nquad2, 1);
374
375 // top singular vertex - (1+c)/2 x (1-b)/2 x (1+a)/2 component
376 Blas::Daxpy(nquad2, inarray[1], base2.data() + nquad2, 1,
377 &tmp[0] + order1 * nquad2, 1);
378 }
379
380 // Perform summation over '1' direction
381 mode = 0;
382 for (i = 0; i < order0; ++i)
383 {
384 Blas::Dgemm('N', 'T', nquad1, nquad2, order1 - i, 1.0,
385 base1.data() + mode * nquad1, nquad1,
386 tmp.data() + mode * nquad2, nquad2, 0.0,
387 tmp1.data() + i * nquad1 * nquad2, nquad1);
388 mode += order1 - i;
389 }
390
391 // fix for modified basis by adding additional split of
392 // top and base singular vertex modes as well as singular
393 // edge
395 {
396 // use tmp to sort out singular vertices and
397 // singular edge components with (1+b)/2 (1+a)/2 form
398 for (i = 0; i < nquad2; ++i)
399 {
400 Blas::Daxpy(nquad1, tmp[nquad2 + i], base1.data() + nquad1, 1,
401 &tmp1[nquad1 * nquad2] + i * nquad1, 1);
402 }
403 }
404
405 // Perform summation over '0' direction
406 Blas::Dgemm('N', 'T', nquad0, nquad1 * nquad2, order0, 1.0, base0.data(),
407 nquad0, tmp1.data(), nquad1 * nquad2, 0.0, outarray.data(),
408 nquad0);
409}
410
411/**
412 * @param inarray array of physical quadrature points to be
413 * transformed.
414 * @param outarray updated array of expansion coefficients.
415 */
417 Array<OneD, NekDouble> &outarray)
418{ // int numMax = nmodes0;
419 v_IProductWRTBase(inarray, outarray);
420
421 // get Mass matrix inverse
422 StdMatrixKey masskey(eInvMass, DetShapeType(), *this);
423 DNekMatSharedPtr matsys = GetStdMatrix(masskey);
424
425 // copy inarray in case inarray == outarray
426 DNekVec in(m_ncoeffs, outarray);
427 DNekVec out(m_ncoeffs, outarray, eWrapper);
428
429 out = (*matsys) * in;
430}
431
432//---------------------------------------
433// Inner product functions
434//---------------------------------------
435
436/**
437 * \f$ \begin{array}{rcl} I_{pqr} = (\phi_{pqr}, u)_{\delta} & = &
438 * \sum_{i=0}^{nq_0} \sum_{j=0}^{nq_1} \sum_{k=0}^{nq_2} \psi_{p}^{a}
439 * (\eta_{1i}) \psi_{pq}^{b} (\eta_{2j}) \psi_{pqr}^{c} (\eta_{3k})
440 * w_i w_j w_k u(\eta_{1,i} \eta_{2,j} \eta_{3,k}) J_{i,j,k}\\ & = &
441 * \sum_{i=0}^{nq_0} \psi_p^a(\eta_{1,i}) \sum_{j=0}^{nq_1}
442 * \psi_{pq}^b(\eta_{2,j}) \sum_{k=0}^{nq_2} \psi_{pqr}^c
443 * u(\eta_{1i},\eta_{2j},\eta_{3k}) J_{i,j,k} \end{array} \f$ \n
444 *
445 * where
446 *
447 * \f$ \phi_{pqr} (\xi_1 , \xi_2 , \xi_3) = \psi_p^a (\eta_1)
448 * \psi_{pq}^b (\eta_2) \psi_{pqr}^c (\eta_3) \f$
449 *
450 * which can be implemented as \n \f$f_{pqr} (\xi_{3k}) =
451 * \sum_{k=0}^{nq_3} \psi_{pqr}^c u(\eta_{1i},\eta_{2j},\eta_{3k})
452 *
453 * J_{i,j,k} = {\bf B_3 U} \f$ \n
454 *
455 * \f$ g_{pq} (\xi_{3k}) = \sum_{j=0}^{nq_1} \psi_{pq}^b (\xi_{2j})
456 * f_{pqr} (\xi_{3k}) = {\bf B_2 F} \f$ \n
457 *
458 * \f$ (\phi_{pqr}, u)_{\delta} = \sum_{k=0}^{nq_0} \psi_{p}^a
459 * (\xi_{3k}) g_{pq} (\xi_{3k}) = {\bf B_1 G} \f$
460 *
461 * @param inarray Function evaluated at physical collocation
462 * points.
463 * @param outarray Inner product with respect to each basis
464 * function over the element.
465 */
467 Array<OneD, NekDouble> &outarray)
468{
471 "Basis[1] is not a general tensor type");
472
475 "Basis[2] is not a general tensor type");
476
477 if (m_base[0]->Collocation() && m_base[1]->Collocation())
478 {
479 MultiplyByQuadratureMetric(inarray, outarray);
480 }
481 else
482 {
483 StdTetExp::v_IProductWRTBase_SumFac(inarray, outarray);
484 }
485}
486
487/**
488 * @param inarray Function evaluated at physical collocation
489 * points.
490 * @param outarray Inner product with respect to each basis
491 * function over the element.
492 */
494 const Array<OneD, const NekDouble> &inarray,
495 Array<OneD, NekDouble> &outarray, bool multiplybyweights)
496{
497 int nquad0 = m_base[0]->GetNumPoints();
498 int nquad1 = m_base[1]->GetNumPoints();
499 int nquad2 = m_base[2]->GetNumPoints();
500 int order0 = m_base[0]->GetNumModes();
501 int order1 = m_base[1]->GetNumModes();
502
503 Array<OneD, NekDouble> wsp(nquad1 * nquad2 * order0 +
504 nquad2 * order0 * (2 * order1 - order0 + 1) / 2);
505
506 if (multiplybyweights)
507 {
508 Array<OneD, NekDouble> tmp(nquad0 * nquad1 * nquad2);
509 MultiplyByQuadratureMetric(inarray, tmp);
510
512 m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
513 tmp, outarray, wsp, true, true, true);
514 }
515 else
516 {
518 m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
519 inarray, outarray, wsp, true, true, true);
520 }
521}
522
524 const Array<OneD, const NekDouble> &base0,
525 const Array<OneD, const NekDouble> &base1,
526 const Array<OneD, const NekDouble> &base2,
527 const Array<OneD, const NekDouble> &inarray,
529 [[maybe_unused]] bool doCheckCollDir0,
530 [[maybe_unused]] bool doCheckCollDir1,
531 [[maybe_unused]] bool doCheckCollDir2)
532{
533 int nquad0 = m_base[0]->GetNumPoints();
534 int nquad1 = m_base[1]->GetNumPoints();
535 int nquad2 = m_base[2]->GetNumPoints();
536
537 int order0 = m_base[0]->GetNumModes();
538 int order1 = m_base[1]->GetNumModes();
539 int order2 = m_base[2]->GetNumModes();
540
541 Array<OneD, NekDouble> tmp1 = wsp;
542 Array<OneD, NekDouble> tmp2 = wsp + nquad1 * nquad2 * order0;
543
544 int i, j, mode, mode1, cnt;
545
546 // Inner product with respect to the '0' direction
547 Blas::Dgemm('T', 'N', nquad1 * nquad2, order0, nquad0, 1.0, inarray.data(),
548 nquad0, base0.data(), nquad0, 0.0, tmp1.data(),
549 nquad1 * nquad2);
550
551 // Inner product with respect to the '1' direction
552 for (mode = i = 0; i < order0; ++i)
553 {
554 Blas::Dgemm('T', 'N', nquad2, order1 - i, nquad1, 1.0,
555 tmp1.data() + i * nquad1 * nquad2, nquad1,
556 base1.data() + mode * nquad1, nquad1, 0.0,
557 tmp2.data() + mode * nquad2, nquad2);
558 mode += order1 - i;
559 }
560
561 // fix for modified basis for base singular vertex
563 {
564 // base singular vertex and singular edge (1+b)/2
565 //(1+a)/2 components (makes tmp[nquad2] entry into (1+b)/2)
566 Blas::Dgemv('T', nquad1, nquad2, 1.0, tmp1.data() + nquad1 * nquad2,
567 nquad1, base1.data() + nquad1, 1, 1.0, tmp2.data() + nquad2,
568 1);
569 }
570
571 // Inner product with respect to the '2' direction
572 mode = mode1 = cnt = 0;
573 for (i = 0; i < order0; ++i)
574 {
575 for (j = 0; j < order1 - i; ++j, ++cnt)
576 {
577 Blas::Dgemv('T', nquad2, order2 - i - j, 1.0,
578 base2.data() + mode * nquad2, nquad2,
579 tmp2.data() + cnt * nquad2, 1, 0.0,
580 outarray.data() + mode1, 1);
581 mode += order2 - i - j;
582 mode1 += order2 - i - j;
583 }
584 // increment mode in case order1!=order2
585 for (j = order1 - i; j < order2 - i; ++j)
586 {
587 mode += order2 - i - j;
588 }
589 }
590
591 // fix for modified basis for top singular vertex component
592 // Already have evaluated (1+c)/2 (1-b)/2 (1-a)/2
594 {
595 // add in (1+c)/2 (1+b)/2 component
596 outarray[1] +=
597 Blas::Ddot(nquad2, base2.data() + nquad2, 1, &tmp2[nquad2], 1);
598
599 // add in (1+c)/2 (1-b)/2 (1+a)/2 component
600 outarray[1] += Blas::Ddot(nquad2, base2.data() + nquad2, 1,
601 &tmp2[nquad2 * order1], 1);
602 }
603}
604
606 const int dir, const Array<OneD, const NekDouble> &inarray,
607 Array<OneD, NekDouble> &outarray)
608{
609 StdTetExp::v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
610}
611
612/**
613 * @param inarray Function evaluated at physical collocation
614 * points.
615 * @param outarray Inner product with respect to each basis
616 * function over the element.
617 */
619 const int dir, const Array<OneD, const NekDouble> &inarray,
620 Array<OneD, NekDouble> &outarray)
621{
622 int i;
623 int nquad0 = m_base[0]->GetNumPoints();
624 int nquad1 = m_base[1]->GetNumPoints();
625 int nquad2 = m_base[2]->GetNumPoints();
626 int nqtot = nquad0 * nquad1 * nquad2;
627 int nmodes0 = m_base[0]->GetNumModes();
628 int nmodes1 = m_base[1]->GetNumModes();
629 int wspsize = nquad0 + nquad1 + nquad2 + max(nqtot, m_ncoeffs) +
630 nquad1 * nquad2 * nmodes0 +
631 nquad2 * nmodes0 * (2 * nmodes1 - nmodes0 + 1) / 2;
632
633 Array<OneD, NekDouble> gfac0(wspsize);
634 Array<OneD, NekDouble> gfac1(gfac0 + nquad0);
635 Array<OneD, NekDouble> gfac2(gfac1 + nquad1);
636 Array<OneD, NekDouble> tmp0(gfac2 + nquad2);
637 Array<OneD, NekDouble> wsp(tmp0 + max(nqtot, m_ncoeffs));
638
639 const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
640 const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
641 const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
642
643 // set up geometric factor: (1+z0)/2
644 for (i = 0; i < nquad0; ++i)
645 {
646 gfac0[i] = 0.5 * (1 + z0[i]);
647 }
648
649 // set up geometric factor: 2/(1-z1)
650 for (i = 0; i < nquad1; ++i)
651 {
652 gfac1[i] = 2.0 / (1 - z1[i]);
653 }
654
655 // Set up geometric factor: 2/(1-z2)
656 for (i = 0; i < nquad2; ++i)
657 {
658 gfac2[i] = 2.0 / (1 - z2[i]);
659 }
660
661 // Derivative in first direction is always scaled as follows
662 for (i = 0; i < nquad1 * nquad2; ++i)
663 {
664 Vmath::Smul(nquad0, gfac1[i % nquad1], &inarray[0] + i * nquad0, 1,
665 &tmp0[0] + i * nquad0, 1);
666 }
667 for (i = 0; i < nquad2; ++i)
668 {
669 Vmath::Smul(nquad0 * nquad1, gfac2[i], &tmp0[0] + i * nquad0 * nquad1,
670 1, &tmp0[0] + i * nquad0 * nquad1, 1);
671 }
672
673 MultiplyByQuadratureMetric(tmp0, tmp0);
674
675 switch (dir)
676 {
677 case 0:
678 {
680 m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
681 m_base[2]->GetBdata(), tmp0, outarray, wsp, false, true, true);
682 }
683 break;
684 case 1:
685 {
687
688 for (i = 0; i < nquad1 * nquad2; ++i)
689 {
690 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp0[0] + i * nquad0, 1,
691 &tmp0[0] + i * nquad0, 1);
692 }
693
695 m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
696 m_base[2]->GetBdata(), tmp0, tmp3, wsp, false, true, true);
697
698 for (i = 0; i < nquad2; ++i)
699 {
700 Vmath::Smul(nquad0 * nquad1, gfac2[i],
701 &inarray[0] + i * nquad0 * nquad1, 1,
702 &tmp0[0] + i * nquad0 * nquad1, 1);
703 }
704 MultiplyByQuadratureMetric(tmp0, tmp0);
706 m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
707 m_base[2]->GetBdata(), tmp0, outarray, wsp, true, false, true);
708 Vmath::Vadd(m_ncoeffs, &tmp3[0], 1, &outarray[0], 1, &outarray[0],
709 1);
710 }
711 break;
712 case 2:
713 {
716 for (i = 0; i < nquad1; ++i)
717 {
718 gfac1[i] = (1 + z1[i]) / 2;
719 }
720
721 for (i = 0; i < nquad1 * nquad2; ++i)
722 {
723 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp0[0] + i * nquad0, 1,
724 &tmp0[0] + i * nquad0, 1);
725 }
727 m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
728 m_base[2]->GetBdata(), tmp0, tmp3, wsp, false, true, true);
729
730 for (i = 0; i < nquad2; ++i)
731 {
732 Vmath::Smul(nquad0 * nquad1, gfac2[i],
733 &inarray[0] + i * nquad0 * nquad1, 1,
734 &tmp0[0] + i * nquad0 * nquad1, 1);
735 }
736 for (i = 0; i < nquad1 * nquad2; ++i)
737 {
738 Vmath::Smul(nquad0, gfac1[i % nquad1], &tmp0[0] + i * nquad0, 1,
739 &tmp0[0] + i * nquad0, 1);
740 }
741 MultiplyByQuadratureMetric(tmp0, tmp0);
743 m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
744 m_base[2]->GetBdata(), tmp0, tmp4, wsp, true, false, true);
745
746 MultiplyByQuadratureMetric(inarray, tmp0);
748 m_base[0]->GetBdata(), m_base[1]->GetBdata(),
749 m_base[2]->GetDbdata(), tmp0, outarray, wsp, true, true, false);
750
751 Vmath::Vadd(m_ncoeffs, &tmp3[0], 1, &outarray[0], 1, &outarray[0],
752 1);
753 Vmath::Vadd(m_ncoeffs, &tmp4[0], 1, &outarray[0], 1, &outarray[0],
754 1);
755 }
756 break;
757 default:
758 {
759 ASSERTL1(false, "input dir is out of range");
760 }
761 break;
762 }
763}
764
765//---------------------------------------
766// Evaluation functions
767//---------------------------------------
768
771{
772 NekDouble d2 = 1.0 - xi[2];
773 NekDouble d12 = -xi[1] - xi[2];
774 if (fabs(d2) < NekConstants::kNekZeroTol)
775 {
776 if (d2 >= 0.)
777 {
779 }
780 else
781 {
783 }
784 }
785 if (fabs(d12) < NekConstants::kNekZeroTol)
786 {
787 if (d12 >= 0.)
788 {
790 }
791 else
792 {
794 }
795 }
796 eta[0] = 2.0 * (1.0 + xi[0]) / d12 - 1.0;
797 eta[1] = 2.0 * (1.0 + xi[1]) / d2 - 1.0;
798 eta[2] = xi[2];
799}
800
803{
804 xi[2] = eta[2];
805 xi[1] = (1.0 + eta[1]) * (1.0 - xi[2]) * 0.5 - 1.0;
806 xi[0] = (1.0 + eta[0]) * (-xi[1] - xi[2]) * 0.5 - 1.0;
807}
808
809void StdTetExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
810{
812 tmp[mode] = 1.0;
813 StdTetExp::v_BwdTrans(tmp, outarray);
814}
815
817 const Array<OneD, const NekDouble> &coords, int mode)
818{
820 LocCoordToLocCollapsed(coords, coll);
821
822 const int nm1 = m_base[1]->GetNumModes();
823 const int nm2 = m_base[2]->GetNumModes();
824
825 const int b = 2 * nm2 + 1;
826 const int mode0 = floor(0.5 * (b - sqrt(b * b - 8.0 * mode / nm1)));
827 const int tmp =
828 mode - nm1 * (mode0 * (nm2 - 1) + 1 - (mode0 - 2) * (mode0 - 1) / 2);
829 const int mode1 = tmp / (nm2 - mode0);
830 const int mode2 = tmp % (nm2 - mode0);
831
833 {
834 // Handle the collapsed vertices and edges in the modified
835 // basis.
836 if (mode == 1)
837 {
838 // Collapsed top vertex
839 return StdExpansion::BaryEvaluateBasis<2>(coll[2], 1);
840 }
841 else if (mode0 == 0 && mode2 == 1)
842 {
843 return StdExpansion::BaryEvaluateBasis<1>(coll[1], 0) *
844 StdExpansion::BaryEvaluateBasis<2>(coll[2], 1);
845 }
846 else if (mode0 == 1 && mode1 == 1 && mode2 == 0)
847 {
848 return StdExpansion::BaryEvaluateBasis<0>(coll[0], 0) *
849 StdExpansion::BaryEvaluateBasis<1>(coll[1], 1);
850 }
851 }
852
853 return StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
854 StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
855 StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
856}
857
859 const Array<OneD, NekDouble> &coord,
860 const Array<OneD, const NekDouble> &inarray,
861 std::array<NekDouble, 3> &firstOrderDerivs)
862{
863 // Collapse coordinates
864 Array<OneD, NekDouble> coll(3, 0.0);
865 LocCoordToLocCollapsed(coord, coll);
866
867 // If near singularity do the old interpolation matrix method
868 if ((1 - coll[1]) < 1e-5 || (1 - coll[2]) < 1e-5)
869 {
870 int totPoints = GetTotPoints();
871 Array<OneD, NekDouble> EphysDeriv0(totPoints), EphysDeriv1(totPoints),
872 EphysDeriv2(totPoints);
873 PhysDeriv(inarray, EphysDeriv0, EphysDeriv1, EphysDeriv2);
874
876 I[0] = GetBase()[0]->GetI(coll);
877 I[1] = GetBase()[1]->GetI(coll + 1);
878 I[2] = GetBase()[2]->GetI(coll + 2);
879
880 firstOrderDerivs[0] = PhysEvaluate(I, EphysDeriv0);
881 firstOrderDerivs[1] = PhysEvaluate(I, EphysDeriv1);
882 firstOrderDerivs[2] = PhysEvaluate(I, EphysDeriv2);
883 return PhysEvaluate(I, inarray);
884 }
885
886 std::array<NekDouble, 3> interDeriv;
887 NekDouble val = BaryTensorDeriv(coll, inarray, interDeriv);
888
889 // calculate 2.0/((1-eta_1)(1-eta_2)) * Out_dEta0
890 NekDouble temp = 2.0 / ((1 - coll[1]) * (1 - coll[2]));
891 interDeriv[0] *= temp;
892
893 // out_dxi0 = 4.0/((1-eta_1)(1-eta_2)) * Out_dEta0
894 firstOrderDerivs[0] = 2 * interDeriv[0];
895
896 // fac0 = 1 + eta_0
897 NekDouble fac0;
898 fac0 = 1 + coll[0];
899
900 // calculate 2.0*(1+eta_0)/((1-eta_1)(1-eta_2)) * Out_dEta0
901 interDeriv[0] *= fac0;
902
903 // calculate 2/(1.0-eta_2) * out_dEta1
904 fac0 = 2 / (1 - coll[2]);
905 interDeriv[1] *= fac0;
906
907 // calculate out_dxi1 = 2.0(1+eta_0)/((1-eta_1)(1-eta_2))
908 // * Out_dEta0 + 2/(1.0-eta_2) out_dEta1
909 firstOrderDerivs[1] = interDeriv[0] + interDeriv[1];
910
911 // calculate (1 + eta_1)/(1 -eta_2)*out_dEta1
912 fac0 = (1 + coll[1]) / 2;
913 interDeriv[1] *= fac0;
914
915 // calculate out_dxi2 =
916 // 2.0(1+eta_0)/((1-eta_1)(1-eta_2)) Out_dEta0 +
917 // (1 + eta_1)/(1 -eta_2)*out_dEta1 + out_dEta2
918 firstOrderDerivs[2] = interDeriv[0] + interDeriv[1] + interDeriv[2];
919
920 return val;
921}
922
923void StdTetExp::v_GetTraceNumModes(const int fid, int &numModes0,
924 int &numModes1,
925 [[maybe_unused]] Orientation faceOrient)
926{
927 int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
928 m_base[2]->GetNumModes()};
929 switch (fid)
930 {
931 case 0:
932 {
933 numModes0 = nummodes[0];
934 numModes1 = nummodes[1];
935 }
936 break;
937 case 1:
938 {
939 numModes0 = nummodes[0];
940 numModes1 = nummodes[2];
941 }
942 break;
943 case 2:
944 case 3:
945 {
946 numModes0 = nummodes[1];
947 numModes1 = nummodes[2];
948 }
949 break;
950 }
951}
952
953//---------------------------
954// Helper functions
955//---------------------------
956
958{
959 return 4;
960}
961
963{
964 return 6;
965}
966
968{
969 return 4;
970}
971
973{
975}
976
978{
981 "BasisType is not a boundary interior form");
984 "BasisType is not a boundary interior form");
987 "BasisType is not a boundary interior form");
988
989 int P = m_base[0]->GetNumModes();
990 int Q = m_base[1]->GetNumModes();
991 int R = m_base[2]->GetNumModes();
992
994}
995
997{
1000 "BasisType is not a boundary interior form");
1003 "BasisType is not a boundary interior form");
1006 "BasisType is not a boundary interior form");
1007
1008 int P = m_base[0]->GetNumModes() - 1;
1009 int Q = m_base[1]->GetNumModes() - 1;
1010 int R = m_base[2]->GetNumModes() - 1;
1011
1012 return (Q + 1) + P * (1 + 2 * Q - P) / 2 // base face
1013 + (R + 1) + P * (1 + 2 * R - P) / 2 // front face
1014 + 2 * (R + 1) + Q * (1 + 2 * R - Q); // back two faces
1015}
1016
1017int StdTetExp::v_GetTraceNcoeffs(const int i) const
1018{
1019 ASSERTL2((i >= 0) && (i <= 3), "face id is out of range");
1020 int nFaceCoeffs = 0;
1021 int nummodesA, nummodesB, P, Q;
1022 if (i == 0)
1023 {
1024 nummodesA = GetBasisNumModes(0);
1025 nummodesB = GetBasisNumModes(1);
1026 }
1027 else if ((i == 1) || (i == 2))
1028 {
1029 nummodesA = GetBasisNumModes(0);
1030 nummodesB = GetBasisNumModes(2);
1031 }
1032 else
1033 {
1034 nummodesA = GetBasisNumModes(1);
1035 nummodesB = GetBasisNumModes(2);
1036 }
1037 P = nummodesA - 1;
1038 Q = nummodesB - 1;
1039 nFaceCoeffs = Q + 1 + (P * (1 + 2 * Q - P)) / 2;
1040 return nFaceCoeffs;
1041}
1042
1044{
1045 ASSERTL2((i >= 0) && (i <= 3), "face id is out of range");
1046 int Pi = m_base[0]->GetNumModes() - 2;
1047 int Qi = m_base[1]->GetNumModes() - 2;
1048 int Ri = m_base[2]->GetNumModes() - 2;
1049
1050 if ((i == 0))
1051 {
1052 return Pi * (2 * Qi - Pi - 1) / 2;
1053 }
1054 else if ((i == 1))
1055 {
1056 return Pi * (2 * Ri - Pi - 1) / 2;
1057 }
1058 else
1059 {
1060 return Qi * (2 * Ri - Qi - 1) / 2;
1061 }
1062}
1063
1065{
1066 ASSERTL2(i >= 0 && i <= 3, "face id is out of range");
1067
1068 if (i == 0)
1069 {
1070 return m_base[0]->GetNumPoints() * m_base[1]->GetNumPoints();
1071 }
1072 else if (i == 1)
1073 {
1074 return m_base[0]->GetNumPoints() * m_base[2]->GetNumPoints();
1075 }
1076 else
1077 {
1078 return m_base[1]->GetNumPoints() * m_base[2]->GetNumPoints();
1079 }
1080}
1081
1082int StdTetExp::v_GetEdgeNcoeffs(const int i) const
1083{
1084 ASSERTL2((i >= 0) && (i <= 5), "edge id is out of range");
1085 int P = m_base[0]->GetNumModes();
1086 int Q = m_base[1]->GetNumModes();
1087 int R = m_base[2]->GetNumModes();
1088
1089 if (i == 0)
1090 {
1091 return P;
1092 }
1093 else if (i == 1 || i == 2)
1094 {
1095 return Q;
1096 }
1097 else
1098 {
1099 return R;
1100 }
1101}
1102
1104 const int j) const
1105{
1106 ASSERTL2(i >= 0 && i <= 3, "face id is out of range");
1107 ASSERTL2(j == 0 || j == 1, "face direction is out of range");
1108
1109 if (i == 0)
1110 {
1111 return m_base[j]->GetPointsKey();
1112 }
1113 else if (i == 1)
1114 {
1115 return m_base[2 * j]->GetPointsKey();
1116 }
1117 else
1118 {
1119 return m_base[j + 1]->GetPointsKey();
1120 }
1121}
1122
1124 const std::vector<unsigned int> &nummodes, int &modes_offset)
1125{
1127 nummodes[modes_offset], nummodes[modes_offset + 1],
1128 nummodes[modes_offset + 2]);
1129 modes_offset += 3;
1130
1131 return nmodes;
1132}
1133
1135 const int k,
1136 bool UseGLL) const
1137{
1138 ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
1139 ASSERTL2(k == 0 || k == 1, "face direction out of range");
1140
1141 int dir = k;
1142 switch (i)
1143 {
1144 case 0:
1145 dir = k;
1146 break;
1147 case 1:
1148 dir = 2 * k;
1149 break;
1150 case 2:
1151 case 3:
1152 dir = k + 1;
1153 break;
1154 }
1155
1157 m_base[dir]->GetNumPoints(),
1158 m_base[dir]->GetNumModes(), UseGLL);
1159}
1160
1164{
1165 Array<OneD, const NekDouble> eta_x = m_base[0]->GetZ();
1166 Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
1167 Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
1168 int Qx = GetNumPoints(0);
1169 int Qy = GetNumPoints(1);
1170 int Qz = GetNumPoints(2);
1171
1172 // Convert collapsed coordinates into cartesian coordinates: eta
1173 // --> xi
1174 for (int k = 0; k < Qz; ++k)
1175 {
1176 for (int j = 0; j < Qy; ++j)
1177 {
1178 for (int i = 0; i < Qx; ++i)
1179 {
1180 int s = i + Qx * (j + Qy * k);
1181 xi_x[s] =
1182 (eta_x[i] + 1.0) * (1.0 - eta_y[j]) * (1.0 - eta_z[k]) / 4 -
1183 1.0;
1184 xi_y[s] = (eta_y[j] + 1.0) * (1.0 - eta_z[k]) / 2 - 1.0;
1185 xi_z[s] = eta_z[k];
1186 }
1187 }
1188 }
1189}
1190
1192{
1193 return (m_base[0]->GetBasisType() == LibUtilities::eModified_A) &&
1194 (m_base[1]->GetBasisType() == LibUtilities::eModified_B) &&
1196}
1197
1198//--------------------------
1199// Mappings
1200//--------------------------
1201int StdTetExp::v_GetVertexMap(const int localVertexId, bool useCoeffPacking)
1202{
1206 "Mapping not defined for this type of basis");
1207
1208 int localDOF = 0;
1209 if (useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1210 {
1211 switch (localVertexId)
1212 {
1213 case 0:
1214 {
1215 localDOF = GetMode(0, 0, 0);
1216 break;
1217 }
1218 case 1:
1219 {
1220 localDOF = GetMode(0, 0, 1);
1221 break;
1222 }
1223 case 2:
1224 {
1225 localDOF = GetMode(0, 1, 0);
1226 break;
1227 }
1228 case 3:
1229 {
1230 localDOF = GetMode(1, 0, 0);
1231 break;
1232 }
1233 default:
1234 {
1235 ASSERTL0(false, "Vertex ID must be between 0 and 3");
1236 break;
1237 }
1238 }
1239 }
1240 else
1241 {
1242 switch (localVertexId)
1243 {
1244 case 0:
1245 {
1246 localDOF = GetMode(0, 0, 0);
1247 break;
1248 }
1249 case 1:
1250 {
1251 localDOF = GetMode(1, 0, 0);
1252 break;
1253 }
1254 case 2:
1255 {
1256 localDOF = GetMode(0, 1, 0);
1257 break;
1258 }
1259 case 3:
1260 {
1261 localDOF = GetMode(0, 0, 1);
1262 break;
1263 }
1264 default:
1265 {
1266 ASSERTL0(false, "Vertex ID must be between 0 and 3");
1267 break;
1268 }
1269 }
1270 }
1271
1272 return localDOF;
1273}
1274
1275/**
1276 * Maps interior modes of an edge to the elemental modes.
1277 */
1278
1279/**
1280 * List of all interior modes in the expansion.
1281 */
1283{
1286 "BasisType is not a boundary interior form");
1289 "BasisType is not a boundary interior form");
1292 "BasisType is not a boundary interior form");
1293
1294 int P = m_base[0]->GetNumModes();
1295 int Q = m_base[1]->GetNumModes();
1296 int R = m_base[2]->GetNumModes();
1297
1298 int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1299
1300 if (outarray.size() != nIntCoeffs)
1301 {
1302 outarray = Array<OneD, unsigned int>(nIntCoeffs);
1303 }
1304
1305 int idx = 0;
1306 for (int i = 2; i < P; ++i)
1307 {
1308 for (int j = 1; j < Q - i - 1; ++j)
1309 {
1310 for (int k = 1; k < R - i - j; ++k)
1311 {
1312 outarray[idx++] = GetMode(i, j, k);
1313 }
1314 }
1315 }
1316}
1317
1318/**
1319 * List of all boundary modes in the the expansion.
1320 */
1322{
1325 "BasisType is not a boundary interior form");
1328 "BasisType is not a boundary interior form");
1331 "BasisType is not a boundary interior form");
1332
1333 int P = m_base[0]->GetNumModes();
1334 int Q = m_base[1]->GetNumModes();
1335 int R = m_base[2]->GetNumModes();
1336
1337 int i, j, k;
1338 int idx = 0;
1339
1340 int nBnd = NumBndryCoeffs();
1341
1342 if (outarray.size() != nBnd)
1343 {
1344 outarray = Array<OneD, unsigned int>(nBnd);
1345 }
1346
1347 for (i = 0; i < P; ++i)
1348 {
1349 // First two Q-R planes are entirely boundary modes
1350 if (i < 2)
1351 {
1352 for (j = 0; j < Q - i; j++)
1353 {
1354 for (k = 0; k < R - i - j; ++k)
1355 {
1356 outarray[idx++] = GetMode(i, j, k);
1357 }
1358 }
1359 }
1360 // Remaining Q-R planes contain boundary modes on bottom and
1361 // left edge.
1362 else
1363 {
1364 for (k = 0; k < R - i; ++k)
1365 {
1366 outarray[idx++] = GetMode(i, 0, k);
1367 }
1368 for (j = 1; j < Q - i; ++j)
1369 {
1370 outarray[idx++] = GetMode(i, j, 0);
1371 }
1372 }
1373 }
1374}
1375
1376void StdTetExp::v_GetTraceCoeffMap(const unsigned int fid,
1377 Array<OneD, unsigned int> &maparray)
1378{
1379 int i, j, k;
1380 int P = 0, Q = 0, idx = 0;
1381 int nFaceCoeffs = 0;
1382
1383 switch (fid)
1384 {
1385 case 0:
1386 P = m_base[0]->GetNumModes();
1387 Q = m_base[1]->GetNumModes();
1388 break;
1389 case 1:
1390 P = m_base[0]->GetNumModes();
1391 Q = m_base[2]->GetNumModes();
1392 break;
1393 case 2:
1394 case 3:
1395 P = m_base[1]->GetNumModes();
1396 Q = m_base[2]->GetNumModes();
1397 break;
1398 default:
1399 ASSERTL0(false, "fid must be between 0 and 3");
1400 }
1401
1402 nFaceCoeffs = P * (2 * Q - P + 1) / 2;
1403
1404 if (maparray.size() != nFaceCoeffs)
1405 {
1406 maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1407 }
1408
1409 switch (fid)
1410 {
1411 case 0:
1412 idx = 0;
1413 for (i = 0; i < P; ++i)
1414 {
1415 for (j = 0; j < Q - i; ++j)
1416 {
1417 maparray[idx++] = GetMode(i, j, 0);
1418 }
1419 }
1420 break;
1421 case 1:
1422 idx = 0;
1423 for (i = 0; i < P; ++i)
1424 {
1425 for (k = 0; k < Q - i; ++k)
1426 {
1427 maparray[idx++] = GetMode(i, 0, k);
1428 }
1429 }
1430 break;
1431 case 2:
1432 idx = 0;
1433 for (j = 0; j < P - 1; ++j)
1434 {
1435 for (k = 0; k < Q - 1 - j; ++k)
1436 {
1437 maparray[idx++] = GetMode(1, j, k);
1438 // Incorporate modes from zeroth plane where needed.
1439 if (j == 0 && k == 0)
1440 {
1441 maparray[idx++] = GetMode(0, 0, 1);
1442 }
1443 if (j == 0 && k == Q - 2)
1444 {
1445 for (int r = 0; r < Q - 1; ++r)
1446 {
1447 maparray[idx++] = GetMode(0, 1, r);
1448 }
1449 }
1450 }
1451 }
1452 break;
1453 case 3:
1454 idx = 0;
1455 for (j = 0; j < P; ++j)
1456 {
1457 for (k = 0; k < Q - j; ++k)
1458 {
1459 maparray[idx++] = GetMode(0, j, k);
1460 }
1461 }
1462 break;
1463 default:
1464 ASSERTL0(false, "Element map not available.");
1465 }
1466}
1467
1468void StdTetExp::v_GetElmtTraceToTraceMap(const unsigned int fid,
1469 Array<OneD, unsigned int> &maparray,
1470 Array<OneD, int> &signarray,
1471 Orientation faceOrient, int P, int Q)
1472{
1473 int nummodesA = 0, nummodesB = 0, i, j, k, idx;
1474
1476 "Method only implemented for Modified_A BasisType (x "
1477 "direction), Modified_B BasisType (y direction), and "
1478 "Modified_C BasisType(z direction)");
1479
1480 int nFaceCoeffs = 0;
1481
1482 switch (fid)
1483 {
1484 case 0:
1485 nummodesA = m_base[0]->GetNumModes();
1486 nummodesB = m_base[1]->GetNumModes();
1487 break;
1488 case 1:
1489 nummodesA = m_base[0]->GetNumModes();
1490 nummodesB = m_base[2]->GetNumModes();
1491 break;
1492 case 2:
1493 case 3:
1494 nummodesA = m_base[1]->GetNumModes();
1495 nummodesB = m_base[2]->GetNumModes();
1496 break;
1497 default:
1498 ASSERTL0(false, "fid must be between 0 and 3");
1499 }
1500
1501 if (P == -1)
1502 {
1503 P = nummodesA;
1504 Q = nummodesB;
1505 }
1506
1507 nFaceCoeffs = P * (2 * Q - P + 1) / 2;
1508
1509 // Allocate the map array and sign array; set sign array to ones (+)
1510 if (maparray.size() != nFaceCoeffs)
1511 {
1512 maparray = Array<OneD, unsigned int>(nFaceCoeffs, 1);
1513 }
1514
1515 if (signarray.size() != nFaceCoeffs)
1516 {
1517 signarray = Array<OneD, int>(nFaceCoeffs, 1);
1518 }
1519 else
1520 {
1521 fill(signarray.data(), signarray.data() + nFaceCoeffs, 1);
1522 }
1523
1524 // zero signmap and set maparray to zero if elemental
1525 // modes are not as large as face modesl
1526 idx = 0;
1527 int cnt = 0;
1528 int minPA = min(nummodesA, P);
1529 int minQB = min(nummodesB, Q);
1530
1531 for (j = 0; j < minPA; ++j)
1532 {
1533 // set maparray
1534 for (k = 0; k < minQB - j; ++k, ++cnt)
1535 {
1536 maparray[idx++] = cnt;
1537 }
1538
1539 cnt += nummodesB - minQB;
1540
1541 for (k = nummodesB - j; k < Q - j; ++k)
1542 {
1543 signarray[idx] = 0.0;
1544 maparray[idx++] = maparray[0];
1545 }
1546 }
1547
1548 for (j = nummodesA; j < P; ++j)
1549 {
1550 for (k = 0; k < Q - j; ++k)
1551 {
1552 signarray[idx] = 0.0;
1553 maparray[idx++] = maparray[0];
1554 }
1555 }
1556
1557 if (faceOrient == eDir1BwdDir1_Dir2FwdDir2)
1558 {
1559 idx = 0;
1560 for (i = 0; i < P; ++i)
1561 {
1562 for (j = 0; j < Q - i; ++j, idx++)
1563 {
1564 if (i > 1)
1565 {
1566 signarray[idx] = (i % 2 ? -1 : 1);
1567 }
1568 }
1569 }
1570
1571 swap(maparray[0], maparray[Q]);
1572
1573 for (i = 1; i < Q - 1; ++i)
1574 {
1575 swap(maparray[i + 1], maparray[Q + i]);
1576 }
1577 }
1578}
1579
1580/**
1581 * Maps interior modes of an edge to the elemental modes.
1582 */
1584 const int eid, Array<OneD, unsigned int> &maparray,
1585 Array<OneD, int> &signarray, const Orientation edgeOrient)
1586{
1587 int i;
1588 const int P = m_base[0]->GetNumModes();
1589 const int Q = m_base[1]->GetNumModes();
1590 const int R = m_base[2]->GetNumModes();
1591
1592 const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid) - 2;
1593
1594 if (maparray.size() != nEdgeIntCoeffs)
1595 {
1596 maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1597 }
1598 else
1599 {
1600 fill(maparray.data(), maparray.data() + nEdgeIntCoeffs, 0);
1601 }
1602
1603 if (signarray.size() != nEdgeIntCoeffs)
1604 {
1605 signarray = Array<OneD, int>(nEdgeIntCoeffs, 1);
1606 }
1607 else
1608 {
1609 fill(signarray.data(), signarray.data() + nEdgeIntCoeffs, 1);
1610 }
1611
1612 switch (eid)
1613 {
1614 case 0:
1615 for (i = 0; i < P - 2; ++i)
1616 {
1617 maparray[i] = GetMode(i + 2, 0, 0);
1618 }
1619 if (edgeOrient == eBackwards)
1620 {
1621 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1622 {
1623 signarray[i] = -1;
1624 }
1625 }
1626 break;
1627 case 1:
1628 for (i = 0; i < Q - 2; ++i)
1629 {
1630 maparray[i] = GetMode(1, i + 1, 0);
1631 }
1632 if (edgeOrient == eBackwards)
1633 {
1634 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1635 {
1636 signarray[i] = -1;
1637 }
1638 }
1639 break;
1640 case 2:
1641 for (i = 0; i < Q - 2; ++i)
1642 {
1643 maparray[i] = GetMode(0, i + 2, 0);
1644 }
1645 if (edgeOrient == eBackwards)
1646 {
1647 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1648 {
1649 signarray[i] = -1;
1650 }
1651 }
1652 break;
1653 case 3:
1654 for (i = 0; i < R - 2; ++i)
1655 {
1656 maparray[i] = GetMode(0, 0, i + 2);
1657 }
1658 if (edgeOrient == eBackwards)
1659 {
1660 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1661 {
1662 signarray[i] = -1;
1663 }
1664 }
1665 break;
1666 case 4:
1667 for (i = 0; i < R - 2; ++i)
1668 {
1669 maparray[i] = GetMode(1, 0, i + 1);
1670 }
1671 if (edgeOrient == eBackwards)
1672 {
1673 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1674 {
1675 signarray[i] = -1;
1676 }
1677 }
1678 break;
1679 case 5:
1680 for (i = 0; i < R - 2; ++i)
1681 {
1682 maparray[i] = GetMode(0, 1, i + 1);
1683 }
1684 if (edgeOrient == eBackwards)
1685 {
1686 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1687 {
1688 signarray[i] = -1;
1689 }
1690 }
1691 break;
1692 default:
1693 ASSERTL0(false, "Edge not defined.");
1694 break;
1695 }
1696}
1697
1699 const int fid, Array<OneD, unsigned int> &maparray,
1700 Array<OneD, int> &signarray, const Orientation faceOrient)
1701{
1702 int i, j, idx, k;
1703 const int P = m_base[0]->GetNumModes();
1704 const int Q = m_base[1]->GetNumModes();
1705 const int R = m_base[2]->GetNumModes();
1706
1707 const int nFaceIntCoeffs = v_GetTraceIntNcoeffs(fid);
1708
1709 if (maparray.size() != nFaceIntCoeffs)
1710 {
1711 maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1712 }
1713
1714 if (signarray.size() != nFaceIntCoeffs)
1715 {
1716 signarray = Array<OneD, int>(nFaceIntCoeffs, 1);
1717 }
1718 else
1719 {
1720 fill(signarray.data(), signarray.data() + nFaceIntCoeffs, 1);
1721 }
1722
1723 switch (fid)
1724 {
1725 case 0:
1726 idx = 0;
1727 for (i = 2; i < P; ++i)
1728 {
1729 for (j = 1; j < Q - i; ++j)
1730 {
1731 if ((int)faceOrient == 7)
1732 {
1733 signarray[idx] = (i % 2 ? -1 : 1);
1734 }
1735 maparray[idx++] = GetMode(i, j, 0);
1736 }
1737 }
1738 break;
1739 case 1:
1740 idx = 0;
1741 for (i = 2; i < P; ++i)
1742 {
1743 for (k = 1; k < R - i; ++k)
1744 {
1745 if ((int)faceOrient == 7)
1746 {
1747 signarray[idx] = (i % 2 ? -1 : 1);
1748 }
1749 maparray[idx++] = GetMode(i, 0, k);
1750 }
1751 }
1752 break;
1753 case 2:
1754 idx = 0;
1755 for (j = 1; j < Q - 1; ++j)
1756 {
1757 for (k = 1; k < R - 1 - j; ++k)
1758 {
1759 if ((int)faceOrient == 7)
1760 {
1761 signarray[idx] = ((j + 1) % 2 ? -1 : 1);
1762 }
1763 maparray[idx++] = GetMode(1, j, k);
1764 }
1765 }
1766 break;
1767 case 3:
1768 idx = 0;
1769 for (j = 2; j < Q; ++j)
1770 {
1771 for (k = 1; k < R - j; ++k)
1772 {
1773 if ((int)faceOrient == 7)
1774 {
1775 signarray[idx] = (j % 2 ? -1 : 1);
1776 }
1777 maparray[idx++] = GetMode(0, j, k);
1778 }
1779 }
1780 break;
1781 default:
1782 ASSERTL0(false, "Face interior map not available.");
1783 break;
1784 }
1785}
1786//---------------------------------------
1787// Wrapper functions
1788//---------------------------------------
1790{
1791
1792 MatrixType mtype = mkey.GetMatrixType();
1793
1794 DNekMatSharedPtr Mat;
1795
1796 switch (mtype)
1797 {
1799 {
1800 int nq0 = m_base[0]->GetNumPoints();
1801 int nq1 = m_base[1]->GetNumPoints();
1802 int nq2 = m_base[2]->GetNumPoints();
1803 int nq;
1804
1805 // take definition from key
1807 {
1808 nq = (int)mkey.GetConstFactor(eFactorConst);
1809 }
1810 else
1811 {
1812 nq = max(nq0, max(nq1, nq2));
1813 }
1814
1815 int neq =
1818 Array<OneD, NekDouble> coll(3);
1820 Array<OneD, NekDouble> tmp(nq0);
1821
1822 Mat =
1823 MemoryManager<DNekMat>::AllocateSharedPtr(neq, nq0 * nq1 * nq2);
1824 int cnt = 0;
1825
1826 for (int i = 0; i < nq; ++i)
1827 {
1828 for (int j = 0; j < nq - i; ++j)
1829 {
1830 for (int k = 0; k < nq - i - j; ++k, ++cnt)
1831 {
1832 coords[cnt] = Array<OneD, NekDouble>(3);
1833 coords[cnt][0] = -1.0 + 2 * k / (NekDouble)(nq - 1);
1834 coords[cnt][1] = -1.0 + 2 * j / (NekDouble)(nq - 1);
1835 coords[cnt][2] = -1.0 + 2 * i / (NekDouble)(nq - 1);
1836 }
1837 }
1838 }
1839
1840 for (int i = 0; i < neq; ++i)
1841 {
1842 LocCoordToLocCollapsed(coords[i], coll);
1843
1844 I[0] = m_base[0]->GetI(coll);
1845 I[1] = m_base[1]->GetI(coll + 1);
1846 I[2] = m_base[2]->GetI(coll + 2);
1847
1848 // interpolate first coordinate direction
1849 NekDouble fac;
1850 for (int k = 0; k < nq2; ++k)
1851 {
1852 for (int j = 0; j < nq1; ++j)
1853 {
1854
1855 fac = (I[1]->GetPtr())[j] * (I[2]->GetPtr())[k];
1856 Vmath::Smul(nq0, fac, I[0]->GetPtr(), 1, tmp, 1);
1857
1858 Vmath::Vcopy(nq0, &tmp[0], 1,
1859 Mat->GetRawPtr() + k * nq0 * nq1 * neq +
1860 j * nq0 * neq + i,
1861 neq);
1862 }
1863 }
1864 }
1865 }
1866 break;
1867 default:
1868 {
1870 }
1871 break;
1872 }
1873
1874 return Mat;
1875}
1876
1878{
1879 return v_GenMatrix(mkey);
1880}
1881
1882//---------------------------------------
1883// Private helper functions
1884//---------------------------------------
1885
1886/**
1887 * @brief Compute the mode number in the expansion for a particular
1888 * tensorial combination.
1889 *
1890 * Modes are numbered with the r index travelling fastest, followed by
1891 * q and then p, and each q-r plane is of size
1892 * (Q+1)*(Q+2)/2+max(0,R-Q-p)*Q. For example, when P=2, Q=3 and R=4 (nm0=3, nm1
1893 * = 4, nm2 = 5) the indexing inside each q-r plane (with r increasing upwards
1894 * and q to the right) is:
1895 *
1896 * 4
1897 * 3 8 17
1898 * 2 7 11 16 20 25
1899 * 1 6 10 13 15 19 22 24 27
1900 * 0 5 9 12 14 18 21 23 26
1901 *
1902 * Geometrically they can be interpreted as
1903 * p = 0: p = 2: p = 1:
1904 * ----------------------------------
1905 * 1
1906 * 4 8 17
1907 * 3 11 7 25 16 20
1908 * 2 10 13 6 24 27 15 19 22
1909 * 0 9 12 5 23 26 14 18 21
1910 *
1911 * so we have the following breakdown
1912 *
1913 * Vertices V[0,1,2,3] = [0, 14, 5, 1]
1914 * Edges E[0,1,2,3,4,5,6] =[[23],[18, 21],
1915 * [9, 12], [2,3,4], [15, 16, 17], [6,7,8]]
1916 * Faces F[0.1,2,3] = [[26], [24,25], [19, 22, 20], [10, 13, 11]
1917 * Interior [27]
1918 * Note that in this element, we must have that \f$ P \leq Q \leq
1919 * R\f$.
1920 */
1921int StdTetExp::GetMode(const int I, const int J, const int K)
1922{
1923 const int Q = m_base[1]->GetNumModes();
1924 const int R = m_base[2]->GetNumModes();
1925
1926 int i, j, q_hat, k_hat;
1927 int cnt = 0;
1928
1929 // Traverse to q-r plane number I
1930 for (i = 0; i < I; ++i)
1931 {
1932 // Size of triangle part
1933 q_hat = Q - i;
1934 // Size of rectangle part
1935 k_hat = R - Q;
1936 cnt += q_hat * (q_hat + 1) / 2 + k_hat * (Q - i);
1937 }
1938
1939 // Traverse to q column J
1940 q_hat = R - I;
1941 for (j = 0; j < J; ++j)
1942 {
1943 cnt += q_hat;
1944 q_hat--;
1945 }
1946
1947 // Traverse up stacks to K
1948 cnt += K;
1949
1950 return cnt;
1951}
1952
1954 const Array<OneD, const NekDouble> &inarray,
1955 Array<OneD, NekDouble> &outarray)
1956{
1957 int i, j;
1958
1959 int nquad0 = m_base[0]->GetNumPoints();
1960 int nquad1 = m_base[1]->GetNumPoints();
1961 int nquad2 = m_base[2]->GetNumPoints();
1962
1963 const Array<OneD, const NekDouble> &w0 = m_base[0]->GetW();
1964 const Array<OneD, const NekDouble> &w1 = m_base[1]->GetW();
1965 const Array<OneD, const NekDouble> &w2 = m_base[2]->GetW();
1966
1967 const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
1968 const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
1969
1970 // multiply by integration constants
1971 for (i = 0; i < nquad1 * nquad2; ++i)
1972 {
1973 Vmath::Vmul(nquad0, (NekDouble *)&inarray[0] + i * nquad0, 1, w0.data(),
1974 1, &outarray[0] + i * nquad0, 1);
1975 }
1976
1977 switch (m_base[1]->GetPointsType())
1978 {
1979 // (1,0) Jacobi Inner product.
1980 case LibUtilities::eGaussRadauMAlpha1Beta0:
1981 for (j = 0; j < nquad2; ++j)
1982 {
1983 for (i = 0; i < nquad1; ++i)
1984 {
1985 Blas::Dscal(nquad0, 0.5 * w1[i],
1986 &outarray[0] + i * nquad0 + j * nquad0 * nquad1,
1987 1);
1988 }
1989 }
1990 break;
1991
1992 default:
1993 for (j = 0; j < nquad2; ++j)
1994 {
1995 for (i = 0; i < nquad1; ++i)
1996 {
1997 Blas::Dscal(nquad0, 0.5 * (1 - z1[i]) * w1[i],
1998 &outarray[0] + i * nquad0 + j * nquad0 * nquad1,
1999 1);
2000 }
2001 }
2002 break;
2003 }
2004
2005 switch (m_base[2]->GetPointsType())
2006 {
2007 // (2,0) Jacobi inner product.
2008 case LibUtilities::eGaussRadauMAlpha2Beta0:
2009 for (i = 0; i < nquad2; ++i)
2010 {
2011 Blas::Dscal(nquad0 * nquad1, 0.25 * w2[i],
2012 &outarray[0] + i * nquad0 * nquad1, 1);
2013 }
2014 break;
2015 // (1,0) Jacobi inner product.
2016 case LibUtilities::eGaussRadauMAlpha1Beta0:
2017 for (i = 0; i < nquad2; ++i)
2018 {
2019 Blas::Dscal(nquad0 * nquad1, 0.25 * (1 - z2[i]) * w2[i],
2020 &outarray[0] + i * nquad0 * nquad1, 1);
2021 }
2022 break;
2023 default:
2024 for (i = 0; i < nquad2; ++i)
2025 {
2026 Blas::Dscal(nquad0 * nquad1,
2027 0.25 * (1 - z2[i]) * (1 - z2[i]) * w2[i],
2028 &outarray[0] + i * nquad0 * nquad1, 1);
2029 }
2030 break;
2031 }
2032}
2033
2035 const StdMatrixKey &mkey)
2036{
2037 // To do : 1) add a test to ensure 0 \leq SvvCutoff \leq 1.
2038 // 2) check if the transfer function needs an analytical
2039 // Fourier transform.
2040 // 3) if it doesn't : find a transfer function that renders
2041 // the if( cutoff_a ...) useless to reduce computational
2042 // cost.
2043 // 4) add SVVDiffCoef to both models!!
2044
2045 int qa = m_base[0]->GetNumPoints();
2046 int qb = m_base[1]->GetNumPoints();
2047 int qc = m_base[2]->GetNumPoints();
2048 int nmodes_a = m_base[0]->GetNumModes();
2049 int nmodes_b = m_base[1]->GetNumModes();
2050 int nmodes_c = m_base[2]->GetNumModes();
2051
2052 // Declare orthogonal basis.
2056
2060
2061 StdTetExp OrthoExp(Ba, Bb, Bc);
2062
2063 Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2064 int i, j, k, cnt = 0;
2065
2066 // project onto physical space.
2067 OrthoExp.FwdTrans(array, orthocoeffs);
2068
2070 {
2071 // Rodrigo's power kernel
2073 NekDouble SvvDiffCoeff =
2076
2077 for (i = 0; i < nmodes_a; ++i)
2078 {
2079 for (j = 0; j < nmodes_b - j; ++j)
2080 {
2081 NekDouble fac1 = std::max(
2082 pow((1.0 * i) / (nmodes_a - 1), cutoff * nmodes_a),
2083 pow((1.0 * j) / (nmodes_b - 1), cutoff * nmodes_b));
2084
2085 for (k = 0; k < nmodes_c - i - j; ++k)
2086 {
2087 NekDouble fac =
2088 std::max(fac1, pow((1.0 * k) / (nmodes_c - 1),
2089 cutoff * nmodes_c));
2090
2091 orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2092 cnt++;
2093 }
2094 }
2095 }
2096 }
2097 else if (mkey.ConstFactorExists(
2098 eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
2099 {
2102
2103 int max_abc = max(nmodes_a - kSVVDGFiltermodesmin,
2104 nmodes_b - kSVVDGFiltermodesmin);
2105 max_abc = max(max_abc, nmodes_c - kSVVDGFiltermodesmin);
2106 // clamp max_abc
2107 max_abc = max(max_abc, 0);
2108 max_abc = min(max_abc, kSVVDGFiltermodesmax - kSVVDGFiltermodesmin);
2109
2110 for (i = 0; i < nmodes_a; ++i)
2111 {
2112 for (j = 0; j < nmodes_b - j; ++j)
2113 {
2114 int maxij = max(i, j);
2115
2116 for (k = 0; k < nmodes_c - i - j; ++k)
2117 {
2118 int maxijk = max(maxij, k);
2119 maxijk = min(maxijk, kSVVDGFiltermodesmax - 1);
2120
2121 orthocoeffs[cnt] *=
2122 SvvDiffCoeff * kSVVDGFilter[max_abc][maxijk];
2123 cnt++;
2124 }
2125 }
2126 }
2127 }
2128 else
2129 {
2130
2131 // SVV filter paramaters (how much added diffusion
2132 // relative to physical one and fraction of modes from
2133 // which you start applying this added diffusion)
2134
2135 NekDouble SvvDiffCoeff =
2137 NekDouble SVVCutOff =
2139
2140 // Defining the cut of mode
2141 int cutoff_a = (int)(SVVCutOff * nmodes_a);
2142 int cutoff_b = (int)(SVVCutOff * nmodes_b);
2143 int cutoff_c = (int)(SVVCutOff * nmodes_c);
2144 int nmodes = min(min(nmodes_a, nmodes_b), nmodes_c);
2145 NekDouble cutoff = min(min(cutoff_a, cutoff_b), cutoff_c);
2146 NekDouble epsilon = 1;
2147
2148 //------"New" Version August 22nd '13--------------------
2149 for (i = 0; i < nmodes_a; ++i)
2150 {
2151 for (j = 0; j < nmodes_b - i; ++j)
2152 {
2153 for (k = 0; k < nmodes_c - i - j; ++k)
2154 {
2155 if (i + j + k >= cutoff)
2156 {
2157 orthocoeffs[cnt] *= ((SvvDiffCoeff)*exp(
2158 -(i + j + k - nmodes) * (i + j + k - nmodes) /
2159 ((NekDouble)((i + j + k - cutoff + epsilon) *
2160 (i + j + k - cutoff + epsilon)))));
2161 }
2162 else
2163 {
2164 orthocoeffs[cnt] *= 0.0;
2165 }
2166 cnt++;
2167 }
2168 }
2169 }
2170 }
2171
2172 // backward transform to physical space
2173 OrthoExp.BwdTrans(orthocoeffs, array);
2174}
2175
2177 const Array<OneD, const NekDouble> &inarray,
2178 Array<OneD, NekDouble> &outarray)
2179{
2180 int nquad0 = m_base[0]->GetNumPoints();
2181 int nquad1 = m_base[1]->GetNumPoints();
2182 int nquad2 = m_base[2]->GetNumPoints();
2183 int nqtot = nquad0 * nquad1 * nquad2;
2184 int nmodes0 = m_base[0]->GetNumModes();
2185 int nmodes1 = m_base[1]->GetNumModes();
2186 int nmodes2 = m_base[2]->GetNumModes();
2187 int numMax = nmodes0;
2188
2190 Array<OneD, NekDouble> coeff_tmp1(m_ncoeffs, 0.0);
2191 Array<OneD, NekDouble> coeff_tmp2(m_ncoeffs, 0.0);
2192 Array<OneD, NekDouble> phys_tmp(nqtot, 0.0);
2193 Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
2194
2195 Vmath::Vcopy(m_ncoeffs, inarray, 1, coeff_tmp2, 1);
2196
2197 const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
2198 const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
2199 const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
2200
2201 LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A, nmodes0, Pkey0);
2202 LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_B, nmodes1, Pkey1);
2203 LibUtilities::BasisKey bortho2(LibUtilities::eOrtho_C, nmodes2, Pkey2);
2204
2205 Vmath::Zero(m_ncoeffs, coeff_tmp2, 1);
2206
2209 bortho0, bortho1, bortho2);
2210
2211 BwdTrans(inarray, phys_tmp);
2212 OrthoTetExp->FwdTrans(phys_tmp, coeff);
2213
2214 Vmath::Zero(m_ncoeffs, outarray, 1);
2215
2216 // filtering
2217 int cnt = 0;
2218 for (int u = 0; u < numMin; ++u)
2219 {
2220 for (int i = 0; i < numMin - u; ++i)
2221 {
2222 Vmath::Vcopy(numMin - u - i, tmp = coeff + cnt, 1,
2223 tmp2 = coeff_tmp1 + cnt, 1);
2224 cnt += numMax - u - i;
2225 }
2226 for (int i = numMin; i < numMax - u; ++i)
2227 {
2228 cnt += numMax - u - i;
2229 }
2230 }
2231
2232 OrthoTetExp->BwdTrans(coeff_tmp1, phys_tmp);
2233 FwdTrans(phys_tmp, outarray);
2234}
2235
2237 Array<OneD, int> &conn, [[maybe_unused]] bool standard)
2238{
2239 int np0 = m_base[0]->GetNumPoints();
2240 int np1 = m_base[1]->GetNumPoints();
2241 int np2 = m_base[2]->GetNumPoints();
2242 int np = max(np0, max(np1, np2));
2243
2244 conn = Array<OneD, int>(4 * (np - 1) * (np - 1) * (np - 1));
2245
2246 int row = 0;
2247 int rowp1 = 0;
2248 int plane = 0;
2249 int row1 = 0;
2250 int row1p1 = 0;
2251 int planep1 = 0;
2252 int cnt = 0;
2253 for (int i = 0; i < np - 1; ++i)
2254 {
2255 planep1 += (np - i) * (np - i + 1) / 2;
2256 row = 0; // current plane row offset
2257 rowp1 = 0; // current plane row plus one offset
2258 row1 = 0; // next plane row offset
2259 row1p1 = 0; // nex plane row plus one offset
2260 for (int j = 0; j < np - i - 1; ++j)
2261 {
2262 rowp1 += np - i - j;
2263 row1p1 += np - i - j - 1;
2264 for (int k = 0; k < np - i - j - 2; ++k)
2265 {
2266 conn[cnt++] = plane + row + k + 1;
2267 conn[cnt++] = plane + row + k;
2268 conn[cnt++] = plane + rowp1 + k;
2269 conn[cnt++] = planep1 + row1 + k;
2270
2271 conn[cnt++] = plane + row + k + 1;
2272 conn[cnt++] = plane + rowp1 + k + 1;
2273 conn[cnt++] = planep1 + row1 + k + 1;
2274 conn[cnt++] = planep1 + row1 + k;
2275
2276 conn[cnt++] = plane + rowp1 + k + 1;
2277 conn[cnt++] = plane + row + k + 1;
2278 conn[cnt++] = plane + rowp1 + k;
2279 conn[cnt++] = planep1 + row1 + k;
2280
2281 conn[cnt++] = planep1 + row1 + k;
2282 conn[cnt++] = planep1 + row1p1 + k;
2283 conn[cnt++] = plane + rowp1 + k;
2284 conn[cnt++] = plane + rowp1 + k + 1;
2285
2286 conn[cnt++] = planep1 + row1 + k;
2287 conn[cnt++] = planep1 + row1p1 + k;
2288 conn[cnt++] = planep1 + row1 + k + 1;
2289 conn[cnt++] = plane + rowp1 + k + 1;
2290
2291 if (k < np - i - j - 3)
2292 {
2293 conn[cnt++] = plane + rowp1 + k + 1;
2294 conn[cnt++] = planep1 + row1p1 + k + 1;
2295 conn[cnt++] = planep1 + row1 + k + 1;
2296 conn[cnt++] = planep1 + row1p1 + k;
2297 }
2298 }
2299
2300 conn[cnt++] = plane + row + np - i - j - 1;
2301 conn[cnt++] = plane + row + np - i - j - 2;
2302 conn[cnt++] = plane + rowp1 + np - i - j - 2;
2303 conn[cnt++] = planep1 + row1 + np - i - j - 2;
2304
2305 row += np - i - j;
2306 row1 += np - i - j - 1;
2307 }
2308 plane += (np - i) * (np - i + 1) / 2;
2309 }
2310}
2311
2312} // namespace Nektar::StdRegions
#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
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:265
Describes the specification for a Basis.
Definition: Basis.h:45
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:74
Defines a specification for a set of points.
Definition: Points.h:50
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
void BwdTrans_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 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)
NekDouble BaryTensorDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs)
void PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d1, Array< OneD, NekDouble > &outarray_d2, Array< OneD, NekDouble > &outarray_d3)
Calculate the 3D derivative in the local tensor/collapsed coordinate at the physical points.
The base class for all shapes.
Definition: StdExpansion.h:65
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
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:612
void LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Convert local cartesian coordinate xi into local collapsed coordinates eta.
const Array< OneD, const LibUtilities::BasisSharedPtr > & GetBase() const
This function gets the shared point to basis.
Definition: StdExpansion.h:100
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
NekDouble PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals)
This function evaluates the expansion at a single (arbitrary) point of the domain.
Definition: StdExpansion.h:928
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
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
Definition: StdExpansion.h:169
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
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:83
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:124
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:133
int v_GetNtraces() const override
Definition: StdTetExp.cpp:967
int v_GetTraceNcoeffs(const int i) const override
Definition: StdTetExp.cpp:1017
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
Definition: StdTetExp.cpp:858
void v_GetTraceNumModes(const int fid, int &numModes0, int &numModes1, Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
Definition: StdTetExp.cpp:923
int v_GetTraceIntNcoeffs(const int i) const override
Definition: StdTetExp.cpp:1043
void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray) override
Definition: StdTetExp.cpp:809
void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdTetExp.cpp:416
void v_GetInteriorMap(Array< OneD, unsigned int > &outarray) override
Definition: StdTetExp.cpp:1282
int v_GetNedges() const override
Definition: StdTetExp.cpp:962
int GetMode(const int i, const int j, const int k)
Compute the mode number in the expansion for a particular tensorial combination.
Definition: StdTetExp.cpp:1921
int v_GetEdgeNcoeffs(const int i) const override
Definition: StdTetExp.cpp:1082
bool v_IsBoundaryInteriorExpansion() const override
Definition: StdTetExp.cpp:1191
LibUtilities::PointsKey v_GetTracePointsKey(const int i, const int j) const override
Definition: StdTetExp.cpp:1103
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) final
Definition: StdTetExp.cpp:816
void v_GetEdgeInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
Definition: StdTetExp.cpp:1583
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdTetExp.cpp:1953
StdTetExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc)
Definition: StdTetExp.cpp:42
int v_GetTraceNumPoints(const int i) const override
Definition: StdTetExp.cpp:1064
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdTetExp.cpp:618
void v_BwdTrans_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) override
Definition: StdTetExp.cpp:322
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey) override
Definition: StdTetExp.cpp:2034
DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey) override
Definition: StdTetExp.cpp:1877
void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta) override
Definition: StdTetExp.cpp:769
void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray) override
Definition: StdTetExp.cpp:1321
void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true) override
Definition: StdTetExp.cpp:2236
void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdTetExp.cpp:466
int v_NumBndryCoeffs() const override
Definition: StdTetExp.cpp:977
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
Definition: StdTetExp.cpp:493
void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdTetExp.cpp:266
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdTetExp.cpp:605
void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdTetExp.cpp:2176
void v_StdPhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2) override
Definition: StdTetExp.cpp:226
void v_GetCoords(Array< OneD, NekDouble > &coords_x, Array< OneD, NekDouble > &coords_y, Array< OneD, NekDouble > &coords_z) override
Definition: StdTetExp.cpp:1161
void v_GetTraceCoeffMap(const unsigned int fid, Array< OneD, unsigned int > &maparray) override
Definition: StdTetExp.cpp:1376
void v_GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
Definition: StdTetExp.cpp:1698
void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdTetExp.cpp:293
void v_GetElmtTraceToTraceMap(const unsigned int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards, int P=-1, int Q=-1) override
Definition: StdTetExp.cpp:1468
int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false) override
Definition: StdTetExp.cpp:1201
int v_GetNverts() const override
Definition: StdTetExp.cpp:957
void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_dx, Array< OneD, NekDouble > &out_dy, Array< OneD, NekDouble > &out_dz) override
Calculate the derivative of the physical points.
Definition: StdTetExp.cpp:83
const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int k, bool UseGLL=false) const override
Definition: StdTetExp.cpp:1134
LibUtilities::ShapeType v_DetShapeType() const override
Definition: StdTetExp.cpp:972
void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi) override
Definition: StdTetExp.cpp:801
LibUtilities::ShapeType DetShapeType() const
Definition: StdTetExp.h:56
int v_NumDGBndryCoeffs() const override
Definition: StdTetExp.cpp:996
void v_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) override
Definition: StdTetExp.cpp:523
int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset) override
Definition: StdTetExp.cpp:1123
DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey) override
Definition: StdTetExp.cpp:1789
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 = alpha A x plus beta y where A[m x n].
Definition: Blas.hpp:211
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:149
static double Ddot(const int &n, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: output = .
Definition: Blas.hpp:163
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition: Blas.hpp:383
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
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:210
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:187
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
@ P
Monomial polynomials .
Definition: BasisType.h:62
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:42
@ eModified_C
Principle Modified Functions .
Definition: BasisType.h:50
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:56
@ eOrtho_C
Principle Orthogonal Functions .
Definition: BasisType.h:46
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:44
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
static const NekDouble kNekZeroTol
std::shared_ptr< StdTetExp > StdTetExpSharedPtr
Definition: StdTetExp.h:227
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes, bool UseGLL)
const int kSVVDGFiltermodesmin
Definition: StdRegions.hpp:500
const int kSVVDGFiltermodesmax
Definition: StdRegions.hpp:501
const NekDouble kSVVDGFilter[9][11]
Definition: StdRegions.hpp:503
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
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 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 Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.hpp:194
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