Nektar++
StdPyrExp.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: StdPyrExp.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: pyramidic routines built upon StdExpansion3D
32//
33///////////////////////////////////////////////////////////////////////////////
34
37#include <iomanip>
38
39using namespace std;
40
41namespace Nektar::StdRegions
42{
44 const LibUtilities::BasisKey &Bb,
45 const LibUtilities::BasisKey &Bc)
46 : StdExpansion(LibUtilities::StdPyrData::getNumberOfCoefficients(
47 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
48 3, Ba, Bb, Bc),
49 StdExpansion3D(LibUtilities::StdPyrData::getNumberOfCoefficients(
50 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
51 Ba, Bb, Bc)
52{
53
54 ASSERTL0(Ba.GetNumModes() <= Bc.GetNumModes(),
55 "order in 'a' direction is higher "
56 "than order in 'c' direction");
57 ASSERTL0(Bb.GetNumModes() <= Bc.GetNumModes(),
58 "order in 'b' direction is higher "
59 "than order in 'c' direction");
62 "Expected basis type in 'c' direction to be ModifiedPyr_C or "
63 "OrthoPyr_C");
64}
65
66//---------------------------------------
67// Differentiation/integration Methods
68//---------------------------------------
69
70/**
71 * \brief Calculate the derivative of the physical points
72 *
73 * The derivative is evaluated at the nodal physical points.
74 * Derivatives with respect to the local Cartesian coordinates.
75 *
76 * \f$\begin{Bmatrix} \frac {\partial} {\partial \xi_1} \\ \frac
77 * {\partial} {\partial \xi_2} \\ \frac {\partial} {\partial \xi_3}
78 * \end{Bmatrix} = \begin{Bmatrix} \frac 2 {(1-\eta_3)} \frac \partial
79 * {\partial \bar \eta_1} \\ \frac {\partial} {\partial \xi_2} \ \
80 * \frac {(1 + \bar \eta_1)} {(1 - \eta_3)} \frac \partial {\partial
81 * \bar \eta_1} + \frac {\partial} {\partial \eta_3} \end{Bmatrix}\f$
82 */
84 Array<OneD, NekDouble> &out_dxi1,
85 Array<OneD, NekDouble> &out_dxi2,
86 Array<OneD, NekDouble> &out_dxi3)
87{
88 // PhysDerivative implementation based on Spen's book page 152.
89 int Qx = m_base[0]->GetNumPoints();
90 int Qy = m_base[1]->GetNumPoints();
91 int Qz = m_base[2]->GetNumPoints();
92
93 Array<OneD, NekDouble> dEta_bar1(Qx * Qy * Qz, 0.0);
94 Array<OneD, NekDouble> dXi2(Qx * Qy * Qz, 0.0);
95 Array<OneD, NekDouble> dEta3(Qx * Qy * Qz, 0.0);
96 PhysTensorDeriv(u_physical, dEta_bar1, dXi2, dEta3);
97
98 Array<OneD, const NekDouble> eta_x, eta_y, eta_z;
99 eta_x = m_base[0]->GetZ();
100 eta_y = m_base[1]->GetZ();
101 eta_z = m_base[2]->GetZ();
102
103 int i, j, k, n;
104
105 if (out_dxi1.size() > 0)
106 {
107 for (k = 0, n = 0; k < Qz; ++k)
108 {
109 NekDouble fac = 2.0 / (1.0 - eta_z[k]);
110 for (j = 0; j < Qy; ++j)
111 {
112 for (i = 0; i < Qx; ++i, ++n)
113 {
114 out_dxi1[n] = fac * dEta_bar1[n];
115 }
116 }
117 }
118 }
119
120 if (out_dxi2.size() > 0)
121 {
122 for (k = 0, n = 0; k < Qz; ++k)
123 {
124 NekDouble fac = 2.0 / (1.0 - eta_z[k]);
125 for (j = 0; j < Qy; ++j)
126 {
127 for (i = 0; i < Qx; ++i, ++n)
128 {
129 out_dxi2[n] = fac * dXi2[n];
130 }
131 }
132 }
133 }
134
135 if (out_dxi3.size() > 0)
136 {
137 for (k = 0, n = 0; k < Qz; ++k)
138 {
139 NekDouble fac = 1.0 / (1.0 - eta_z[k]);
140 for (j = 0; j < Qy; ++j)
141 {
142 NekDouble fac1 = (1.0 + eta_y[j]);
143 for (i = 0; i < Qx; ++i, ++n)
144 {
145 out_dxi3[n] = (1.0 + eta_x[i]) * fac * dEta_bar1[n] +
146 fac1 * fac * dXi2[n] + dEta3[n];
147 }
148 }
149 }
150 }
151}
152
153void StdPyrExp::v_PhysDeriv(const int dir,
154 const Array<OneD, const NekDouble> &inarray,
155 Array<OneD, NekDouble> &outarray)
156{
157 switch (dir)
158 {
159 case 0:
160 {
161 v_PhysDeriv(inarray, outarray, NullNekDouble1DArray,
163 break;
164 }
165
166 case 1:
167 {
168 v_PhysDeriv(inarray, NullNekDouble1DArray, outarray,
170 break;
171 }
172
173 case 2:
174 {
176 outarray);
177 break;
178 }
179
180 default:
181 {
182 ASSERTL1(false, "input dir is out of range");
183 }
184 break;
185 }
186}
187
192{
193 StdPyrExp::v_PhysDeriv(inarray, out_d0, out_d1, out_d2);
194}
195
196void StdPyrExp::v_StdPhysDeriv(const int dir,
197 const Array<OneD, const NekDouble> &inarray,
198 Array<OneD, NekDouble> &outarray)
199{
200 StdPyrExp::v_PhysDeriv(dir, inarray, outarray);
201}
202
203//---------------------------------------
204// Transforms
205//---------------------------------------
206
207/**
208 * \brief Backward transformation is evaluated at the quadrature
209 * points.
210 *
211 * \f$ u^{\delta} (\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{m(pqr)} \hat
212 * u_{pqr} \phi_{pqr} (\xi_{1i}, \xi_{2j}, \xi_{3k})\f$
213 *
214 * Backward transformation is three dimensional tensorial expansion
215 *
216 * \f$ u (\xi_{1i}, \xi_{2j}, \xi_{3k}) = \sum_{p=0}^{Q_x} \psi_p^a
217 * (\xi_{1i}) \lbrace { \sum_{q=0}^{Q_y} \psi_{q}^a (\xi_{2j})
218 * \lbrace { \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{pqr}^c (\xi_{3k})
219 * \rbrace} \rbrace}. \f$
220 *
221 * And sumfactorizing step of the form is as:\ \ \f$ f_{pqr}
222 * (\xi_{3k}) = \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{pqr}^c
223 * (\xi_{3k}),\\ g_{p} (\xi_{2j}, \xi_{3k}) = \sum_{r=0}^{Q_y}
224 * \psi_{p}^a (\xi_{2j}) f_{pqr} (\xi_{3k}),\\ u(\xi_{1i}, \xi_{2j},
225 * \xi_{3k}) = \sum_{p=0}^{Q_x} \psi_{p}^a (\xi_{1i}) g_{p}
226 * (\xi_{2j}, \xi_{3k}). \f$
227 **/
229 Array<OneD, NekDouble> &outarray)
230{
231 if (m_base[0]->Collocation() && m_base[1]->Collocation() &&
232 m_base[2]->Collocation())
233 {
235 m_base[2]->GetNumPoints(),
236 inarray, 1, outarray, 1);
237 }
238 else
239 {
240 StdPyrExp::v_BwdTrans_SumFac(inarray, outarray);
241 }
242}
243
244/**
245 * Sum-factorisation implementation of the BwdTrans operation.
246 */
248 Array<OneD, NekDouble> &outarray)
249{
250 int nquad0 = m_base[0]->GetNumPoints();
251 int nquad1 = m_base[1]->GetNumPoints();
252 int nquad2 = m_base[2]->GetNumPoints();
253 int order0 = m_base[0]->GetNumModes();
254 int order1 = m_base[1]->GetNumModes();
255
256 Array<OneD, NekDouble> wsp(nquad2 * order0 * order1 +
257 nquad2 * nquad1 * nquad0);
258
259 v_BwdTrans_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(),
260 m_base[2]->GetBdata(), inarray, outarray, wsp, true,
261 true, true);
262}
263
265 const Array<OneD, const NekDouble> &base0,
266 const Array<OneD, const NekDouble> &base1,
267 const Array<OneD, const NekDouble> &base2,
268 const Array<OneD, const NekDouble> &inarray,
270 [[maybe_unused]] bool doCheckCollDir0,
271 [[maybe_unused]] bool doCheckCollDir1,
272 [[maybe_unused]] bool doCheckCollDir2)
273{
274 int nquad0 = m_base[0]->GetNumPoints();
275 int nquad1 = m_base[1]->GetNumPoints();
276 int nquad2 = m_base[2]->GetNumPoints();
277
278 int order0 = m_base[0]->GetNumModes();
279 int order1 = m_base[1]->GetNumModes();
280 int order2 = m_base[2]->GetNumModes();
281
282 Array<OneD, NekDouble> tmp = wsp;
283 Array<OneD, NekDouble> tmp1 = tmp + nquad2 * order0 * order1;
284
285 int i, j, mode, mode1, cnt;
286
287 // Perform summation over '2' direction
288 mode = mode1 = cnt = 0;
289 for (i = 0; i < order0; ++i)
290 {
291 for (j = 0; j < order1; ++j, ++cnt)
292 {
293 int ijmax = max(i, j);
294 Blas::Dgemv('N', nquad2, order2 - ijmax, 1.0,
295 base2.data() + mode * nquad2, nquad2,
296 inarray.data() + mode1, 1, 0.0,
297 tmp.data() + cnt * nquad2, 1);
298 mode += order2 - ijmax;
299 mode1 += order2 - ijmax;
300 }
301 // increment mode in case order1!=order2
302 for (j = order1; j < order2; ++j)
303 {
304 int ijmax = max(i, j);
305 mode += order2 - ijmax;
306 }
307 }
308
309 // fix for modified basis by adding split of top singular
310 // vertex mode - currently (1+c)/2 x (1-b)/2 x (1-a)/2
311 // component is evaluated
313 {
314
315 // Not sure why we could not use basis as 1.0
316 // top singular vertex - (1+c)/2 x (1+b)/2 x (1-a)/2 component
317 Blas::Daxpy(nquad2, inarray[1], base2.data() + nquad2, 1,
318 &tmp[0] + nquad2, 1);
319
320 // top singular vertex - (1+c)/2 x (1-b)/2 x (1+a)/2 component
321 Blas::Daxpy(nquad2, inarray[1], base2.data() + nquad2, 1,
322 &tmp[0] + order1 * nquad2, 1);
323
324 // top singular vertex - (1+c)/2 x (1+b)/2 x (1+a)/2 component
325 Blas::Daxpy(nquad2, inarray[1], base2.data() + nquad2, 1,
326 &tmp[0] + order1 * nquad2 + nquad2, 1);
327 }
328
329 // Perform summation over '1' direction
330 mode = 0;
331 for (i = 0; i < order0; ++i)
332 {
333 Blas::Dgemm('N', 'T', nquad1, nquad2, order1, 1.0, base1.data(), nquad1,
334 tmp.data() + mode * nquad2, nquad2, 0.0,
335 tmp1.data() + i * nquad1 * nquad2, nquad1);
336 mode += order1;
337 }
338
339 // Perform summation over '0' direction
340 Blas::Dgemm('N', 'T', nquad0, nquad1 * nquad2, order0, 1.0, base0.data(),
341 nquad0, tmp1.data(), nquad1 * nquad2, 0.0, outarray.data(),
342 nquad0);
343}
344
345/** \brief Forward transform from physical quadrature space
346 stored in \a inarray and evaluate the expansion coefficients and
347 store in \a outarray
348
349 Inputs:\n
350
351 - \a inarray: array of physical quadrature points to be transformed
352
353 Outputs:\n
354
355 - \a outarray: updated array of expansion coefficients.
356
357*/
359 Array<OneD, NekDouble> &outarray)
360{
361 v_IProductWRTBase(inarray, outarray);
362
363 // get Mass matrix inverse
364 StdMatrixKey imasskey(eInvMass, DetShapeType(), *this);
365 DNekMatSharedPtr imatsys = GetStdMatrix(imasskey);
366
367 // copy inarray in case inarray == outarray
368 DNekVec in(m_ncoeffs, outarray);
369 DNekVec out(m_ncoeffs, outarray, eWrapper);
370
371 out = (*imatsys) * in;
372}
373
374//---------------------------------------
375// Inner product functions
376//---------------------------------------
377
378/** \brief Inner product of \a inarray over region with respect to the
379 expansion basis m_base[0]->GetBdata(),m_base[1]->GetBdata(),
380 m_base[2]->GetBdata() and return in \a outarray
381
382 Wrapper call to StdPyrExp::IProductWRTBase
383
384 Input:\n
385
386 - \a inarray: array of function evaluated at the physical collocation points
387
388 Output:\n
389
390 - \a outarray: array of inner product with respect to each basis over region
391
392*/
394 Array<OneD, NekDouble> &outarray)
395{
396 if (m_base[0]->Collocation() && m_base[1]->Collocation() &&
397 m_base[2]->Collocation())
398 {
399 v_MultiplyByStdQuadratureMetric(inarray, outarray);
400 }
401 else
402 {
403 StdPyrExp::v_IProductWRTBase_SumFac(inarray, outarray);
404 }
405}
406
408 const Array<OneD, const NekDouble> &inarray,
409 Array<OneD, NekDouble> &outarray, bool multiplybyweights)
410{
411
412 int nquad1 = m_base[1]->GetNumPoints();
413 int nquad2 = m_base[2]->GetNumPoints();
414 int order0 = m_base[0]->GetNumModes();
415 int order1 = m_base[1]->GetNumModes();
416
417 Array<OneD, NekDouble> wsp(order0 * nquad2 * (nquad1 + order1));
418
419 if (multiplybyweights)
420 {
421 Array<OneD, NekDouble> tmp(inarray.size());
422
424
426 m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
427 tmp, outarray, wsp, true, true, true);
428 }
429 else
430 {
432 m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
433 inarray, outarray, wsp, true, true, true);
434 }
435}
436
438 const Array<OneD, const NekDouble> &base0,
439 const Array<OneD, const NekDouble> &base1,
440 const Array<OneD, const NekDouble> &base2,
441 const Array<OneD, const NekDouble> &inarray,
443 [[maybe_unused]] bool doCheckCollDir0,
444 [[maybe_unused]] bool doCheckCollDir1,
445 [[maybe_unused]] bool doCheckCollDir2)
446{
447 int nquad0 = m_base[0]->GetNumPoints();
448 int nquad1 = m_base[1]->GetNumPoints();
449 int nquad2 = m_base[2]->GetNumPoints();
450
451 int order0 = m_base[0]->GetNumModes();
452 int order1 = m_base[1]->GetNumModes();
453 int order2 = m_base[2]->GetNumModes();
454
455 ASSERTL1(wsp.size() >= nquad1 * nquad2 * order0 + nquad2 * order0 * order1,
456 "Insufficient workspace size");
457
458 Array<OneD, NekDouble> tmp1 = wsp;
459 Array<OneD, NekDouble> tmp2 = wsp + nquad1 * nquad2 * order0;
460
461 int i, j, mode, mode1, cnt;
462
463 // Inner product with respect to the '0' direction
464 Blas::Dgemm('T', 'N', nquad1 * nquad2, order0, nquad0, 1.0, inarray.data(),
465 nquad0, base0.data(), nquad0, 0.0, tmp1.data(),
466 nquad1 * nquad2);
467
468 // Inner product with respect to the '1' direction
469 for (mode = i = 0; i < order0; ++i)
470 {
471 Blas::Dgemm('T', 'N', nquad2, order1, nquad1, 1.0,
472 tmp1.data() + i * nquad1 * nquad2, nquad1, base1.data(),
473 nquad1, 0.0, tmp2.data() + mode * nquad2, nquad2);
474 mode += order1;
475 }
476
477 // Inner product with respect to the '2' direction
478 mode = mode1 = cnt = 0;
479 for (i = 0; i < order0; ++i)
480 {
481 for (j = 0; j < order1; ++j, ++cnt)
482 {
483 int ijmax = max(i, j);
484
485 Blas::Dgemv('T', nquad2, order2 - ijmax, 1.0,
486 base2.data() + mode * nquad2, nquad2,
487 tmp2.data() + cnt * nquad2, 1, 0.0,
488 outarray.data() + mode1, 1);
489 mode += order2 - ijmax;
490 mode1 += order2 - ijmax;
491 }
492
493 // increment mode in case order1!=order2
494 for (j = order1; j < order2; ++j)
495 {
496 int ijmax = max(i, j);
497 mode += order2 - ijmax;
498 }
499 }
500
501 // fix for modified basis for top singular vertex component
502 // Already have evaluated (1+c)/2 (1-b)/2 (1-a)/2
504 {
505 // add in (1+c)/2 (1+b)/2 (1-a)/2 component
506 outarray[1] +=
507 Blas::Ddot(nquad2, base2.data() + nquad2, 1, &tmp2[nquad2], 1);
508
509 // add in (1+c)/2 (1-b)/2 (1+a)/2 component
510 outarray[1] += Blas::Ddot(nquad2, base2.data() + nquad2, 1,
511 &tmp2[nquad2 * order1], 1);
512
513 // add in (1+c)/2 (1+b)/2 (1+a)/2 component
514 outarray[1] += Blas::Ddot(nquad2, base2.data() + nquad2, 1,
515 &tmp2[nquad2 * order1 + nquad2], 1);
516 }
517}
518
520 const int dir, const Array<OneD, const NekDouble> &inarray,
521 Array<OneD, NekDouble> &outarray)
522{
523 StdPyrExp::v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
524}
525
526/**
527 * @param inarray Function evaluated at physical collocation
528 * points.
529 * @param outarray Inner product with respect to each basis
530 * function over the element.
531 */
533 const int dir, const Array<OneD, const NekDouble> &inarray,
534 Array<OneD, NekDouble> &outarray)
535{
536 int i;
537 int nquad0 = m_base[0]->GetNumPoints();
538 int nquad1 = m_base[1]->GetNumPoints();
539 int nquad2 = m_base[2]->GetNumPoints();
540 int nqtot = nquad0 * nquad1 * nquad2;
541
542 Array<OneD, NekDouble> gfac0(nquad0);
543 Array<OneD, NekDouble> gfac1(nquad1);
544 Array<OneD, NekDouble> gfac2(nquad2);
545 Array<OneD, NekDouble> tmp0(nqtot);
546
547 int order0 = m_base[0]->GetNumModes();
548 int order1 = m_base[1]->GetNumModes();
549
550 Array<OneD, NekDouble> wsp(nquad1 * nquad2 * order0 +
551 nquad2 * order0 * order1);
552
553 const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
554 const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
555 const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
556
557 // set up geometric factor: (1+z0)/2
558 for (i = 0; i < nquad0; ++i)
559 {
560 gfac0[i] = 0.5 * (1 + z0[i]);
561 }
562
563 // set up geometric factor: (1+z1)/2
564 for (i = 0; i < nquad1; ++i)
565 {
566 gfac1[i] = 0.5 * (1 + z1[i]);
567 }
568
569 // Set up geometric factor: 2/(1-z2)
570 for (i = 0; i < nquad2; ++i)
571 {
572 gfac2[i] = 2.0 / (1 - z2[i]);
573 }
574
575 // Derivative in first/second direction is always scaled as follows
576 const int nq01 = nquad0 * nquad1;
577 for (i = 0; i < nquad2; ++i)
578 {
579 Vmath::Smul(nq01, gfac2[i], &inarray[0] + i * nq01, 1,
580 &tmp0[0] + i * nq01, 1);
581 }
582
584
585 switch (dir)
586 {
587 case 0:
588 {
590 m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
591 m_base[2]->GetBdata(), tmp0, outarray, wsp, false, true, true);
592 break;
593 }
594 case 1:
595 {
597 m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
598 m_base[2]->GetBdata(), tmp0, outarray, wsp, true, false, true);
599 break;
600 }
601 case 2:
602 {
605
606 // Scale eta_1 derivative by gfac0
607 for (i = 0; i < nquad1 * nquad2; ++i)
608 {
609 Vmath::Vmul(nquad0, tmp0.data() + i * nquad0, 1, gfac0.data(),
610 1, tmp0.data() + i * nquad0, 1);
611 }
613 m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
614 m_base[2]->GetBdata(), tmp0, tmp3, wsp, false, true, true);
615
616 // Scale eta_2 derivative by gfac1*gfac2
617 for (i = 0; i < nquad2; ++i)
618 {
619 Vmath::Smul(nq01, gfac2[i], &inarray[0] + i * nq01, 1,
620 &tmp0[0] + i * nq01, 1);
621 }
622 for (i = 0; i < nquad1 * nquad2; ++i)
623 {
624 Vmath::Smul(nquad0, gfac1[i % nquad1], &tmp0[0] + i * nquad0, 1,
625 &tmp0[0] + i * nquad0, 1);
626 }
627
630 m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
631 m_base[2]->GetBdata(), tmp0, tmp4, wsp, true, false, true);
632
633 v_MultiplyByStdQuadratureMetric(inarray, tmp0);
635 m_base[0]->GetBdata(), m_base[1]->GetBdata(),
636 m_base[2]->GetDbdata(), tmp0, outarray, wsp, true, true, false);
637
638 Vmath::Vadd(m_ncoeffs, &tmp3[0], 1, &outarray[0], 1, &outarray[0],
639 1);
640 Vmath::Vadd(m_ncoeffs, &tmp4[0], 1, &outarray[0], 1, &outarray[0],
641 1);
642 break;
643 }
644 default:
645 {
646 ASSERTL1(false, "input dir is out of range");
647 break;
648 }
649 }
650}
651
652//---------------------------------------
653// Evaluation functions
654//---------------------------------------
655
658{
659 NekDouble d2 = 1.0 - xi[2];
660 if (fabs(d2) < NekConstants::kNekZeroTol)
661 {
662 if (d2 >= 0.)
663 {
665 }
666 else
667 {
669 }
670 }
671 eta[2] = xi[2]; // eta_z = xi_z
672 eta[1] = 2.0 * (1.0 + xi[1]) / d2 - 1.0;
673 eta[0] = 2.0 * (1.0 + xi[0]) / d2 - 1.0;
674}
675
678{
679 xi[0] = (1.0 + eta[0]) * (1.0 - eta[2]) * 0.5 - 1.0;
680 xi[1] = (1.0 + eta[1]) * (1.0 - eta[2]) * 0.5 - 1.0;
681 xi[2] = eta[2];
682}
683
687{
688 Array<OneD, const NekDouble> etaBar_x = m_base[0]->GetZ();
689 Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
690 Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
691 int Qx = GetNumPoints(0);
692 int Qy = GetNumPoints(1);
693 int Qz = GetNumPoints(2);
694
695 // Convert collapsed coordinates into cartesian coordinates: eta --> xi
696 for (int k = 0; k < Qz; ++k)
697 {
698 for (int j = 0; j < Qy; ++j)
699 {
700 for (int i = 0; i < Qx; ++i)
701 {
702 int s = i + Qx * (j + Qy * k);
703
704 xi_z[s] = eta_z[k];
705 xi_y[s] = (1.0 + eta_y[j]) * (1.0 - eta_z[k]) / 2.0 - 1.0;
706 xi_x[s] = (1.0 + etaBar_x[i]) * (1.0 - eta_z[k]) / 2.0 - 1.0;
707 }
708 }
709 }
710}
711
713 const Array<OneD, const NekDouble> &coords, int mode)
714{
716 LocCoordToLocCollapsed(coords, coll);
717
718 const int nm0 = m_base[0]->GetNumModes();
719 const int nm1 = m_base[1]->GetNumModes();
720 const int nm2 = m_base[2]->GetNumModes();
721
722 int mode0 = 0, mode1 = 0, mode2 = 0, cnt = 0;
723
724 bool found = false;
725 for (mode0 = 0; mode0 < nm0; ++mode0)
726 {
727 for (mode1 = 0; mode1 < nm1; ++mode1)
728 {
729 int maxpq = max(mode0, mode1);
730 for (mode2 = 0; mode2 < nm2 - maxpq; ++mode2, ++cnt)
731 {
732 if (cnt == mode)
733 {
734 found = true;
735 break;
736 }
737 }
738
739 if (found)
740 {
741 break;
742 }
743 }
744
745 if (found)
746 {
747 break;
748 }
749
750 for (int j = nm1; j < nm2; ++j)
751 {
752 int ijmax = max(mode0, j);
753 mode2 += nm2 - ijmax;
754 }
755 }
756
757 if (mode == 1 && m_base[0]->GetBasisType() == LibUtilities::eModified_A)
758 {
759 return StdExpansion::BaryEvaluateBasis<2>(coll[2], 1);
760 }
761 else
762 {
763 return StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
764 StdExpansion::BaryEvaluateBasis<1>(coll[1], mode1) *
765 StdExpansion::BaryEvaluateBasis<2>(coll[2], mode2);
766 }
767}
768
770 const Array<OneD, NekDouble> &coord,
771 const Array<OneD, const NekDouble> &inarray,
772 std::array<NekDouble, 3> &firstOrderDerivs)
773{
774 // Collapse coordinates
775 Array<OneD, NekDouble> coll(3, 0.0);
776 LocCoordToLocCollapsed(coord, coll);
777
778 // If near singularity do the old interpolation matrix method
779 if ((1 - coll[2]) < 1e-5)
780 {
781 int totPoints = GetTotPoints();
782 Array<OneD, NekDouble> EphysDeriv0(totPoints), EphysDeriv1(totPoints),
783 EphysDeriv2(totPoints);
784 PhysDeriv(inarray, EphysDeriv0, EphysDeriv1, EphysDeriv2);
785
787 I[0] = GetBase()[0]->GetI(coll);
788 I[1] = GetBase()[1]->GetI(coll + 1);
789 I[2] = GetBase()[2]->GetI(coll + 2);
790
791 firstOrderDerivs[0] = PhysEvaluate(I, EphysDeriv0);
792 firstOrderDerivs[1] = PhysEvaluate(I, EphysDeriv1);
793 firstOrderDerivs[2] = PhysEvaluate(I, EphysDeriv2);
794 return PhysEvaluate(I, inarray);
795 }
796
797 std::array<NekDouble, 3> interDeriv;
798 NekDouble val = StdExpansion3D::BaryTensorDeriv(coll, inarray, interDeriv);
799
800 NekDouble fac = 2.0 / (1.0 - coll[2]);
801
802 firstOrderDerivs[0] = fac * interDeriv[0];
803 firstOrderDerivs[1] = fac * interDeriv[1];
804 firstOrderDerivs[2] = ((1.0 + coll[0]) / (1.0 - coll[2])) * interDeriv[0] +
805 ((1.0 + coll[1]) / (1.0 - coll[2])) * interDeriv[1] +
806 interDeriv[2];
807
808 return val;
809}
810
811void StdPyrExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
812{
814 tmp[mode] = 1.0;
815 v_BwdTrans(tmp, outarray);
816}
817
818void StdPyrExp::v_GetTraceNumModes(const int fid, int &numModes0,
819 int &numModes1, Orientation faceOrient)
820{
821 int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
822 m_base[2]->GetNumModes()};
823 switch (fid)
824 {
825 // quad
826 case 0:
827 {
828 numModes0 = nummodes[0];
829 numModes1 = nummodes[1];
830 }
831 break;
832 case 1:
833 case 3:
834 {
835 numModes0 = nummodes[0];
836 numModes1 = nummodes[2];
837 }
838 break;
839 case 2:
840 case 4:
841 {
842 numModes0 = nummodes[1];
843 numModes1 = nummodes[2];
844 }
845 break;
846 }
847
848 if (faceOrient >= 9)
849 {
850 std::swap(numModes0, numModes1);
851 }
852}
853
854//---------------------------------------
855// Helper functions
856//---------------------------------------
857
859{
860 return 5;
861}
862
864{
865 return 8;
866}
867
869{
870 return 5;
871}
872
874{
876}
877
879{
882 "BasisType is not a boundary interior form");
885 "BasisType is not a boundary interior form");
888 "BasisType is not a boundary interior form");
889
890 int P = m_base[0]->GetNumModes();
891 int Q = m_base[1]->GetNumModes();
892 int R = m_base[2]->GetNumModes();
893
895}
896
898{
901 "BasisType is not a boundary interior form");
904 "BasisType is not a boundary interior form");
907 "BasisType is not a boundary interior form");
908
909 int P = m_base[0]->GetNumModes() - 1;
910 int Q = m_base[1]->GetNumModes() - 1;
911 int R = m_base[2]->GetNumModes() - 1;
912
913 return (P + 1) * (Q + 1) // 1 rect. face on base
914 + 2 * (R + 1) + P * (1 + 2 * R - P) // 2 tri. (P,R) faces
915 + 2 * (R + 1) + Q * (1 + 2 * R - Q); // 2 tri. (Q,R) faces
916}
917
918int StdPyrExp::v_GetTraceNcoeffs(const int i) const
919{
920 ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
921
922 if (i == 0)
923 {
924 return GetBasisNumModes(0) * GetBasisNumModes(1);
925 }
926 else if (i == 1 || i == 3)
927 {
928 int P = GetBasisNumModes(0) - 1, Q = GetBasisNumModes(2) - 1;
929 return Q + 1 + (P * (1 + 2 * Q - P)) / 2;
930 }
931 else
932 {
933 int P = GetBasisNumModes(1) - 1, Q = GetBasisNumModes(2) - 1;
934 return Q + 1 + (P * (1 + 2 * Q - P)) / 2;
935 }
936}
937
939{
940 ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
941
942 int P = m_base[0]->GetNumModes() - 1;
943 int Q = m_base[1]->GetNumModes() - 1;
944 int R = m_base[2]->GetNumModes() - 1;
945
946 if (i == 0)
947 {
948 return (P - 1) * (Q - 1);
949 }
950 else if (i == 1 || i == 3)
951 {
952 return (P - 1) * (2 * (R - 1) - (P - 1) - 1) / 2;
953 }
954 else
955 {
956 return (Q - 1) * (2 * (R - 1) - (Q - 1) - 1) / 2;
957 }
958}
959
960int StdPyrExp::v_GetTraceNumPoints(const int i) const
961{
962 ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
963
964 if (i == 0)
965 {
966 return m_base[0]->GetNumPoints() * m_base[1]->GetNumPoints();
967 }
968 else if (i == 1 || i == 3)
969 {
970 return m_base[0]->GetNumPoints() * m_base[2]->GetNumPoints();
971 }
972 else
973 {
974 return m_base[1]->GetNumPoints() * m_base[2]->GetNumPoints();
975 }
976}
977
978int StdPyrExp::v_GetEdgeNcoeffs(const int i) const
979{
980 ASSERTL2(i >= 0 && i <= 7, "edge id is out of range");
981
982 if (i == 0 || i == 2)
983 {
984 return GetBasisNumModes(0);
985 }
986 else if (i == 1 || i == 3)
987 {
988 return GetBasisNumModes(1);
989 }
990 else
991 {
992 return GetBasisNumModes(2);
993 }
994}
995
997 const int k,
998 bool UseGLL) const
999{
1000 ASSERTL2(i >= 0 && i <= 4, "face id is out of range");
1001 ASSERTL2(k >= 0 && k <= 1, "basis key id is out of range");
1002
1003 switch (i)
1004 {
1005 case 0:
1006 {
1008 m_base[k]->GetNumPoints(),
1009 m_base[k]->GetNumModes());
1010 }
1011 case 1:
1012 case 3:
1013 {
1015 k, m_base[2 * k]->GetBasisType(), m_base[2 * k]->GetNumPoints(),
1016 m_base[2 * k]->GetNumModes(), UseGLL);
1017 }
1018 case 2:
1019 case 4:
1020 {
1022 k, m_base[k + 1]->GetBasisType(), m_base[k + 1]->GetNumPoints(),
1023 m_base[k + 1]->GetNumModes(), UseGLL);
1024 }
1025 }
1026
1027 // Should never get here.
1029}
1030
1032 const std::vector<unsigned int> &nummodes, int &modes_offset)
1033{
1035 nummodes[modes_offset], nummodes[modes_offset + 1],
1036 nummodes[modes_offset + 2]);
1037
1038 modes_offset += 3;
1039 return nmodes;
1040}
1041
1042int StdPyrExp::v_GetVertexMap(int vId, bool useCoeffPacking)
1043{
1047 "Mapping not defined for this type of basis");
1048
1049 int l = 0;
1050
1051 if (useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
1052 {
1053 switch (vId)
1054 {
1055 case 0:
1056 l = GetMode(0, 0, 0);
1057 break;
1058 case 1:
1059 l = GetMode(0, 0, 1);
1060 break;
1061 case 2:
1062 l = GetMode(0, 1, 0);
1063 break;
1064 case 3:
1065 l = GetMode(1, 0, 0);
1066 break;
1067 case 4:
1068 l = GetMode(1, 1, 0);
1069 break;
1070 default:
1071 ASSERTL0(false, "local vertex id must be between 0 and 4");
1072 }
1073 }
1074 else
1075 {
1076 switch (vId)
1077 {
1078 case 0:
1079 l = GetMode(0, 0, 0);
1080 break;
1081 case 1:
1082 l = GetMode(1, 0, 0);
1083 break;
1084 case 2:
1085 l = GetMode(1, 1, 0);
1086 break;
1087 case 3:
1088 l = GetMode(0, 1, 0);
1089 break;
1090 case 4:
1091 l = GetMode(0, 0, 1);
1092 break;
1093 default:
1094 ASSERTL0(false, "local vertex id must be between 0 and 4");
1095 }
1096 }
1097
1098 return l;
1099}
1100
1102{
1105 "BasisType is not a boundary interior form");
1108 "BasisType is not a boundary interior form");
1111 "BasisType is not a boundary interior form");
1112
1113 int P = m_base[0]->GetNumModes() - 1, p;
1114 int Q = m_base[1]->GetNumModes() - 1, q;
1115 int R = m_base[2]->GetNumModes() - 1, r;
1116
1117 int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1118
1119 if (outarray.size() != nIntCoeffs)
1120 {
1121 outarray = Array<OneD, unsigned int>(nIntCoeffs);
1122 }
1123
1124 int idx = 0;
1125
1126 // Loop over all interior modes.
1127 for (p = 2; p <= P; ++p)
1128 {
1129 for (q = 2; q <= Q; ++q)
1130 {
1131 int maxpq = max(p, q);
1132 for (r = 1; r <= R - maxpq; ++r)
1133 {
1134 outarray[idx++] = GetMode(p, q, r);
1135 }
1136 }
1137 }
1138}
1139
1141{
1144 "BasisType is not a boundary interior form");
1147 "BasisType is not a boundary interior form");
1150 "BasisType is not a boundary interior form");
1151
1152 int P = m_base[0]->GetNumModes() - 1, p;
1153 int Q = m_base[1]->GetNumModes() - 1, q;
1154 int R = m_base[2]->GetNumModes() - 1, r;
1155 int idx = 0;
1156
1157 int nBnd = NumBndryCoeffs();
1158
1159 if (maparray.size() != nBnd)
1160 {
1161 maparray = Array<OneD, unsigned int>(nBnd);
1162 }
1163
1164 // Loop over all boundary modes (in ascending order).
1165 for (p = 0; p <= P; ++p)
1166 {
1167 // First two q-r planes are entirely boundary modes.
1168 if (p <= 1)
1169 {
1170 for (q = 0; q <= Q; ++q)
1171 {
1172 int maxpq = max(p, q);
1173 for (r = 0; r <= R - maxpq; ++r)
1174 {
1175 maparray[idx++] = GetMode(p, q, r);
1176 }
1177 }
1178 }
1179 else
1180 {
1181 // Remaining q-r planes contain boundary modes on the two
1182 // front and back sides and edges 0 2.
1183 for (q = 0; q <= Q; ++q)
1184 {
1185 if (q <= 1)
1186 {
1187 for (r = 0; r <= R - p; ++r)
1188 {
1189 maparray[idx++] = GetMode(p, q, r);
1190 }
1191 }
1192 else
1193 {
1194 maparray[idx++] = GetMode(p, q, 0);
1195 }
1196 }
1197 }
1198 }
1199}
1200
1201void StdPyrExp::v_GetTraceCoeffMap(const unsigned int fid,
1202 Array<OneD, unsigned int> &maparray)
1203{
1205 "Method only implemented if BasisType is identical"
1206 "in x and y directions");
1209 "Method only implemented for Modified_A BasisType"
1210 "(x and y direction) and ModifiedPyr_C BasisType (z "
1211 "direction)");
1212
1213 int p, q, r, P = 0, Q = 0, idx = 0;
1214
1215 int order0 = m_base[0]->GetNumModes();
1216 int order1 = m_base[1]->GetNumModes();
1217 int order2 = m_base[2]->GetNumModes();
1218
1219 switch (fid)
1220 {
1221 case 0:
1222 P = order0;
1223 Q = order1;
1224 break;
1225 case 1:
1226 case 3:
1227 P = order0;
1228 Q = order2;
1229 break;
1230 case 2:
1231 case 4:
1232 P = order1;
1233 Q = order2;
1234 break;
1235 default:
1236 ASSERTL0(false, "fid must be between 0 and 4");
1237 }
1238
1239 if (maparray.size() != P * Q)
1240 {
1241 maparray = Array<OneD, unsigned int>(P * Q);
1242 }
1243
1244 // Set up ordering inside each 2D face. Also for triangular faces,
1245 // populate signarray.
1246 switch (fid)
1247 {
1248 case 0: // Bottom quad
1249
1250 for (q = 0; q < Q; ++q)
1251 {
1252 for (p = 0; p < P; ++p)
1253 {
1254 maparray[q * P + p] = GetMode(p, q, 0);
1255 }
1256 }
1257 break;
1258
1259 case 1: // Front triangle
1260 for (p = 0; p < P; ++p)
1261 {
1262 for (r = 0; r < Q - p; ++r)
1263 {
1264 maparray[idx++] = GetMode(p, 0, r);
1265 }
1266 }
1267 break;
1268
1269 case 2: // Right triangle
1270 maparray[idx++] = GetMode(1, 0, 0);
1271 maparray[idx++] = GetMode(0, 0, 1);
1272 for (r = 1; r < Q - 1; ++r)
1273 {
1274 maparray[idx++] = GetMode(1, 0, r);
1275 }
1276
1277 for (q = 1; q < P; ++q)
1278 {
1279 for (r = 0; r < Q - q; ++r)
1280 {
1281 maparray[idx++] = GetMode(1, q, r);
1282 }
1283 }
1284 break;
1285
1286 case 3: // Rear triangle
1287 maparray[idx++] = GetMode(0, 1, 0);
1288 maparray[idx++] = GetMode(0, 0, 1);
1289 for (r = 1; r < Q - 1; ++r)
1290 {
1291 maparray[idx++] = GetMode(0, 1, r);
1292 }
1293
1294 for (p = 1; p < P; ++p)
1295 {
1296 for (r = 0; r < Q - p; ++r)
1297 {
1298 maparray[idx++] = GetMode(p, 1, r);
1299 }
1300 }
1301 break;
1302
1303 case 4: // Left triangle
1304 for (q = 0; q < P; ++q)
1305 {
1306 for (r = 0; r < Q - q; ++r)
1307 {
1308 maparray[idx++] = GetMode(0, q, r);
1309 }
1310 }
1311 break;
1312
1313 default:
1314 ASSERTL0(false, "Face to element map unavailable.");
1315 }
1316}
1317
1318void StdPyrExp::v_GetElmtTraceToTraceMap(const unsigned int fid,
1319 Array<OneD, unsigned int> &maparray,
1320 Array<OneD, int> &signarray,
1321 Orientation faceOrient, int P, int Q)
1322{
1324 "Method only implemented if BasisType is identical"
1325 "in x and y directions");
1328 "Method only implemented for Modified_A BasisType"
1329 "(x and y direction) and ModifiedPyr_C BasisType (z "
1330 "direction)");
1331
1332 int i, j, k, p, r, nFaceCoeffs;
1333 int nummodesA = 0, nummodesB = 0;
1334
1335 int order0 = m_base[0]->GetNumModes();
1336 int order1 = m_base[1]->GetNumModes();
1337 int order2 = m_base[2]->GetNumModes();
1338
1339 switch (fid)
1340 {
1341 case 0:
1342 nummodesA = order0;
1343 nummodesB = order1;
1344 break;
1345 case 1:
1346 case 3:
1347 nummodesA = order0;
1348 nummodesB = order2;
1349 break;
1350 case 2:
1351 case 4:
1352 nummodesA = order1;
1353 nummodesB = order2;
1354 break;
1355 default:
1356 ASSERTL0(false, "fid must be between 0 and 4");
1357 }
1358
1359 if (P == -1)
1360 {
1361 P = nummodesA;
1362 Q = nummodesB;
1363 nFaceCoeffs = GetTraceNcoeffs(fid);
1364 }
1365 else if (fid > 0)
1366 {
1367 nFaceCoeffs = P * (2 * Q - P + 1) / 2;
1368 }
1369 else
1370 {
1371 nFaceCoeffs = P * Q;
1372 }
1373
1374 // Allocate the map array and sign array; set sign array to ones (+)
1375 if (maparray.size() != nFaceCoeffs)
1376 {
1377 maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1378 }
1379
1380 if (signarray.size() != nFaceCoeffs)
1381 {
1382 signarray = Array<OneD, int>(nFaceCoeffs, 1);
1383 }
1384 else
1385 {
1386 fill(signarray.data(), signarray.data() + nFaceCoeffs, 1);
1387 }
1388
1389 // triangular faces
1390 if (fid > 0)
1391 {
1392 // zero signmap and set maparray to zero if elemental
1393 // modes are not as large as face modesl
1394 int idx = 0;
1395 int cnt = 0;
1396 int minPA = min(nummodesA, P);
1397 int minQB = min(nummodesB, Q);
1398
1399 for (j = 0; j < minPA; ++j)
1400 {
1401 // set maparray
1402 for (k = 0; k < minQB - j; ++k, ++cnt)
1403 {
1404 maparray[idx++] = cnt;
1405 }
1406
1407 cnt += nummodesB - minQB;
1408
1409 for (k = nummodesB - j; k < Q - j; ++k)
1410 {
1411 signarray[idx] = 0.0;
1412 maparray[idx++] = maparray[0];
1413 }
1414 }
1415#if 0 // not required?
1416
1417 for (j = minPA; j < nummodesA; ++j)
1418 {
1419 // set maparray
1420 for (k = 0; k < minQB-j; ++k, ++cnt)
1421 {
1422 maparray[idx++] = cnt;
1423 }
1424
1425 cnt += nummodesB-minQB;
1426
1427 for (k = nummodesB-j; k < Q-j; ++k)
1428 {
1429 signarray[idx] = 0.0;
1430 maparray[idx++] = maparray[0];
1431 }
1432 }
1433#endif
1434 for (j = nummodesA; j < P; ++j)
1435 {
1436 for (k = 0; k < Q - j; ++k)
1437 {
1438 signarray[idx] = 0.0;
1439 maparray[idx++] = maparray[0];
1440 }
1441 }
1442
1443 // Triangles only have one possible orientation (base
1444 // direction reversed); swap edge modes.
1445 if (faceOrient == eDir1BwdDir1_Dir2FwdDir2)
1446 {
1447 swap(maparray[0], maparray[Q]);
1448 for (i = 1; i < Q - 1; ++i)
1449 {
1450 swap(maparray[i + 1], maparray[Q + i]);
1451 }
1452
1453 idx = 0;
1454 for (p = 0; p < P; ++p)
1455 {
1456 for (r = 0; r < Q - p; ++r, idx++)
1457 {
1458 if (p > 1)
1459 {
1460 signarray[idx] = p % 2 ? -1 : 1;
1461 }
1462 }
1463 }
1464 }
1465 }
1466 else
1467 {
1468
1469 // Set up an array indexing for quads, since the ordering may need
1470 // to be transposed.
1471 Array<OneD, int> arrayindx(nFaceCoeffs, -1);
1472
1473 for (i = 0; i < Q; i++)
1474 {
1475 for (j = 0; j < P; j++)
1476 {
1477 if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1478 {
1479 arrayindx[i * P + j] = i * P + j;
1480 }
1481 else
1482 {
1483 arrayindx[i * P + j] = j * Q + i;
1484 }
1485 }
1486 }
1487
1488 // zero signmap and set maparray to zero if elemental
1489 // modes are not as large as face modesl
1490 for (j = 0; j < P; ++j)
1491 {
1492 // set up default maparray
1493 for (k = 0; k < Q; k++)
1494 {
1495 maparray[arrayindx[j + k * P]] = j + k * nummodesA;
1496 }
1497
1498 for (k = nummodesB; k < Q; ++k)
1499 {
1500 signarray[arrayindx[j + k * P]] = 0.0;
1501 maparray[arrayindx[j + k * P]] = maparray[0];
1502 }
1503 }
1504
1505 for (j = nummodesA; j < P; ++j)
1506 {
1507 for (k = 0; k < Q; ++k)
1508 {
1509 signarray[arrayindx[j + k * P]] = 0.0;
1510 maparray[arrayindx[j + k * P]] = maparray[0];
1511 }
1512 }
1513
1514 // The code below is exactly the same as that taken from
1515 // StdHexExp and reverses the 'b' and 'a' directions as
1516 // appropriate (1st and 2nd if statements respectively) in
1517 // quadrilateral faces.
1518 if (faceOrient == eDir1FwdDir1_Dir2BwdDir2 ||
1519 faceOrient == eDir1BwdDir1_Dir2BwdDir2 ||
1520 faceOrient == eDir1BwdDir2_Dir2FwdDir1 ||
1521 faceOrient == eDir1BwdDir2_Dir2BwdDir1)
1522 {
1523 if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1524 {
1525 for (i = 3; i < Q; i += 2)
1526 {
1527 for (j = 0; j < P; j++)
1528 {
1529 signarray[arrayindx[i * P + j]] *= -1;
1530 }
1531 }
1532
1533 for (i = 0; i < P; i++)
1534 {
1535 swap(maparray[i], maparray[i + P]);
1536 swap(signarray[i], signarray[i + P]);
1537 }
1538 }
1539 else
1540 {
1541 for (i = 0; i < Q; i++)
1542 {
1543 for (j = 3; j < P; j += 2)
1544 {
1545 signarray[arrayindx[i * P + j]] *= -1;
1546 }
1547 }
1548
1549 for (i = 0; i < Q; i++)
1550 {
1551 swap(maparray[i], maparray[i + Q]);
1552 swap(signarray[i], signarray[i + Q]);
1553 }
1554 }
1555 }
1556
1557 if (faceOrient == eDir1BwdDir1_Dir2FwdDir2 ||
1558 faceOrient == eDir1BwdDir1_Dir2BwdDir2 ||
1559 faceOrient == eDir1FwdDir2_Dir2BwdDir1 ||
1560 faceOrient == eDir1BwdDir2_Dir2BwdDir1)
1561 {
1562 if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1563 {
1564 for (i = 0; i < Q; i++)
1565 {
1566 for (j = 3; j < P; j += 2)
1567 {
1568 signarray[arrayindx[i * P + j]] *= -1;
1569 }
1570 }
1571
1572 for (i = 0; i < Q; i++)
1573 {
1574 swap(maparray[i * P], maparray[i * P + 1]);
1575 swap(signarray[i * P], signarray[i * P + 1]);
1576 }
1577 }
1578 else
1579 {
1580 for (i = 3; i < Q; i += 2)
1581 {
1582 for (j = 0; j < P; j++)
1583 {
1584 signarray[arrayindx[i * P + j]] *= -1;
1585 }
1586 }
1587
1588 for (i = 0; i < P; i++)
1589 {
1590 swap(maparray[i * Q], maparray[i * Q + 1]);
1591 swap(signarray[i * Q], signarray[i * Q + 1]);
1592 }
1593 }
1594 }
1595 }
1596}
1597
1599 const int eid, Array<OneD, unsigned int> &maparray,
1600 Array<OneD, int> &signarray, const Orientation edgeOrient)
1601{
1602 int i;
1603 bool signChange;
1604 const int P = m_base[0]->GetNumModes() - 1;
1605 const int Q = m_base[1]->GetNumModes() - 1;
1606 const int R = m_base[2]->GetNumModes() - 1;
1607 const int nEdgeIntCoeffs = v_GetEdgeNcoeffs(eid) - 2;
1608
1609 if (maparray.size() != nEdgeIntCoeffs)
1610 {
1611 maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1612 }
1613
1614 if (signarray.size() != nEdgeIntCoeffs)
1615 {
1616 signarray = Array<OneD, int>(nEdgeIntCoeffs, 1);
1617 }
1618 else
1619 {
1620 fill(signarray.data(), signarray.data() + nEdgeIntCoeffs, 1);
1621 }
1622
1623 // If edge is oriented backwards, change sign of modes which have
1624 // degree 2n+1, n >= 1.
1625 signChange = edgeOrient == eBackwards;
1626
1627 switch (eid)
1628 {
1629 case 0:
1630 for (i = 2; i <= P; ++i)
1631 {
1632 maparray[i - 2] = GetMode(i, 0, 0);
1633 }
1634 break;
1635
1636 case 1:
1637 for (i = 2; i <= Q; ++i)
1638 {
1639 maparray[i - 2] = GetMode(1, i, 0);
1640 }
1641 break;
1642 case 2:
1643 for (i = 2; i <= P; ++i)
1644 {
1645 maparray[i - 2] = GetMode(i, 1, 0);
1646 }
1647 break;
1648
1649 case 3:
1650 for (i = 2; i <= Q; ++i)
1651 {
1652 maparray[i - 2] = GetMode(0, i, 0);
1653 }
1654 break;
1655 case 4:
1656 for (i = 2; i <= R; ++i)
1657 {
1658 maparray[i - 2] = GetMode(0, 0, i);
1659 }
1660 break;
1661
1662 case 5:
1663 for (i = 1; i <= R - 1; ++i)
1664 {
1665 maparray[i - 1] = GetMode(1, 0, i);
1666 }
1667 break;
1668 case 6:
1669 for (i = 1; i <= R - 1; ++i)
1670 {
1671 maparray[i - 1] = GetMode(1, 1, i);
1672 }
1673 break;
1674
1675 case 7:
1676 for (i = 1; i <= R - 1; ++i)
1677 {
1678 maparray[i - 1] = GetMode(0, 1, i);
1679 }
1680 break;
1681 default:
1682 ASSERTL0(false, "Edge not defined.");
1683 break;
1684 }
1685
1686 if (signChange)
1687 {
1688 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1689 {
1690 signarray[i] = -1;
1691 }
1692 }
1693}
1694
1696 const int fid, Array<OneD, unsigned int> &maparray,
1697 Array<OneD, int> &signarray, const Orientation faceOrient)
1698{
1699 const int P = m_base[0]->GetNumModes() - 1;
1700 const int Q = m_base[1]->GetNumModes() - 1;
1701 const int R = m_base[2]->GetNumModes() - 1;
1702 const int nFaceIntCoeffs = v_GetTraceIntNcoeffs(fid);
1703 int p, q, r, idx = 0;
1704 int nummodesA = 0;
1705 int nummodesB = 0;
1706 int i, j;
1707
1708 if (maparray.size() != nFaceIntCoeffs)
1709 {
1710 maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1711 }
1712
1713 if (signarray.size() != nFaceIntCoeffs)
1714 {
1715 signarray = Array<OneD, int>(nFaceIntCoeffs, 1);
1716 }
1717 else
1718 {
1719 fill(signarray.data(), signarray.data() + nFaceIntCoeffs, 1);
1720 }
1721
1722 // Set up an array indexing for quad faces, since the ordering may
1723 // need to be transposed depending on orientation.
1724 Array<OneD, int> arrayindx(nFaceIntCoeffs);
1725 if (fid == 0)
1726 {
1727 nummodesA = P - 1;
1728 nummodesB = Q - 1;
1729
1730 for (i = 0; i < nummodesB; i++)
1731 {
1732 for (j = 0; j < nummodesA; j++)
1733 {
1734 if (faceOrient < 9)
1735 {
1736 arrayindx[i * nummodesA + j] = i * nummodesA + j;
1737 }
1738 else
1739 {
1740 arrayindx[i * nummodesA + j] = j * nummodesB + i;
1741 }
1742 }
1743 }
1744 }
1745
1746 switch (fid)
1747 {
1748 case 0: // Bottom quad
1749 for (q = 2; q <= Q; ++q)
1750 {
1751 for (p = 2; p <= P; ++p)
1752 {
1753 maparray[arrayindx[(q - 2) * nummodesA + (p - 2)]] =
1754 GetMode(p, q, 0);
1755 }
1756 }
1757 break;
1758 case 1: // Front triangle
1759 for (p = 2; p <= P; ++p)
1760 {
1761 for (r = 1; r <= R - p; ++r)
1762 {
1763 if ((int)faceOrient == 7)
1764 {
1765 signarray[idx] = p % 2 ? -1 : 1;
1766 }
1767 maparray[idx++] = GetMode(p, 0, r);
1768 }
1769 }
1770 break;
1771 case 2: // Right triangle
1772 for (q = 2; q <= Q; ++q)
1773 {
1774 for (r = 1; r <= R - q; ++r)
1775 {
1776 if ((int)faceOrient == 7)
1777 {
1778 signarray[idx] = q % 2 ? -1 : 1;
1779 }
1780 maparray[idx++] = GetMode(1, q, r);
1781 }
1782 }
1783 break;
1784
1785 case 3: // Rear triangle
1786 for (p = 2; p <= P; ++p)
1787 {
1788 for (r = 1; r <= R - p; ++r)
1789 {
1790 if ((int)faceOrient == 7)
1791 {
1792 signarray[idx] = p % 2 ? -1 : 1;
1793 }
1794 maparray[idx++] = GetMode(p, 1, r);
1795 }
1796 }
1797 break;
1798
1799 case 4: // Left triangle
1800 for (q = 2; q <= Q; ++q)
1801 {
1802 for (r = 1; r <= R - q; ++r)
1803 {
1804 if ((int)faceOrient == 7)
1805 {
1806 signarray[idx] = q % 2 ? -1 : 1;
1807 }
1808 maparray[idx++] = GetMode(0, q, r);
1809 }
1810 }
1811 break;
1812 default:
1813 ASSERTL0(false, "Face interior map not available.");
1814 }
1815
1816 // Triangular faces are processed in the above switch loop; for
1817 // remaining quad faces, set up orientation if necessary.
1818 if (fid > 0)
1819 {
1820 return;
1821 }
1822
1823 if (faceOrient == 6 || faceOrient == 8 || faceOrient == 11 ||
1824 faceOrient == 12)
1825 {
1826 if (faceOrient < 9)
1827 {
1828 for (i = 1; i < nummodesB; i += 2)
1829 {
1830 for (j = 0; j < nummodesA; j++)
1831 {
1832 signarray[arrayindx[i * nummodesA + j]] *= -1;
1833 }
1834 }
1835 }
1836 else
1837 {
1838 for (i = 0; i < nummodesB; i++)
1839 {
1840 for (j = 1; j < nummodesA; j += 2)
1841 {
1842 signarray[arrayindx[i * nummodesA + j]] *= -1;
1843 }
1844 }
1845 }
1846 }
1847
1848 if (faceOrient == 7 || faceOrient == 8 || faceOrient == 10 ||
1849 faceOrient == 12)
1850 {
1851 if (faceOrient < 9)
1852 {
1853 for (i = 0; i < nummodesB; i++)
1854 {
1855 for (j = 1; j < nummodesA; j += 2)
1856 {
1857 signarray[arrayindx[i * nummodesA + j]] *= -1;
1858 }
1859 }
1860 }
1861 else
1862 {
1863 for (i = 1; i < nummodesB; i += 2)
1864 {
1865 for (j = 0; j < nummodesA; j++)
1866 {
1867 signarray[arrayindx[i * nummodesA + j]] *= -1;
1868 }
1869 }
1870 }
1871 }
1872}
1873
1874//---------------------------------------
1875// Wrapper functions
1876//---------------------------------------
1877
1879{
1880 return CreateGeneralMatrix(mkey);
1881}
1882
1884{
1885 return v_GenMatrix(mkey);
1886}
1887
1888/**
1889 * @brief Compute the mode number in the expansion for a
1890 * particular tensorial combination.
1891 *
1892 * Modes are numbered with the r index travelling fastest,
1893 * followed by q and then p, and each q-r plane is of size
1894 *
1895 * (R+1-p)*(Q+1) - l(l+1)/2 where l = max(0,Q-p)
1896 *
1897 * For example, when P=2, Q=3 and R=4 the indexing inside each
1898 * q-r plane (with r increasing upwards and q to the right)
1899 * is:
1900 *
1901 * p = 0: p = 1: p = 2:
1902 * ----------------------------------
1903 * 4
1904 * 3 8 17 21
1905 * 2 7 11 16 20 24 29 32 35
1906 * 1 6 10 13 15 19 23 26 28 31 34 37
1907 * 0 5 9 12 14 18 22 25 27 30 33 36
1908 *
1909 * Note that in this element, we must have that \f$ P,Q \leq
1910 * R\f$.
1911 */
1912int StdPyrExp::GetMode(const int I, const int J, const int K)
1913{
1914 const int Q = m_base[1]->GetNumModes() - 1;
1915 const int R = m_base[2]->GetNumModes() - 1;
1916
1917 int i, l;
1918 int cnt = 0;
1919
1920 // Traverse to q-r plane number I
1921 for (i = 0; i < I; ++i)
1922 {
1923 // Size of triangle part
1924 l = max(0, Q - i);
1925
1926 // Size of rectangle part
1927 cnt += (R + 1 - i) * (Q + 1) - l * (l + 1) / 2;
1928 }
1929
1930 // Traverse to q column J (Pretend this is a face of width J)
1931 l = max(0, J - 1 - I);
1932 cnt += (R + 1 - I) * J - l * (l + 1) / 2;
1933
1934 // Traverse up stacks to K
1935 cnt += K;
1936
1937 return cnt;
1938}
1939
1941 const Array<OneD, const NekDouble> &inarray,
1942 Array<OneD, NekDouble> &outarray)
1943{
1944 int i, j;
1945
1946 int nquad0 = m_base[0]->GetNumPoints();
1947 int nquad1 = m_base[1]->GetNumPoints();
1948 int nquad2 = m_base[2]->GetNumPoints();
1949
1950 const Array<OneD, const NekDouble> &w0 = m_base[0]->GetW();
1951 const Array<OneD, const NekDouble> &w1 = m_base[1]->GetW();
1952 const Array<OneD, const NekDouble> &w2 = m_base[2]->GetW();
1953
1954 const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
1955
1956 // Multiply by integration constants in x-direction
1957 for (i = 0; i < nquad1 * nquad2; ++i)
1958 {
1959 Vmath::Vmul(nquad0, inarray.data() + i * nquad0, 1, w0.data(), 1,
1960 outarray.data() + i * nquad0, 1);
1961 }
1962
1963 // Multiply by integration constants in y-direction
1964 for (j = 0; j < nquad2; ++j)
1965 {
1966 for (i = 0; i < nquad1; ++i)
1967 {
1968 Blas::Dscal(nquad0, w1[i],
1969 &outarray[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1970 }
1971 }
1972
1973 // Multiply by integration constants in z-direction; need to
1974 // incorporate factor [(1-eta_3)/2]^2 into weights, but only if
1975 // using GLL quadrature points.
1976 switch (m_base[2]->GetPointsType())
1977 {
1978 // (2,0) Jacobi inner product.
1979 case LibUtilities::eGaussRadauMAlpha2Beta0:
1980 for (i = 0; i < nquad2; ++i)
1981 {
1982 Blas::Dscal(nquad0 * nquad1, 0.25 * w2[i],
1983 &outarray[0] + i * nquad0 * nquad1, 1);
1984 }
1985 break;
1986
1987 default:
1988 for (i = 0; i < nquad2; ++i)
1989 {
1990 Blas::Dscal(nquad0 * nquad1,
1991 0.25 * (1 - z2[i]) * (1 - z2[i]) * w2[i],
1992 &outarray[0] + i * nquad0 * nquad1, 1);
1993 }
1994 break;
1995 }
1996}
1997
1999 const StdMatrixKey &mkey)
2000{
2001 // Generate an orthonogal expansion
2002 int qa = m_base[0]->GetNumPoints();
2003 int qb = m_base[1]->GetNumPoints();
2004 int qc = m_base[2]->GetNumPoints();
2005 int nmodes_a = m_base[0]->GetNumModes();
2006 int nmodes_b = m_base[1]->GetNumModes();
2007 int nmodes_c = m_base[2]->GetNumModes();
2008 // Declare orthogonal basis.
2012
2016 StdPyrExp OrthoExp(Ba, Bb, Bc);
2017
2018 Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2019 int i, j, k, cnt = 0;
2020
2021 // project onto modal space.
2022 OrthoExp.FwdTrans(array, orthocoeffs);
2023
2025 {
2026 // Rodrigo's power kernel
2028 NekDouble SvvDiffCoeff =
2031
2032 for (i = 0; i < nmodes_a; ++i)
2033 {
2034 for (j = 0; j < nmodes_b; ++j)
2035 {
2036 int maxij = max(i, j);
2037 NekDouble fac1 = std::max(
2038 pow((1.0 * i) / (nmodes_a - 1), cutoff * nmodes_a),
2039 pow((1.0 * j) / (nmodes_b - 1), cutoff * nmodes_b));
2040
2041 for (k = 0; k < nmodes_c - maxij; ++k)
2042 {
2043 NekDouble fac =
2044 std::max(fac1, pow((1.0 * k) / (nmodes_c - 1),
2045 cutoff * nmodes_c));
2046
2047 orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2048 cnt++;
2049 }
2050 }
2051 }
2052 }
2053 else if (mkey.ConstFactorExists(
2054 eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
2055 {
2058
2059 int max_abc = max(nmodes_a - kSVVDGFiltermodesmin,
2060 nmodes_b - kSVVDGFiltermodesmin);
2061 max_abc = max(max_abc, nmodes_c - kSVVDGFiltermodesmin);
2062 // clamp max_abc
2063 max_abc = max(max_abc, 0);
2064 max_abc = min(max_abc, kSVVDGFiltermodesmax - kSVVDGFiltermodesmin);
2065
2066 for (i = 0; i < nmodes_a; ++i)
2067 {
2068 for (j = 0; j < nmodes_b; ++j)
2069 {
2070 int maxij = max(i, j);
2071
2072 for (k = 0; k < nmodes_c - maxij; ++k)
2073 {
2074 int maxijk = max(maxij, k);
2075 maxijk = min(maxijk, kSVVDGFiltermodesmax - 1);
2076
2077 orthocoeffs[cnt] *=
2078 SvvDiffCoeff * kSVVDGFilter[max_abc][maxijk];
2079 cnt++;
2080 }
2081 }
2082 }
2083 }
2084 else
2085 {
2086 // SVV filter paramaters (how much added diffusion relative
2087 // to physical one and fraction of modes from which you
2088 // start applying this added diffusion)
2089 //
2090 NekDouble SvvDiffCoeff =
2092 NekDouble SVVCutOff =
2094
2095 // Defining the cut of mode
2096 int cutoff_a = (int)(SVVCutOff * nmodes_a);
2097 int cutoff_b = (int)(SVVCutOff * nmodes_b);
2098 int cutoff_c = (int)(SVVCutOff * nmodes_c);
2099 // To avoid the fac[j] from blowing up
2100 NekDouble epsilon = 1;
2101
2102 int nmodes = min(min(nmodes_a, nmodes_b), nmodes_c);
2103 NekDouble cutoff = min(min(cutoff_a, cutoff_b), cutoff_c);
2104
2105 for (i = 0; i < nmodes_a; ++i) // P
2106 {
2107 for (j = 0; j < nmodes_b; ++j) // Q
2108 {
2109 int maxij = max(i, j);
2110 for (k = 0; k < nmodes_c - maxij; ++k) // R
2111 {
2112 if (j + k >= cutoff || i + k >= cutoff)
2113 {
2114 orthocoeffs[cnt] *=
2115 (SvvDiffCoeff *
2116 exp(-(i + k - nmodes) * (i + k - nmodes) /
2117 ((NekDouble)((i + k - cutoff + epsilon) *
2118 (i + k - cutoff + epsilon)))) *
2119 exp(-(j - nmodes) * (j - nmodes) /
2120 ((NekDouble)((j - cutoff + epsilon) *
2121 (j - cutoff + epsilon)))));
2122 }
2123 else
2124 {
2125 orthocoeffs[cnt] *= 0.0;
2126 }
2127 cnt++;
2128 }
2129 }
2130 }
2131 }
2132
2133 // backward transform to physical space
2134 OrthoExp.BwdTrans(orthocoeffs, array);
2135}
2136
2138 const Array<OneD, const NekDouble> &inarray,
2139 Array<OneD, NekDouble> &outarray)
2140{
2141 int nquad0 = m_base[0]->GetNumPoints();
2142 int nquad1 = m_base[1]->GetNumPoints();
2143 int nquad2 = m_base[2]->GetNumPoints();
2144 int nqtot = nquad0 * nquad1 * nquad2;
2145 int nmodes0 = m_base[0]->GetNumModes();
2146 int nmodes1 = m_base[1]->GetNumModes();
2147 int nmodes2 = m_base[2]->GetNumModes();
2148 int numMax = nmodes0;
2149
2151 Array<OneD, NekDouble> coeff_tmp1(m_ncoeffs, 0.0);
2152 Array<OneD, NekDouble> phys_tmp(nqtot, 0.0);
2153 Array<OneD, NekDouble> tmp, tmp2, tmp3, tmp4;
2154
2155 const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
2156 const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
2157 const LibUtilities::PointsKey Pkey2 = m_base[2]->GetPointsKey();
2158
2159 LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A, nmodes0, Pkey0);
2160 LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_A, nmodes1, Pkey1);
2161 LibUtilities::BasisKey bortho2(LibUtilities::eOrthoPyr_C, nmodes2, Pkey2);
2162
2163 int cnt = 0;
2164 int u = 0;
2165 int i = 0;
2167
2169 bortho0, bortho1, bortho2);
2170
2171 BwdTrans(inarray, phys_tmp);
2172 OrthoPyrExp->FwdTrans(phys_tmp, coeff);
2173
2174 // filtering
2175 for (u = 0; u < numMin; ++u)
2176 {
2177 for (i = 0; i < numMin; ++i)
2178 {
2179
2180 int maxui = max(u, i);
2181 Vmath::Vcopy(numMin - maxui, tmp = coeff + cnt, 1,
2182 tmp2 = coeff_tmp1 + cnt, 1);
2183 cnt += nmodes2 - maxui;
2184 }
2185
2186 for (i = numMin; i < nmodes1; ++i)
2187 {
2188 int maxui = max(u, i);
2189 cnt += numMax - maxui;
2190 }
2191 }
2192
2193 OrthoPyrExp->BwdTrans(coeff_tmp1, phys_tmp);
2194 StdPyrExp::FwdTrans(phys_tmp, outarray);
2195}
2196
2197} // 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
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:131
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 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
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
int GetTraceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th trace.
Definition: StdExpansion.h:261
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
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
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:124
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:133
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdPyrExp.cpp:532
DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey) override
Definition: StdPyrExp.cpp:1878
int v_NumBndryCoeffs() const override
Definition: StdPyrExp.cpp:878
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) final
Definition: StdPyrExp.cpp:712
void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdPyrExp.cpp:2137
int v_NumDGBndryCoeffs() const override
Definition: StdPyrExp.cpp:897
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: StdPyrExp.cpp:358
void v_GetElmtTraceToTraceMap(const unsigned int fid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation faceOrient, int P, int Q) override
Definition: StdPyrExp.cpp:1318
int v_GetTraceNcoeffs(const int i) const override
Definition: StdPyrExp.cpp:918
void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi) override
Definition: StdPyrExp.cpp:676
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey) override
Definition: StdPyrExp.cpp:1998
int v_GetEdgeNcoeffs(const int i) const override
Definition: StdPyrExp.cpp:978
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: StdPyrExp.cpp:264
int v_GetNtraces() const override
Definition: StdPyrExp.cpp:868
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdPyrExp.cpp:519
void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray) override
Definition: StdPyrExp.cpp:811
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
Definition: StdPyrExp.cpp:407
int v_GetNedges() const override
Definition: StdPyrExp.cpp:863
int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset) override
Definition: StdPyrExp.cpp:1031
int GetMode(int I, int J, int K)
Compute the mode number in the expansion for a particular tensorial combination.
Definition: StdPyrExp.cpp:1912
int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false) override
Definition: StdPyrExp.cpp:1042
int v_GetTraceNumPoints(const int i) const override
Definition: StdPyrExp.cpp:960
void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Backward transformation is evaluated at the quadrature points.
Definition: StdPyrExp.cpp:228
void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray) override
Definition: StdPyrExp.cpp:1140
void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Inner product of inarray over region with respect to the expansion basis m_base[0]->GetBdata(),...
Definition: StdPyrExp.cpp:393
void v_GetInteriorMap(Array< OneD, unsigned int > &outarray) override
Definition: StdPyrExp.cpp:1101
void v_GetEdgeInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
Definition: StdPyrExp.cpp:1598
void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta) override
Definition: StdPyrExp.cpp:656
DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey) override
Definition: StdPyrExp.cpp:1883
int v_GetNverts() const override
Definition: StdPyrExp.cpp:858
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdPyrExp.cpp:1940
void v_GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
Definition: StdPyrExp.cpp:1695
void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdPyrExp.cpp:247
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: StdPyrExp.cpp:437
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: StdPyrExp.cpp:188
LibUtilities::ShapeType v_DetShapeType() const override
Definition: StdPyrExp.cpp:873
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: StdPyrExp.cpp:83
void v_GetTraceNumModes(const int fid, int &numModes0, int &numModes1, Orientation faceOrient=eDir1FwdDir1_Dir2FwdDir2) override
Definition: StdPyrExp.cpp:818
void v_GetTraceCoeffMap(const unsigned int fid, Array< OneD, unsigned int > &maparray) override
Definition: StdPyrExp.cpp:1201
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
Definition: StdPyrExp.cpp:769
void v_GetCoords(Array< OneD, NekDouble > &xi_x, Array< OneD, NekDouble > &xi_y, Array< OneD, NekDouble > &xi_z) override
Definition: StdPyrExp.cpp:684
const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int k, bool UseGLL=false) const override
Definition: StdPyrExp.cpp:996
int v_GetTraceIntNcoeffs(const int i) const override
Definition: StdPyrExp.cpp:938
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 getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:232
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:259
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
@ P
Monomial polynomials .
Definition: BasisType.h:62
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:42
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:56
@ eModifiedPyr_C
Principle Modified Functions.
Definition: BasisType.h:53
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
@ eOrthoPyr_C
Principle Orthogonal Functions .
Definition: BasisType.h:51
static const NekDouble kNekZeroTol
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
std::shared_ptr< StdPyrExp > StdPyrExpSharedPtr
Definition: StdPyrExp.h:214
const int kSVVDGFiltermodesmax
Definition: StdRegions.hpp:501
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
const NekDouble kSVVDGFilter[9][11]
Definition: StdRegions.hpp:503
std::vector< double > q(NPUPPER *NPUPPER)
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 Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
STL namespace.