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