Nektar++
StdHexExp.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: StdHexExp.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: Heaxhedral methods
32//
33///////////////////////////////////////////////////////////////////////////////
34
36
37#ifdef max
38#undef max
39#endif
40
41using namespace std;
42
43namespace Nektar::StdRegions
44{
46 const LibUtilities::BasisKey &Bb,
47 const LibUtilities::BasisKey &Bc)
48 : StdExpansion(Ba.GetNumModes() * Bb.GetNumModes() * Bc.GetNumModes(), 3,
49 Ba, Bb, Bc),
50 StdExpansion3D(Ba.GetNumModes() * Bb.GetNumModes() * Bc.GetNumModes(), Ba,
51 Bb, Bc)
52{
53}
54
56{
60 (m_base[0]->GetBasisType() == LibUtilities::eGLL_Lagrange &&
61 m_base[1]->GetBasisType() == LibUtilities::eGLL_Lagrange &&
62 m_base[1]->GetBasisType() == LibUtilities::eGLL_Lagrange);
63}
64
65///////////////////////////////
66/// Differentiation Methods
67///////////////////////////////
68/**
69 * For Hexahedral region can use the PhysTensorDeriv function defined
70 * under StdExpansion. Following tenserproduct:
71 */
76{
77 StdExpansion3D::PhysTensorDeriv(inarray, out_d0, out_d1, out_d2);
78}
79
80/**
81 * @param dir Direction in which to compute derivative.
82 * Valid values are 0, 1, 2.
83 * @param inarray Input array.
84 * @param outarray Output array.
85 */
86void StdHexExp::v_PhysDeriv(const int dir,
87 const Array<OneD, const NekDouble> &inarray,
88 Array<OneD, NekDouble> &outarray)
89{
90 switch (dir)
91 {
92 case 0:
93 {
94 PhysDeriv(inarray, outarray, NullNekDouble1DArray,
96 }
97 break;
98 case 1:
99 {
100 PhysDeriv(inarray, NullNekDouble1DArray, outarray,
102 }
103 break;
104 case 2:
105 {
107 outarray);
108 }
109 break;
110 default:
111 {
112 ASSERTL1(false, "input dir is out of range");
113 }
114 break;
115 }
116}
117
122{
123 StdHexExp::v_PhysDeriv(inarray, out_d0, out_d1, out_d2);
124}
125
126void StdHexExp::v_StdPhysDeriv(const int dir,
127 const Array<OneD, const NekDouble> &inarray,
128 Array<OneD, NekDouble> &outarray)
129{
130 StdHexExp::v_PhysDeriv(dir, inarray, outarray);
131}
132
133/**
134 * Backward transformation is three dimensional tensorial expansion
135 * \f$ u (\xi_{1i}, \xi_{2j}, \xi_{3k})
136 * = \sum_{p=0}^{Q_x} \psi_p^a (\xi_{1i})
137 * \lbrace { \sum_{q=0}^{Q_y} \psi_{q}^a (\xi_{2j})
138 * \lbrace { \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{r}^a (\xi_{3k})
139 * \rbrace}
140 * \rbrace}. \f$
141 * And sumfactorizing step of the form is as:\\
142 * \f$ f_{r} (\xi_{3k})
143 * = \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{r}^a (\xi_{3k}),\\
144 * g_{p} (\xi_{2j}, \xi_{3k})
145 * = \sum_{r=0}^{Q_y} \psi_{p}^a (\xi_{2j}) f_{r} (\xi_{3k}),\\
146 * u(\xi_{1i}, \xi_{2j}, \xi_{3k})
147 * = \sum_{p=0}^{Q_x} \psi_{p}^a (\xi_{1i}) g_{p} (\xi_{2j}, \xi_{3k}).
148 * \f$
149 *
150 * @param inarray ?
151 * @param outarray ?
152 */
154 Array<OneD, NekDouble> &outarray)
155{
158 "Basis[1] is not a general tensor type");
159
162 "Basis[2] is not a general tensor type");
163
164 if (m_base[0]->Collocation() && m_base[1]->Collocation() &&
165 m_base[2]->Collocation())
166 {
168 m_base[2]->GetNumPoints(),
169 inarray, 1, outarray, 1);
170 }
171 else
172 {
173 StdHexExp::BwdTrans_SumFac(inarray, outarray);
174 }
175}
176
177/**
178 *
179 */
181 Array<OneD, NekDouble> &outarray)
182{
184 m_base[0]->GetNumPoints() * m_base[2]->GetNumModes() *
185 (m_base[1]->GetNumModes() + m_base[1]->GetNumPoints())); // FIX THIS
186
187 BwdTrans_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(),
188 m_base[2]->GetBdata(), inarray, outarray, wsp, true,
189 true, true);
190}
191
192/**
193 * @param base0 x-dirn basis matrix
194 * @param base1 y-dirn basis matrix
195 * @param base2 z-dirn basis matrix
196 * @param inarray Input vector of modes.
197 * @param outarray Output vector of physical space data.
198 * @param wsp Workspace of size Q_x*P_z*(P_y+Q_y)
199 * @param doCheckCollDir0 Check for collocation of basis.
200 * @param doCheckCollDir1 Check for collocation of basis.
201 * @param doCheckCollDir2 Check for collocation of basis.
202 * @todo Account for some directions being collocated. See
203 * StdQuadExp as an example.
204 */
206 const Array<OneD, const NekDouble> &base0,
207 const Array<OneD, const NekDouble> &base1,
208 const Array<OneD, const NekDouble> &base2,
209 const Array<OneD, const NekDouble> &inarray,
211 bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
212{
213 int nquad0 = m_base[0]->GetNumPoints();
214 int nquad1 = m_base[1]->GetNumPoints();
215 int nquad2 = m_base[2]->GetNumPoints();
216 int nmodes0 = m_base[0]->GetNumModes();
217 int nmodes1 = m_base[1]->GetNumModes();
218 int nmodes2 = m_base[2]->GetNumModes();
219
220 // Check if using collocation, if requested.
221 bool colldir0 = doCheckCollDir0 ? (m_base[0]->Collocation()) : false;
222 bool colldir1 = doCheckCollDir1 ? (m_base[1]->Collocation()) : false;
223 bool colldir2 = doCheckCollDir2 ? (m_base[2]->Collocation()) : false;
224
225 // If collocation in all directions, Physical values at quadrature
226 // points is just a copy of the modes.
227 if (colldir0 && colldir1 && colldir2)
228 {
229 Vmath::Vcopy(m_ncoeffs, inarray.data(), 1, outarray.data(), 1);
230 }
231 else
232 {
233 // Check sufficiently large workspace.
234 ASSERTL1(wsp.size() >= nquad0 * nmodes2 * (nmodes1 + nquad1),
235 "Workspace size is not sufficient");
236
237 // Assign second half of workspace for 2nd DGEMM operation.
238 Array<OneD, NekDouble> wsp2 = wsp + nquad0 * nmodes1 * nmodes2;
239
240 // BwdTrans in each direction using DGEMM
241 Blas::Dgemm('T', 'T', nmodes1 * nmodes2, nquad0, nmodes0, 1.0,
242 &inarray[0], nmodes0, base0.data(), nquad0, 0.0, &wsp[0],
243 nmodes1 * nmodes2);
244 Blas::Dgemm('T', 'T', nquad0 * nmodes2, nquad1, nmodes1, 1.0, &wsp[0],
245 nmodes1, base1.data(), nquad1, 0.0, &wsp2[0],
246 nquad0 * nmodes2);
247 Blas::Dgemm('T', 'T', nquad0 * nquad1, nquad2, nmodes2, 1.0, &wsp2[0],
248 nmodes2, base2.data(), nquad2, 0.0, &outarray[0],
249 nquad0 * nquad1);
250 }
251}
252
253/**
254 * Solves the system
255 * \f$ \mathbf{B^{\top}WB\hat{u}}=\mathbf{B^{\top}Wu^{\delta}} \f$
256 *
257 * @param inarray array of physical quadrature points to be
258 * transformed, \f$ \mathbf{u^{\delta}} \f$.
259 * @param outarray array of expansion coefficients,
260 * \f$ \mathbf{\hat{u}} \f$.
261 */
263 Array<OneD, NekDouble> &outarray)
264{
265 // If using collocation expansion, coefficients match physical
266 // data points so just do a direct copy.
267 if ((m_base[0]->Collocation()) && (m_base[1]->Collocation()) &&
268 (m_base[2]->Collocation()))
269 {
270 Vmath::Vcopy(GetNcoeffs(), &inarray[0], 1, &outarray[0], 1);
271 }
272 else
273 {
274 // Compute B^TWu
275 IProductWRTBase(inarray, outarray);
276
277 // get Mass matrix inverse
278 StdMatrixKey masskey(eInvMass, DetShapeType(), *this);
279 DNekMatSharedPtr matsys = GetStdMatrix(masskey);
280
281 // copy inarray in case inarray == outarray
282 DNekVec in(m_ncoeffs, outarray);
283 DNekVec out(m_ncoeffs, outarray, eWrapper);
284
285 // Solve for coefficients.
286 out = (*matsys) * in;
287 }
288}
289
290/**
291 * \f$
292 * \begin{array}{rcl}
293 * I_{pqr} = (\phi_{pqr}, u)_{\delta} & = &
294 * \sum_{i=0}^{nq_0} \sum_{j=0}^{nq_1} \sum_{k=0}^{nq_2}
295 * \psi_{p}^{a}(\xi_{1i}) \psi_{q}^{a}(\xi_{2j}) \psi_{r}^{a}(\xi_{3k})
296 * w_i w_j w_k u(\xi_{1,i} \xi_{2,j} \xi_{3,k})
297 *
298 * J_{i,j,k}\\ & = & \sum_{i=0}^{nq_0} \psi_p^a(\xi_{1,i})
299 * \sum_{j=0}^{nq_1} \psi_{q}^a(\xi_{2,j})
300 * \sum_{k=0}^{nq_2} \psi_{r}^a
301 * u(\xi_{1i},\xi_{2j},\xi_{3k}) J_{i,j,k}
302 * \end{array} \f$ \n
303 * where
304 * \f$ \phi_{pqr} (\xi_1 , \xi_2 , \xi_3)
305 * = \psi_p^a( \xi_1) \psi_{q}^a(\xi_2) \psi_{r}^a(\xi_3) \f$ \n
306 * which can be implemented as \n
307 * \f$f_{r} (\xi_{3k})
308 * = \sum_{k=0}^{nq_3} \psi_{r}^a u(\xi_{1i},\xi_{2j}, \xi_{3k})
309 * J_{i,j,k} = {\bf B_3 U} \f$ \n
310 * \f$ g_{q} (\xi_{3k})
311 * = \sum_{j=0}^{nq_1} \psi_{q}^a(\xi_{2j}) f_{r}(\xi_{3k})
312 * = {\bf B_2 F} \f$ \n
313 * \f$ (\phi_{pqr}, u)_{\delta}
314 * = \sum_{k=0}^{nq_0} \psi_{p}^a (\xi_{3k}) g_{q} (\xi_{3k})
315 * = {\bf B_1 G} \f$
316 *
317 * @param inarray ?
318 * @param outarray ?
319 */
321 Array<OneD, NekDouble> &outarray)
322{
323 if (m_base[0]->Collocation() && m_base[1]->Collocation() &&
324 m_base[2]->Collocation())
325 {
326 MultiplyByQuadratureMetric(inarray, outarray);
327 }
328 else
329 {
330 StdHexExp::v_IProductWRTBase_SumFac(inarray, outarray);
331 }
332}
333
334/**
335 * Implementation of the sum-factorization inner product operation.
336 */
338 const Array<OneD, const NekDouble> &inarray,
339 Array<OneD, NekDouble> &outarray, bool multiplybyweights)
340{
341 int nquad0 = m_base[0]->GetNumPoints();
342 int nquad1 = m_base[1]->GetNumPoints();
343 int nquad2 = m_base[2]->GetNumPoints();
344 int order0 = m_base[0]->GetNumModes();
345 int order1 = m_base[1]->GetNumModes();
346
347 Array<OneD, NekDouble> wsp(nquad0 * nquad1 * (nquad2 + order0) +
348 order0 * order1 * nquad2);
349
350 if (multiplybyweights)
351 {
352 Array<OneD, NekDouble> tmp(inarray.size());
353 MultiplyByQuadratureMetric(inarray, tmp);
354
356 m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
357 tmp, outarray, wsp, true, true, true);
358 }
359 else
360 {
362 m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
363 inarray, outarray, wsp, true, true, true);
364 }
365}
366
367/**
368 * Implementation of the sum-factorisation inner product operation.
369 * @todo Implement cases where only some directions are collocated.
370 */
372 const Array<OneD, const NekDouble> &base0,
373 const Array<OneD, const NekDouble> &base1,
374 const Array<OneD, const NekDouble> &base2,
375 const Array<OneD, const NekDouble> &inarray,
377 bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
378{
379 int nquad0 = m_base[0]->GetNumPoints();
380 int nquad1 = m_base[1]->GetNumPoints();
381 int nquad2 = m_base[2]->GetNumPoints();
382 int nmodes0 = m_base[0]->GetNumModes();
383 int nmodes1 = m_base[1]->GetNumModes();
384 int nmodes2 = m_base[2]->GetNumModes();
385
386 bool colldir0 = doCheckCollDir0 ? (m_base[0]->Collocation()) : false;
387 bool colldir1 = doCheckCollDir1 ? (m_base[1]->Collocation()) : false;
388 bool colldir2 = doCheckCollDir2 ? (m_base[2]->Collocation()) : false;
389
390 if (colldir0 && colldir1 && colldir2)
391 {
392 Vmath::Vcopy(m_ncoeffs, inarray.data(), 1, outarray.data(), 1);
393 }
394 else
395 {
396 ASSERTL1(wsp.size() >= nmodes0 * nquad2 * (nquad1 + nmodes1),
397 "Insufficient workspace size");
398
399 Array<OneD, NekDouble> tmp0 = wsp;
400 Array<OneD, NekDouble> tmp1 = wsp + nmodes0 * nquad1 * nquad2;
401
402 if (colldir0)
403 {
404 // reshuffle data for next operation.
405 for (int n = 0; n < nmodes0; ++n)
406 {
407 Vmath::Vcopy(nquad1 * nquad2, inarray.data() + n, nquad0,
408 tmp0.data() + nquad1 * nquad2 * n, 1);
409 }
410 }
411 else
412 {
413 Blas::Dgemm('T', 'N', nquad1 * nquad2, nmodes0, nquad0, 1.0,
414 inarray.data(), nquad0, base0.data(), nquad0, 0.0,
415 tmp0.data(), nquad1 * nquad2);
416 }
417
418 if (colldir1)
419 {
420 // reshuffle data for next operation.
421 for (int n = 0; n < nmodes1; ++n)
422 {
423 Vmath::Vcopy(nquad2 * nmodes0, tmp0.data() + n, nquad1,
424 tmp1.data() + nquad2 * nmodes0 * n, 1);
425 }
426 }
427 else
428 {
429 Blas::Dgemm('T', 'N', nquad2 * nmodes0, nmodes1, nquad1, 1.0,
430 tmp0.data(), nquad1, base1.data(), nquad1, 0.0,
431 tmp1.data(), nquad2 * nmodes0);
432 }
433
434 if (colldir2)
435 {
436 // reshuffle data for next operation.
437 for (int n = 0; n < nmodes2; ++n)
438 {
439 Vmath::Vcopy(nmodes0 * nmodes1, tmp1.data() + n, nquad2,
440 outarray.data() + nmodes0 * nmodes1 * n, 1);
441 }
442 }
443 else
444 {
445 Blas::Dgemm('T', 'N', nmodes0 * nmodes1, nmodes2, nquad2, 1.0,
446 tmp1.data(), nquad2, base2.data(), nquad2, 0.0,
447 outarray.data(), nmodes0 * nmodes1);
448 }
449 }
450}
451
453 const int dir, const Array<OneD, const NekDouble> &inarray,
454 Array<OneD, NekDouble> &outarray)
455{
456 StdHexExp::IProductWRTDerivBase_SumFac(dir, inarray, outarray);
457}
458
460 const int dir, const Array<OneD, const NekDouble> &inarray,
461 Array<OneD, NekDouble> &outarray)
462{
463 ASSERTL0((dir == 0) || (dir == 1) || (dir == 2),
464 "input dir is out of range");
465
466 int nquad1 = m_base[1]->GetNumPoints();
467 int nquad2 = m_base[2]->GetNumPoints();
468 int order0 = m_base[0]->GetNumModes();
469 int order1 = m_base[1]->GetNumModes();
470
471 // If outarray > inarray then no need for temporary storage.
472 Array<OneD, NekDouble> tmp = outarray;
473 if (outarray.size() < inarray.size())
474 {
475 tmp = Array<OneD, NekDouble>(inarray.size());
476 }
477
478 // Need workspace for sumfackernel though
479 Array<OneD, NekDouble> wsp(order0 * nquad2 * (nquad1 + order1));
480
481 // multiply by integration constants
482 MultiplyByQuadratureMetric(inarray, tmp);
483
484 // perform sum-factorisation
485 switch (dir)
486 {
487 case 0:
489 m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
490 m_base[2]->GetBdata(), tmp, outarray, wsp, false, true, true);
491 break;
492 case 1:
494 m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
495 m_base[2]->GetBdata(), tmp, outarray, wsp, true, false, true);
496 break;
497 case 2:
499 m_base[0]->GetBdata(), m_base[1]->GetBdata(),
500 m_base[2]->GetDbdata(), tmp, outarray, wsp, true, true, false);
501 break;
502 }
503}
504
507{
508 eta[0] = xi[0];
509 eta[1] = xi[1];
510 eta[2] = xi[2];
511}
512
515{
516 xi[0] = eta[0];
517 xi[1] = eta[1];
518 xi[2] = eta[2];
519}
520
521/**
522 * @note for hexahedral expansions _base[0] (i.e. p) modes run fastest.
523 */
524void StdHexExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
525{
526 int nquad0 = m_base[0]->GetNumPoints();
527 int nquad1 = m_base[1]->GetNumPoints();
528 int nquad2 = m_base[2]->GetNumPoints();
529
530 Array<OneD, const NekDouble> base0 = m_base[0]->GetBdata();
531 Array<OneD, const NekDouble> base1 = m_base[1]->GetBdata();
532 Array<OneD, const NekDouble> base2 = m_base[2]->GetBdata();
533
534 int btmp0 = m_base[0]->GetNumModes();
535 int btmp1 = m_base[1]->GetNumModes();
536 int mode2 = mode / (btmp0 * btmp1);
537 int mode1 = (mode - mode2 * btmp0 * btmp1) / btmp0;
538 int mode0 = (mode - mode2 * btmp0 * btmp1) % btmp0;
539
540 ASSERTL2(mode == mode2 * btmp0 * btmp1 + mode1 * btmp0 + mode0,
541 "Mode lookup failed.");
542 ASSERTL2(mode < m_ncoeffs,
543 "Calling argument mode is larger than total expansion "
544 "order");
545
546 for (int i = 0; i < nquad1 * nquad2; ++i)
547 {
548 Vmath::Vcopy(nquad0, (NekDouble *)(base0.data() + mode0 * nquad0), 1,
549 &outarray[0] + i * nquad0, 1);
550 }
551
552 for (int j = 0; j < nquad2; ++j)
553 {
554 for (int i = 0; i < nquad0; ++i)
555 {
556 Vmath::Vmul(nquad1, (NekDouble *)(base1.data() + mode1 * nquad1), 1,
557 &outarray[0] + i + j * nquad0 * nquad1, nquad0,
558 &outarray[0] + i + j * nquad0 * nquad1, nquad0);
559 }
560 }
561
562 for (int i = 0; i < nquad2; i++)
563 {
564 Blas::Dscal(nquad0 * nquad1, base2[mode2 * nquad2 + i],
565 &outarray[0] + i * nquad0 * nquad1, 1);
566 }
567}
568
570 const Array<OneD, const NekDouble> &coords, int mode)
571{
572 ASSERTL2(coords[0] > -1 - NekConstants::kNekZeroTol, "coord[0] < -1");
573 ASSERTL2(coords[0] < 1 + NekConstants::kNekZeroTol, "coord[0] > 1");
574 ASSERTL2(coords[1] > -1 - NekConstants::kNekZeroTol, "coord[1] < -1");
575 ASSERTL2(coords[1] < 1 + NekConstants::kNekZeroTol, "coord[1] > 1");
576 ASSERTL2(coords[2] > -1 - NekConstants::kNekZeroTol, "coord[2] < -1");
577 ASSERTL2(coords[2] < 1 + NekConstants::kNekZeroTol, "coord[2] > 1");
578
579 const int nm0 = m_base[0]->GetNumModes();
580 const int nm1 = m_base[1]->GetNumModes();
581 const int mode2 = mode / (nm0 * nm1);
582 const int mode1 = (mode - mode2 * nm0 * nm1) / nm0;
583 const int mode0 = (mode - mode2 * nm0 * nm1) % nm0;
584
585 return StdExpansion::BaryEvaluateBasis<0>(coords[0], mode0) *
586 StdExpansion::BaryEvaluateBasis<1>(coords[1], mode1) *
587 StdExpansion::BaryEvaluateBasis<2>(coords[2], mode2);
588}
589
591{
592 return 8;
593}
594
596{
597 return 12;
598}
599
601{
602 return 6;
603}
604
606{
608}
609
611{
614 "BasisType is not a boundary interior form");
617 "BasisType is not a boundary interior form");
620 "BasisType is not a boundary interior form");
621
622 int nmodes0 = m_base[0]->GetNumModes();
623 int nmodes1 = m_base[1]->GetNumModes();
624 int nmodes2 = m_base[2]->GetNumModes();
625
626 return (2 * (nmodes0 * nmodes1 + nmodes0 * nmodes2 + nmodes1 * nmodes2) -
627 4 * (nmodes0 + nmodes1 + nmodes2) + 8);
628}
629
631{
634 "BasisType is not a boundary interior form");
637 "BasisType is not a boundary interior form");
640 "BasisType is not a boundary interior form");
641
642 int nmodes0 = m_base[0]->GetNumModes();
643 int nmodes1 = m_base[1]->GetNumModes();
644 int nmodes2 = m_base[2]->GetNumModes();
645
646 return 2 * (nmodes0 * nmodes1 + nmodes0 * nmodes2 + nmodes1 * nmodes2);
647}
648
649int StdHexExp::v_GetTraceNcoeffs(const int i) const
650{
651 ASSERTL2((i >= 0) && (i <= 5), "face id is out of range");
652 if ((i == 0) || (i == 5))
653 {
654 return GetBasisNumModes(0) * GetBasisNumModes(1);
655 }
656 else if ((i == 1) || (i == 3))
657 {
658 return GetBasisNumModes(0) * GetBasisNumModes(2);
659 }
660 else
661 {
662 return GetBasisNumModes(1) * GetBasisNumModes(2);
663 }
664}
665
667{
668 ASSERTL2((i >= 0) && (i <= 5), "face id is out of range");
669 if ((i == 0) || (i == 5))
670 {
671 return (GetBasisNumModes(0) - 2) * (GetBasisNumModes(1) - 2);
672 }
673 else if ((i == 1) || (i == 3))
674 {
675 return (GetBasisNumModes(0) - 2) * (GetBasisNumModes(2) - 2);
676 }
677 else
678 {
679 return (GetBasisNumModes(1) - 2) * (GetBasisNumModes(2) - 2);
680 }
681}
682
683int StdHexExp::v_GetTraceNumPoints(const int i) const
684{
685 ASSERTL2(i >= 0 && i <= 5, "face id is out of range");
686
687 if (i == 0 || i == 5)
688 {
689 return m_base[0]->GetNumPoints() * m_base[1]->GetNumPoints();
690 }
691 else if (i == 1 || i == 3)
692 {
693 return m_base[0]->GetNumPoints() * m_base[2]->GetNumPoints();
694 }
695 else
696 {
697 return m_base[1]->GetNumPoints() * m_base[2]->GetNumPoints();
698 }
699}
700
702 const int j) const
703{
704 ASSERTL2(i >= 0 && i <= 5, "face id is out of range");
705 ASSERTL2(j == 0 || j == 1, "face direction is out of range");
706
707 if (i == 0 || i == 5)
708 {
709 return m_base[j]->GetPointsKey();
710 }
711 else if (i == 1 || i == 3)
712 {
713 return m_base[2 * j]->GetPointsKey();
714 }
715 else
716 {
717 return m_base[j + 1]->GetPointsKey();
718 }
719}
720
722 const std::vector<unsigned int> &nummodes, int &modes_offset)
723{
724 int nmodes = nummodes[modes_offset] * nummodes[modes_offset + 1] *
725 nummodes[modes_offset + 2];
726 modes_offset += 3;
727
728 return nmodes;
729}
730
732 const int i, const int k, [[maybe_unused]] bool UseGLL) const
733{
734 ASSERTL2(i >= 0 && i <= 5, "face id is out of range");
735 ASSERTL2(k >= 0 && k <= 1, "basis key id is out of range");
736
737 int dir = k;
738 switch (i)
739 {
740 case 0:
741 case 5:
742 dir = k;
743 break;
744 case 1:
745 case 3:
746 dir = 2 * k;
747 break;
748 case 2:
749 case 4:
750 dir = k + 1;
751 break;
752 }
753
755 m_base[dir]->GetNumPoints(),
756 m_base[dir]->GetNumModes());
757}
758
762{
763 Array<OneD, const NekDouble> eta_x = m_base[0]->GetZ();
764 Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
765 Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
766 int Qx = GetNumPoints(0);
767 int Qy = GetNumPoints(1);
768 int Qz = GetNumPoints(2);
769
770 // Convert collapsed coordinates into cartesian coordinates:
771 // eta --> xi
772 for (int k = 0; k < Qz; ++k)
773 {
774 for (int j = 0; j < Qy; ++j)
775 {
776 for (int i = 0; i < Qx; ++i)
777 {
778 int s = i + Qx * (j + Qy * k);
779 xi_x[s] = eta_x[i];
780 xi_y[s] = eta_y[j];
781 xi_z[s] = eta_z[k];
782 }
783 }
784 }
785}
786
787void StdHexExp::v_GetTraceNumModes(const int fid, int &numModes0,
788 int &numModes1, Orientation faceOrient)
789{
790 int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
791 m_base[2]->GetNumModes()};
792 switch (fid)
793 {
794 case 0:
795 case 5:
796 {
797 numModes0 = nummodes[0];
798 numModes1 = nummodes[1];
799 }
800 break;
801 case 1:
802 case 3:
803 {
804 numModes0 = nummodes[0];
805 numModes1 = nummodes[2];
806 }
807 break;
808 case 2:
809 case 4:
810 {
811 numModes0 = nummodes[1];
812 numModes1 = nummodes[2];
813 }
814 break;
815 default:
816 {
817 ASSERTL0(false, "fid out of range");
818 }
819 break;
820 }
821
822 if (faceOrient >= eDir1FwdDir2_Dir2FwdDir1)
823 {
824 std::swap(numModes0, numModes1);
825 }
826}
827
828/**
829 * Expansions in each of the three dimensions must be of type
830 * LibUtilities#eModified_A or LibUtilities#eGLL_Lagrange.
831 *
832 * @param localVertexId ID of vertex (0..7)
833 * @returns Position of vertex in local numbering scheme.
834 */
835int StdHexExp::v_GetVertexMap(const int localVertexId, bool useCoeffPacking)
836{
839 "BasisType is not a boundary interior form");
842 "BasisType is not a boundary interior form");
845 "BasisType is not a boundary interior form");
846
847 ASSERTL1((localVertexId >= 0) && (localVertexId < 8),
848 "local vertex id must be between 0 and 7");
849
850 int p = 0;
851 int q = 0;
852 int r = 0;
853
854 // Retrieve the number of modes in each dimension.
855 int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
856 m_base[2]->GetNumModes()};
857
858 if (useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
859 {
860 if (localVertexId > 3)
861 {
863 {
864 r = nummodes[2] - 1;
865 }
866 else
867 {
868 r = 1;
869 }
870 }
871
872 switch (localVertexId % 4)
873 {
874 case 0:
875 break;
876 case 1:
877 {
879 {
880 p = nummodes[0] - 1;
881 }
882 else
883 {
884 p = 1;
885 }
886 }
887 break;
888 case 2:
889 {
891 {
892 q = nummodes[1] - 1;
893 }
894 else
895 {
896 q = 1;
897 }
898 }
899 break;
900 case 3:
901 {
903 {
904 p = nummodes[0] - 1;
905 q = nummodes[1] - 1;
906 }
907 else
908 {
909 p = 1;
910 q = 1;
911 }
912 }
913 break;
914 }
915 }
916 else
917 {
918 // Right face (vertices 1,2,5,6)
919 if ((localVertexId % 4) % 3 > 0)
920 {
922 {
923 p = nummodes[0] - 1;
924 }
925 else
926 {
927 p = 1;
928 }
929 }
930 // Back face (vertices 2,3,6,7)
931 if (localVertexId % 4 > 1)
932 {
934 {
935 q = nummodes[1] - 1;
936 }
937 else
938 {
939 q = 1;
940 }
941 }
942
943 // Top face (vertices 4,5,6,7)
944 if (localVertexId > 3)
945 {
947 {
948 r = nummodes[2] - 1;
949 }
950 else
951 {
952 r = 1;
953 }
954 }
955 }
956 // Compute the local number.
957 return r * nummodes[0] * nummodes[1] + q * nummodes[0] + p;
958}
959
960/**
961 * @param outarray Storage area for computed map.
962 */
964{
967 "BasisType is not a boundary interior form");
970 "BasisType is not a boundary interior form");
973 "BasisType is not a boundary interior form");
974
975 int i;
976 int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
977 m_base[2]->GetNumModes()};
978
979 int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
980
981 if (outarray.size() != nIntCoeffs)
982 {
983 outarray = Array<OneD, unsigned int>(nIntCoeffs);
984 }
985
986 const LibUtilities::BasisType Btype[3] = {GetBasisType(0), GetBasisType(1),
987 GetBasisType(2)};
988
989 int p, q, r;
990 int cnt = 0;
991
992 int IntIdx[3][2];
993
994 for (i = 0; i < 3; i++)
995 {
996 if (Btype[i] == LibUtilities::eModified_A)
997 {
998 IntIdx[i][0] = 2;
999 IntIdx[i][1] = nummodes[i];
1000 }
1001 else
1002 {
1003 IntIdx[i][0] = 1;
1004 IntIdx[i][1] = nummodes[i] - 1;
1005 }
1006 }
1007
1008 for (r = IntIdx[2][0]; r < IntIdx[2][1]; r++)
1009 {
1010 for (q = IntIdx[1][0]; q < IntIdx[1][1]; q++)
1011 {
1012 for (p = IntIdx[0][0]; p < IntIdx[0][1]; p++)
1013 {
1014 outarray[cnt++] =
1015 r * nummodes[0] * nummodes[1] + q * nummodes[0] + p;
1016 }
1017 }
1018 }
1019}
1020
1021/**
1022 * @param outarray Storage for computed map.
1023 */
1025{
1028 "BasisType is not a boundary interior form");
1031 "BasisType is not a boundary interior form");
1034 "BasisType is not a boundary interior form");
1035
1036 int i;
1037 int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
1038 m_base[2]->GetNumModes()};
1039
1040 int nBndCoeffs = NumBndryCoeffs();
1041
1042 if (outarray.size() != nBndCoeffs)
1043 {
1044 outarray = Array<OneD, unsigned int>(nBndCoeffs);
1045 }
1046
1047 const LibUtilities::BasisType Btype[3] = {GetBasisType(0), GetBasisType(1),
1048 GetBasisType(2)};
1049
1050 int p, q, r;
1051 int cnt = 0;
1052
1053 int BndIdx[3][2];
1054 int IntIdx[3][2];
1055
1056 for (i = 0; i < 3; i++)
1057 {
1058 BndIdx[i][0] = 0;
1059
1060 if (Btype[i] == LibUtilities::eModified_A)
1061 {
1062 BndIdx[i][1] = 1;
1063 IntIdx[i][0] = 2;
1064 IntIdx[i][1] = nummodes[i];
1065 }
1066 else
1067 {
1068 BndIdx[i][1] = nummodes[i] - 1;
1069 IntIdx[i][0] = 1;
1070 IntIdx[i][1] = nummodes[i] - 1;
1071 }
1072 }
1073
1074 for (i = 0; i < 2; i++)
1075 {
1076 r = BndIdx[2][i];
1077 for (q = 0; q < nummodes[1]; q++)
1078 {
1079 for (p = 0; p < nummodes[0]; p++)
1080 {
1081 outarray[cnt++] =
1082 r * nummodes[0] * nummodes[1] + q * nummodes[0] + p;
1083 }
1084 }
1085 }
1086
1087 for (r = IntIdx[2][0]; r < IntIdx[2][1]; r++)
1088 {
1089 for (i = 0; i < 2; i++)
1090 {
1091 q = BndIdx[1][i];
1092 for (p = 0; p < nummodes[0]; p++)
1093 {
1094 outarray[cnt++] =
1095 r * nummodes[0] * nummodes[1] + q * nummodes[0] + p;
1096 }
1097 }
1098
1099 for (q = IntIdx[1][0]; q < IntIdx[1][1]; q++)
1100 {
1101 for (i = 0; i < 2; i++)
1102 {
1103 p = BndIdx[0][i];
1104 outarray[cnt++] =
1105 r * nummodes[0] * nummodes[1] + q * nummodes[0] + p;
1106 }
1107 }
1108 }
1109
1110 sort(outarray.data(), outarray.data() + nBndCoeffs);
1111}
1112
1114 const Array<OneD, NekDouble> &coord,
1115 const Array<OneD, const NekDouble> &inarray,
1116 std::array<NekDouble, 3> &firstOrderDerivs)
1117{
1118 return BaryTensorDeriv(coord, inarray, firstOrderDerivs);
1119}
1120
1121/**
1122 * Only for basis type Modified_A or GLL_LAGRANGE in all directions.
1123 */
1124void StdHexExp::v_GetTraceCoeffMap(const unsigned int fid,
1125 Array<OneD, unsigned int> &maparray)
1126{
1127 int i, j;
1128 int nummodesA = 0, nummodesB = 0;
1129
1131 GetBasisType(0) == GetBasisType(2),
1132 "Method only implemented if BasisType is indentical in "
1133 "all directions");
1136 "Method only implemented for Modified_A or "
1137 "GLL_Lagrange BasisType");
1138
1139 const int nummodes0 = m_base[0]->GetNumModes();
1140 const int nummodes1 = m_base[1]->GetNumModes();
1141 const int nummodes2 = m_base[2]->GetNumModes();
1142
1143 switch (fid)
1144 {
1145 case 0:
1146 case 5:
1147 nummodesA = nummodes0;
1148 nummodesB = nummodes1;
1149 break;
1150 case 1:
1151 case 3:
1152 nummodesA = nummodes0;
1153 nummodesB = nummodes2;
1154 break;
1155 case 2:
1156 case 4:
1157 nummodesA = nummodes1;
1158 nummodesB = nummodes2;
1159 break;
1160 default:
1161 ASSERTL0(false, "fid must be between 0 and 5");
1162 }
1163
1164 int nFaceCoeffs = nummodesA * nummodesB;
1165
1166 if (maparray.size() != nFaceCoeffs)
1167 {
1168 maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1169 }
1170
1171 bool modified = (GetBasisType(0) == LibUtilities::eModified_A);
1172
1173 int offset = 0;
1174 int jump1 = 1;
1175 int jump2 = 1;
1176
1177 switch (fid)
1178 {
1179 case 5:
1180 {
1181 if (modified)
1182 {
1183 offset = nummodes0 * nummodes1;
1184 }
1185 else
1186 {
1187 offset = (nummodes2 - 1) * nummodes0 * nummodes1;
1188 jump1 = nummodes0;
1189 }
1190 }
1191 /* Falls through. */
1192 case 0:
1193 {
1194 jump1 = nummodes0;
1195 break;
1196 }
1197 case 3:
1198 {
1199 if (modified)
1200 {
1201 offset = nummodes0;
1202 }
1203 else
1204 {
1205 offset = nummodes0 * (nummodes1 - 1);
1206 jump1 = nummodes0 * nummodes1;
1207 }
1208 }
1209 /* Falls through. */
1210 case 1:
1211 {
1212 jump1 = nummodes0 * nummodes1;
1213 break;
1214 }
1215 case 2:
1216 {
1217 if (modified)
1218 {
1219 offset = 1;
1220 }
1221 else
1222 {
1223 offset = nummodes0 - 1;
1224 jump1 = nummodes0 * nummodes1;
1225 jump2 = nummodes0;
1226 }
1227 }
1228 /* Falls through. */
1229 case 4:
1230 {
1231 jump1 = nummodes0 * nummodes1;
1232 jump2 = nummodes0;
1233 break;
1234 }
1235 default:
1236 ASSERTL0(false, "fid must be between 0 and 5");
1237 }
1238
1239 for (i = 0; i < nummodesB; i++)
1240 {
1241 for (j = 0; j < nummodesA; j++)
1242 {
1243 maparray[i * nummodesA + j] = i * jump1 + j * jump2 + offset;
1244 }
1245 }
1246}
1247
1248void StdHexExp::v_GetElmtTraceToTraceMap(const unsigned int fid,
1249 Array<OneD, unsigned int> &maparray,
1250 Array<OneD, int> &signarray,
1251 Orientation faceOrient, int P, int Q)
1252{
1253 int i, j;
1254 int nummodesA = 0, nummodesB = 0;
1255
1257 GetBasisType(0) == GetBasisType(2),
1258 "Method only implemented if BasisType is indentical in "
1259 "all directions");
1262 "Method only implemented for Modified_A or "
1263 "GLL_Lagrange BasisType");
1264
1265 const int nummodes0 = m_base[0]->GetNumModes();
1266 const int nummodes1 = m_base[1]->GetNumModes();
1267 const int nummodes2 = m_base[2]->GetNumModes();
1268
1269 switch (fid)
1270 {
1271 case 0:
1272 case 5:
1273 nummodesA = nummodes0;
1274 nummodesB = nummodes1;
1275 break;
1276 case 1:
1277 case 3:
1278 nummodesA = nummodes0;
1279 nummodesB = nummodes2;
1280 break;
1281 case 2:
1282 case 4:
1283 nummodesA = nummodes1;
1284 nummodesB = nummodes2;
1285 break;
1286 default:
1287 ASSERTL0(false, "fid must be between 0 and 5");
1288 }
1289
1290 if (P == -1)
1291 {
1292 P = nummodesA;
1293 Q = nummodesB;
1294 }
1295
1296 bool modified = (GetBasisType(0) == LibUtilities::eModified_A);
1297
1298 // check that
1299 if (modified == false)
1300 {
1301 ASSERTL1((P == nummodesA) || (Q == nummodesB),
1302 "Different trace space face dimention "
1303 "and element face dimention not possible for "
1304 "GLL-Lagrange bases");
1305 }
1306
1307 int nFaceCoeffs = P * Q;
1308
1309 if (maparray.size() != nFaceCoeffs)
1310 {
1311 maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1312 }
1313
1314 // fill default mapping as increasing index
1315 for (i = 0; i < nFaceCoeffs; ++i)
1316 {
1317 maparray[i] = i;
1318 }
1319
1320 if (signarray.size() != nFaceCoeffs)
1321 {
1322 signarray = Array<OneD, int>(nFaceCoeffs, 1);
1323 }
1324 else
1325 {
1326 fill(signarray.data(), signarray.data() + nFaceCoeffs, 1);
1327 }
1328
1329 // setup indexing to manage transpose directions
1330 Array<OneD, int> arrayindx(nFaceCoeffs);
1331 for (i = 0; i < Q; i++)
1332 {
1333 for (j = 0; j < P; j++)
1334 {
1335 if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1336 {
1337 arrayindx[i * P + j] = i * P + j;
1338 }
1339 else
1340 {
1341 arrayindx[i * P + j] = j * Q + i;
1342 }
1343 }
1344 }
1345
1346 // zero signmap and set maparray to zero if elemental
1347 // modes are not as large as face models
1348 for (i = 0; i < nummodesB; i++)
1349 {
1350 for (j = nummodesA; j < P; j++)
1351 {
1352 signarray[arrayindx[i * P + j]] = 0.0;
1353 maparray[arrayindx[i * P + j]] = maparray[0];
1354 }
1355 }
1356
1357 for (i = nummodesB; i < Q; i++)
1358 {
1359 for (j = 0; j < P; j++)
1360 {
1361 signarray[arrayindx[i * P + j]] = 0.0;
1362 maparray[arrayindx[i * P + j]] = maparray[0];
1363 }
1364 }
1365
1366 // zero signmap and set maparray to zero entry if
1367 // elemental modes are not as large as face modesl
1368 for (i = 0; i < Q; i++)
1369 {
1370 // fill values into map array of trace size
1371 // for element face index
1372 for (j = 0; j < P; j++)
1373 {
1374 maparray[arrayindx[i * P + j]] = i * nummodesA + j;
1375 }
1376
1377 // zero values if P > numModesA
1378 for (j = nummodesA; j < P; j++)
1379 {
1380 signarray[arrayindx[i * P + j]] = 0.0;
1381 maparray[arrayindx[i * P + j]] = maparray[0];
1382 }
1383 }
1384
1385 // zero values if Q > numModesB
1386 for (i = nummodesB; i < Q; i++)
1387 {
1388 for (j = 0; j < P; j++)
1389 {
1390 signarray[arrayindx[i * P + j]] = 0.0;
1391 maparray[arrayindx[i * P + j]] = maparray[0];
1392 }
1393 }
1394
1395 // Now reorientate indices accordign to orientation
1396 if ((faceOrient == eDir1FwdDir1_Dir2BwdDir2) ||
1397 (faceOrient == eDir1BwdDir1_Dir2BwdDir2) ||
1398 (faceOrient == eDir1BwdDir2_Dir2FwdDir1) ||
1399 (faceOrient == eDir1BwdDir2_Dir2BwdDir1))
1400 {
1401 if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1402 {
1403 if (modified)
1404 {
1405 for (i = 3; i < Q; i += 2)
1406 {
1407 for (j = 0; j < P; j++)
1408 {
1409 signarray[arrayindx[i * P + j]] *= -1;
1410 }
1411 }
1412
1413 for (i = 0; i < P; i++)
1414 {
1415 swap(maparray[i], maparray[i + P]);
1416 swap(signarray[i], signarray[i + P]);
1417 }
1418 }
1419 else
1420 {
1421 for (i = 0; i < P; i++)
1422 {
1423 for (j = 0; j < Q / 2; j++)
1424 {
1425 swap(maparray[i + j * P],
1426 maparray[i + P * Q - P - j * P]);
1427 swap(signarray[i + j * P],
1428 signarray[i + P * Q - P - j * P]);
1429 }
1430 }
1431 }
1432 }
1433 else
1434 {
1435 if (modified)
1436 {
1437 for (i = 0; i < Q; i++)
1438 {
1439 for (j = 3; j < P; j += 2)
1440 {
1441 signarray[arrayindx[i * P + j]] *= -1;
1442 }
1443 }
1444
1445 for (i = 0; i < Q; i++)
1446 {
1447 swap(maparray[i], maparray[i + Q]);
1448 swap(signarray[i], signarray[i + Q]);
1449 }
1450 }
1451 else
1452 {
1453 for (i = 0; i < P; i++)
1454 {
1455 for (j = 0; j < Q / 2; j++)
1456 {
1457 swap(maparray[i * Q + j], maparray[i * Q + Q - 1 - j]);
1458 swap(signarray[i * Q + j],
1459 signarray[i * Q + Q - 1 - j]);
1460 }
1461 }
1462 }
1463 }
1464 }
1465
1466 if ((faceOrient == eDir1BwdDir1_Dir2FwdDir2) ||
1467 (faceOrient == eDir1BwdDir1_Dir2BwdDir2) ||
1468 (faceOrient == eDir1FwdDir2_Dir2BwdDir1) ||
1469 (faceOrient == eDir1BwdDir2_Dir2BwdDir1))
1470 {
1471 if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1472 {
1473 if (modified)
1474 {
1475 for (i = 0; i < Q; i++)
1476 {
1477 for (j = 3; j < P; j += 2)
1478 {
1479 signarray[arrayindx[i * P + j]] *= -1;
1480 }
1481 }
1482
1483 for (i = 0; i < Q; i++)
1484 {
1485 swap(maparray[i * P], maparray[i * P + 1]);
1486 swap(signarray[i * P], signarray[i * P + 1]);
1487 }
1488 }
1489 else
1490 {
1491 for (i = 0; i < Q; i++)
1492 {
1493 for (j = 0; j < P / 2; j++)
1494 {
1495 swap(maparray[i * P + j], maparray[i * P + P - 1 - j]);
1496 swap(signarray[i * P + j],
1497 signarray[i * P + P - 1 - j]);
1498 }
1499 }
1500 }
1501 }
1502 else
1503 {
1504 if (modified)
1505 {
1506 for (i = 3; i < Q; i += 2)
1507 {
1508 for (j = 0; j < P; j++)
1509 {
1510 signarray[arrayindx[i * P + j]] *= -1;
1511 }
1512 }
1513
1514 for (i = 0; i < P; i++)
1515 {
1516 swap(maparray[i * Q], maparray[i * Q + 1]);
1517 swap(signarray[i * Q], signarray[i * Q + 1]);
1518 }
1519 }
1520 else
1521 {
1522 for (i = 0; i < Q; i++)
1523 {
1524 for (j = 0; j < P / 2; j++)
1525 {
1526 swap(maparray[i + j * Q],
1527 maparray[i + P * Q - Q - j * Q]);
1528 swap(signarray[i + j * Q],
1529 signarray[i + P * Q - Q - j * Q]);
1530 }
1531 }
1532 }
1533 }
1534 }
1535}
1536
1537/**
1538 * @param eid The edge to compute the numbering for.
1539 * @param edgeOrient Orientation of the edge.
1540 * @param maparray Storage for computed mapping array.
1541 * @param signarray ?
1542 */
1544 const int eid, Array<OneD, unsigned int> &maparray,
1545 Array<OneD, int> &signarray, const Orientation edgeOrient)
1546{
1549 "BasisType is not a boundary interior form");
1552 "BasisType is not a boundary interior form");
1555 "BasisType is not a boundary interior form");
1556
1557 ASSERTL1((eid >= 0) && (eid < 12),
1558 "local edge id must be between 0 and 11");
1559
1560 int nEdgeIntCoeffs = GetEdgeNcoeffs(eid) - 2;
1561
1562 if (maparray.size() != nEdgeIntCoeffs)
1563 {
1564 maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1565 }
1566
1567 if (signarray.size() != nEdgeIntCoeffs)
1568 {
1569 signarray = Array<OneD, int>(nEdgeIntCoeffs, 1);
1570 }
1571 else
1572 {
1573 fill(signarray.data(), signarray.data() + nEdgeIntCoeffs, 1);
1574 }
1575
1576 int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
1577 m_base[2]->GetNumModes()};
1578
1579 const LibUtilities::BasisType bType[3] = {GetBasisType(0), GetBasisType(1),
1580 GetBasisType(2)};
1581
1582 bool reverseOrdering = false;
1583 bool signChange = false;
1584
1585 int IdxRange[3][2] = {{0, 0}, {0, 0}, {0, 0}};
1586
1587 switch (eid)
1588 {
1589 case 0:
1590 case 1:
1591 case 2:
1592 case 3:
1593 {
1594 IdxRange[2][0] = 0;
1595 IdxRange[2][1] = 1;
1596 }
1597 break;
1598 case 8:
1599 case 9:
1600 case 10:
1601 case 11:
1602 {
1603 if (bType[2] == LibUtilities::eGLL_Lagrange)
1604 {
1605 IdxRange[2][0] = nummodes[2] - 1;
1606 IdxRange[2][1] = nummodes[2];
1607 }
1608 else
1609 {
1610 IdxRange[2][0] = 1;
1611 IdxRange[2][1] = 2;
1612 }
1613 }
1614 break;
1615 case 4:
1616 case 5:
1617 case 6:
1618 case 7:
1619 {
1620 if (bType[2] == LibUtilities::eGLL_Lagrange)
1621 {
1622 IdxRange[2][0] = 1;
1623 IdxRange[2][1] = nummodes[2] - 1;
1624
1625 if (edgeOrient == eBackwards)
1626 {
1627 reverseOrdering = true;
1628 }
1629 }
1630 else
1631 {
1632 IdxRange[2][0] = 2;
1633 IdxRange[2][1] = nummodes[2];
1634
1635 if (edgeOrient == eBackwards)
1636 {
1637 signChange = true;
1638 }
1639 }
1640 }
1641 break;
1642 }
1643
1644 switch (eid)
1645 {
1646 case 0:
1647 case 4:
1648 case 5:
1649 case 8:
1650 {
1651 IdxRange[1][0] = 0;
1652 IdxRange[1][1] = 1;
1653 }
1654 break;
1655 case 2:
1656 case 6:
1657 case 7:
1658 case 10:
1659 {
1660 if (bType[1] == LibUtilities::eGLL_Lagrange)
1661 {
1662 IdxRange[1][0] = nummodes[1] - 1;
1663 IdxRange[1][1] = nummodes[1];
1664 }
1665 else
1666 {
1667 IdxRange[1][0] = 1;
1668 IdxRange[1][1] = 2;
1669 }
1670 }
1671 break;
1672 case 1:
1673 case 9:
1674 {
1675 if (bType[1] == LibUtilities::eGLL_Lagrange)
1676 {
1677 IdxRange[1][0] = 1;
1678 IdxRange[1][1] = nummodes[1] - 1;
1679
1680 if (edgeOrient == eBackwards)
1681 {
1682 reverseOrdering = true;
1683 }
1684 }
1685 else
1686 {
1687 IdxRange[1][0] = 2;
1688 IdxRange[1][1] = nummodes[1];
1689
1690 if (edgeOrient == eBackwards)
1691 {
1692 signChange = true;
1693 }
1694 }
1695 }
1696 break;
1697 case 3:
1698 case 11:
1699 {
1700 if (bType[1] == LibUtilities::eGLL_Lagrange)
1701 {
1702 IdxRange[1][0] = 1;
1703 IdxRange[1][1] = nummodes[1] - 1;
1704
1705 if (edgeOrient == eForwards)
1706 {
1707 reverseOrdering = true;
1708 }
1709 }
1710 else
1711 {
1712 IdxRange[1][0] = 2;
1713 IdxRange[1][1] = nummodes[1];
1714
1715 if (edgeOrient == eForwards)
1716 {
1717 signChange = true;
1718 }
1719 }
1720 }
1721 break;
1722 }
1723
1724 switch (eid)
1725 {
1726 case 3:
1727 case 4:
1728 case 7:
1729 case 11:
1730 {
1731 IdxRange[0][0] = 0;
1732 IdxRange[0][1] = 1;
1733 }
1734 break;
1735 case 1:
1736 case 5:
1737 case 6:
1738 case 9:
1739 {
1740 if (bType[0] == LibUtilities::eGLL_Lagrange)
1741 {
1742 IdxRange[0][0] = nummodes[0] - 1;
1743 IdxRange[0][1] = nummodes[0];
1744 }
1745 else
1746 {
1747 IdxRange[0][0] = 1;
1748 IdxRange[0][1] = 2;
1749 }
1750 }
1751 break;
1752 case 0:
1753 case 8:
1754 {
1755 if (bType[0] == LibUtilities::eGLL_Lagrange)
1756 {
1757 IdxRange[0][0] = 1;
1758 IdxRange[0][1] = nummodes[0] - 1;
1759
1760 if (edgeOrient == eBackwards)
1761 {
1762 reverseOrdering = true;
1763 }
1764 }
1765 else
1766 {
1767 IdxRange[0][0] = 2;
1768 IdxRange[0][1] = nummodes[0];
1769
1770 if (edgeOrient == eBackwards)
1771 {
1772 signChange = true;
1773 }
1774 }
1775 }
1776 break;
1777 case 2:
1778 case 10:
1779 {
1780 if (bType[0] == LibUtilities::eGLL_Lagrange)
1781 {
1782 IdxRange[0][0] = 1;
1783 IdxRange[0][1] = nummodes[0] - 1;
1784
1785 if (edgeOrient == eForwards)
1786 {
1787 reverseOrdering = true;
1788 }
1789 }
1790 else
1791 {
1792 IdxRange[0][0] = 2;
1793 IdxRange[0][1] = nummodes[0];
1794
1795 if (edgeOrient == eForwards)
1796 {
1797 signChange = true;
1798 }
1799 }
1800 }
1801 break;
1802 }
1803
1804 int cnt = 0;
1805
1806 for (int r = IdxRange[2][0]; r < IdxRange[2][1]; r++)
1807 {
1808 for (int q = IdxRange[1][0]; q < IdxRange[1][1]; q++)
1809 {
1810 for (int p = IdxRange[0][0]; p < IdxRange[0][1]; p++)
1811 {
1812 maparray[cnt++] =
1813 r * nummodes[0] * nummodes[1] + q * nummodes[0] + p;
1814 }
1815 }
1816 }
1817
1818 if (reverseOrdering)
1819 {
1820 reverse(maparray.data(), maparray.data() + nEdgeIntCoeffs);
1821 }
1822
1823 if (signChange)
1824 {
1825 for (int p = 1; p < nEdgeIntCoeffs; p += 2)
1826 {
1827 signarray[p] = -1;
1828 }
1829 }
1830}
1831
1832/**
1833 * Generate mapping describing which elemental modes lie on the
1834 * interior of a given face. Accounts for face orientation.
1835 */
1837 const int fid, Array<OneD, unsigned int> &maparray,
1838 Array<OneD, int> &signarray, const Orientation faceOrient)
1839{
1842 "BasisType is not a boundary interior form");
1845 "BasisType is not a boundary interior form");
1848 "BasisType is not a boundary interior form");
1849
1850 ASSERTL1((fid >= 0) && (fid < 6), "local face id must be between 0 and 5");
1851
1852 int nFaceIntCoeffs = v_GetTraceIntNcoeffs(fid);
1853
1854 if (maparray.size() != nFaceIntCoeffs)
1855 {
1856 maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1857 }
1858
1859 if (signarray.size() != nFaceIntCoeffs)
1860 {
1861 signarray = Array<OneD, int>(nFaceIntCoeffs, 1);
1862 }
1863 else
1864 {
1865 fill(signarray.data(), signarray.data() + nFaceIntCoeffs, 1);
1866 }
1867
1868 int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
1869 m_base[2]->GetNumModes()};
1870
1871 const LibUtilities::BasisType bType[3] = {GetBasisType(0), GetBasisType(1),
1872 GetBasisType(2)};
1873
1874 int nummodesA = 0;
1875 int nummodesB = 0;
1876
1877 // Determine the number of modes in face directions A & B based
1878 // on the face index given.
1879 switch (fid)
1880 {
1881 case 0:
1882 case 5:
1883 {
1884 nummodesA = nummodes[0];
1885 nummodesB = nummodes[1];
1886 }
1887 break;
1888 case 1:
1889 case 3:
1890 {
1891 nummodesA = nummodes[0];
1892 nummodesB = nummodes[2];
1893 }
1894 break;
1895 case 2:
1896 case 4:
1897 {
1898 nummodesA = nummodes[1];
1899 nummodesB = nummodes[2];
1900 }
1901 }
1902
1903 Array<OneD, int> arrayindx(nFaceIntCoeffs);
1904
1905 // Create a mapping array to account for transposition of the
1906 // coordinates due to face orientation.
1907 for (int i = 0; i < (nummodesB - 2); i++)
1908 {
1909 for (int j = 0; j < (nummodesA - 2); j++)
1910 {
1911 if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1912 {
1913 arrayindx[i * (nummodesA - 2) + j] = i * (nummodesA - 2) + j;
1914 }
1915 else
1916 {
1917 arrayindx[i * (nummodesA - 2) + j] = j * (nummodesB - 2) + i;
1918 }
1919 }
1920 }
1921
1922 int IdxRange[3][2];
1923 int Incr[3];
1924
1925 Array<OneD, int> sign0(nummodes[0], 1);
1926 Array<OneD, int> sign1(nummodes[1], 1);
1927 Array<OneD, int> sign2(nummodes[2], 1);
1928
1929 // Set the upper and lower bounds, and increment for the faces
1930 // involving the first coordinate direction.
1931 switch (fid)
1932 {
1933 case 0: // bottom face
1934 {
1935 IdxRange[2][0] = 0;
1936 IdxRange[2][1] = 1;
1937 Incr[2] = 1;
1938 }
1939 break;
1940 case 5: // top face
1941 {
1942 if (bType[2] == LibUtilities::eGLL_Lagrange)
1943 {
1944 IdxRange[2][0] = nummodes[2] - 1;
1945 IdxRange[2][1] = nummodes[2];
1946 Incr[2] = 1;
1947 }
1948 else
1949 {
1950 IdxRange[2][0] = 1;
1951 IdxRange[2][1] = 2;
1952 Incr[2] = 1;
1953 }
1954 }
1955 break;
1956 default: // all other faces
1957 {
1958 if (bType[2] == LibUtilities::eGLL_Lagrange)
1959 {
1960 if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 2)
1961 {
1962 IdxRange[2][0] = nummodes[2] - 2;
1963 IdxRange[2][1] = 0;
1964 Incr[2] = -1;
1965 }
1966 else
1967 {
1968 IdxRange[2][0] = 1;
1969 IdxRange[2][1] = nummodes[2] - 1;
1970 Incr[2] = 1;
1971 }
1972 }
1973 else
1974 {
1975 IdxRange[2][0] = 2;
1976 IdxRange[2][1] = nummodes[2];
1977 Incr[2] = 1;
1978
1979 if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 2)
1980 {
1981 for (int i = 3; i < nummodes[2]; i += 2)
1982 {
1983 sign2[i] = -1;
1984 }
1985 }
1986 }
1987 }
1988 }
1989
1990 // Set the upper and lower bounds, and increment for the faces
1991 // involving the second coordinate direction.
1992 switch (fid)
1993 {
1994 case 1:
1995 {
1996 IdxRange[1][0] = 0;
1997 IdxRange[1][1] = 1;
1998 Incr[1] = 1;
1999 }
2000 break;
2001 case 3:
2002 {
2003 if (bType[1] == LibUtilities::eGLL_Lagrange)
2004 {
2005 IdxRange[1][0] = nummodes[1] - 1;
2006 IdxRange[1][1] = nummodes[1];
2007 Incr[1] = 1;
2008 }
2009 else
2010 {
2011 IdxRange[1][0] = 1;
2012 IdxRange[1][1] = 2;
2013 Incr[1] = 1;
2014 }
2015 }
2016 break;
2017 case 0:
2018 case 5:
2019 {
2020 if (bType[1] == LibUtilities::eGLL_Lagrange)
2021 {
2022 if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 2)
2023 {
2024 IdxRange[1][0] = nummodes[1] - 2;
2025 IdxRange[1][1] = 0;
2026 Incr[1] = -1;
2027 }
2028 else
2029 {
2030 IdxRange[1][0] = 1;
2031 IdxRange[1][1] = nummodes[1] - 1;
2032 Incr[1] = 1;
2033 }
2034 }
2035 else
2036 {
2037 IdxRange[1][0] = 2;
2038 IdxRange[1][1] = nummodes[1];
2039 Incr[1] = 1;
2040
2041 if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 2)
2042 {
2043 for (int i = 3; i < nummodes[1]; i += 2)
2044 {
2045 sign1[i] = -1;
2046 }
2047 }
2048 }
2049 }
2050 break;
2051 default: // case2: case4:
2052 {
2053 if (bType[1] == LibUtilities::eGLL_Lagrange)
2054 {
2055 if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1)
2056 {
2057 IdxRange[1][0] = nummodes[1] - 2;
2058 IdxRange[1][1] = 0;
2059 Incr[1] = -1;
2060 }
2061 else
2062 {
2063 IdxRange[1][0] = 1;
2064 IdxRange[1][1] = nummodes[1] - 1;
2065 Incr[1] = 1;
2066 }
2067 }
2068 else
2069 {
2070 IdxRange[1][0] = 2;
2071 IdxRange[1][1] = nummodes[1];
2072 Incr[1] = 1;
2073
2074 if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1)
2075 {
2076 for (int i = 3; i < nummodes[1]; i += 2)
2077 {
2078 sign1[i] = -1;
2079 }
2080 }
2081 }
2082 }
2083 }
2084
2085 switch (fid)
2086 {
2087 case 4:
2088 {
2089 IdxRange[0][0] = 0;
2090 IdxRange[0][1] = 1;
2091 Incr[0] = 1;
2092 }
2093 break;
2094 case 2:
2095 {
2096 if (bType[0] == LibUtilities::eGLL_Lagrange)
2097 {
2098 IdxRange[0][0] = nummodes[0] - 1;
2099 IdxRange[0][1] = nummodes[0];
2100 Incr[0] = 1;
2101 }
2102 else
2103 {
2104 IdxRange[0][0] = 1;
2105 IdxRange[0][1] = 2;
2106 Incr[0] = 1;
2107 }
2108 }
2109 break;
2110 default:
2111 {
2112 if (bType[0] == LibUtilities::eGLL_Lagrange)
2113 {
2114 if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1)
2115 {
2116 IdxRange[0][0] = nummodes[0] - 2;
2117 IdxRange[0][1] = 0;
2118 Incr[0] = -1;
2119 }
2120 else
2121 {
2122 IdxRange[0][0] = 1;
2123 IdxRange[0][1] = nummodes[0] - 1;
2124 Incr[0] = 1;
2125 }
2126 }
2127 else
2128 {
2129 IdxRange[0][0] = 2;
2130 IdxRange[0][1] = nummodes[0];
2131 Incr[0] = 1;
2132
2133 if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1)
2134 {
2135 for (int i = 3; i < nummodes[0]; i += 2)
2136 {
2137 sign0[i] = -1;
2138 }
2139 }
2140 }
2141 }
2142 }
2143
2144 int cnt = 0;
2145
2146 for (int r = IdxRange[2][0]; r != IdxRange[2][1]; r += Incr[2])
2147 {
2148 for (int q = IdxRange[1][0]; q != IdxRange[1][1]; q += Incr[1])
2149 {
2150 for (int p = IdxRange[0][0]; p != IdxRange[0][1]; p += Incr[0])
2151 {
2152 maparray[arrayindx[cnt]] =
2153 r * nummodes[0] * nummodes[1] + q * nummodes[0] + p;
2154 signarray[arrayindx[cnt++]] = sign0[p] * sign1[q] * sign2[r];
2155 }
2156 }
2157 }
2158}
2159
2160int StdHexExp::v_GetEdgeNcoeffs(const int i) const
2161{
2162 ASSERTL2((i >= 0) && (i <= 11), "edge id is out of range");
2163
2164 if ((i == 0) || (i == 2) || (i == 8) || (i == 10))
2165 {
2166 return GetBasisNumModes(0);
2167 }
2168 else if ((i == 1) || (i == 3) || (i == 9) || (i == 11))
2169 {
2170 return GetBasisNumModes(1);
2171 }
2172 else
2173 {
2174 return GetBasisNumModes(2);
2175 }
2176}
2177
2179{
2180
2181 MatrixType mtype = mkey.GetMatrixType();
2182
2183 DNekMatSharedPtr Mat;
2184
2185 switch (mtype)
2186 {
2188 {
2189 int nq0 = m_base[0]->GetNumPoints();
2190 int nq1 = m_base[1]->GetNumPoints();
2191 int nq2 = m_base[2]->GetNumPoints();
2192 int nq;
2193
2194 // take definition from key
2196 {
2197 nq = (int)mkey.GetConstFactor(eFactorConst);
2198 }
2199 else
2200 {
2201 nq = max(nq0, max(nq1, nq2));
2202 }
2203
2204 int neq =
2207 Array<OneD, NekDouble> coll(3);
2209 Array<OneD, NekDouble> tmp(nq0);
2210
2211 Mat =
2212 MemoryManager<DNekMat>::AllocateSharedPtr(neq, nq0 * nq1 * nq2);
2213 int cnt = 0;
2214
2215 for (int i = 0; i < nq; ++i)
2216 {
2217 for (int j = 0; j < nq; ++j)
2218 {
2219 for (int k = 0; k < nq; ++k, ++cnt)
2220 {
2221 coords[cnt] = Array<OneD, NekDouble>(3);
2222 coords[cnt][0] = -1.0 + 2 * k / (NekDouble)(nq - 1);
2223 coords[cnt][1] = -1.0 + 2 * j / (NekDouble)(nq - 1);
2224 coords[cnt][2] = -1.0 + 2 * i / (NekDouble)(nq - 1);
2225 }
2226 }
2227 }
2228
2229 for (int i = 0; i < neq; ++i)
2230 {
2231 LocCoordToLocCollapsed(coords[i], coll);
2232
2233 I[0] = m_base[0]->GetI(coll);
2234 I[1] = m_base[1]->GetI(coll + 1);
2235 I[2] = m_base[2]->GetI(coll + 2);
2236
2237 // interpolate first coordinate direction
2238 NekDouble fac;
2239 for (int k = 0; k < nq2; ++k)
2240 {
2241 for (int j = 0; j < nq1; ++j)
2242 {
2243
2244 fac = (I[1]->GetPtr())[j] * (I[2]->GetPtr())[k];
2245 Vmath::Smul(nq0, fac, I[0]->GetPtr(), 1, tmp, 1);
2246
2247 Vmath::Vcopy(nq0, &tmp[0], 1,
2248 Mat->GetRawPtr() + k * nq0 * nq1 * neq +
2249 j * nq0 * neq + i,
2250 neq);
2251 }
2252 }
2253 }
2254 }
2255 break;
2256 default:
2257 {
2259 }
2260 break;
2261 }
2262
2263 return Mat;
2264}
2265
2267{
2268 return v_GenMatrix(mkey);
2269}
2270
2272 Array<OneD, NekDouble> &outarray,
2273 const StdMatrixKey &mkey)
2274{
2275 StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
2276}
2277
2279 Array<OneD, NekDouble> &outarray,
2280 const StdMatrixKey &mkey)
2281{
2282 StdHexExp::v_LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
2283}
2284
2285void StdHexExp::v_LaplacianMatrixOp(const int k1, const int k2,
2286 const Array<OneD, const NekDouble> &inarray,
2287 Array<OneD, NekDouble> &outarray,
2288 const StdMatrixKey &mkey)
2289{
2290 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
2291}
2292
2294 const Array<OneD, const NekDouble> &inarray,
2295 Array<OneD, NekDouble> &outarray,
2296 const StdMatrixKey &mkey)
2297{
2298 StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
2299}
2300
2302 Array<OneD, NekDouble> &outarray,
2303 const StdMatrixKey &mkey)
2304{
2305 StdHexExp::v_HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
2306}
2307
2309 const Array<OneD, const NekDouble> &inarray,
2310 Array<OneD, NekDouble> &outarray)
2311{
2312 int nquad0 = m_base[0]->GetNumPoints();
2313 int nquad1 = m_base[1]->GetNumPoints();
2314 int nquad2 = m_base[2]->GetNumPoints();
2315 int nq01 = nquad0 * nquad1;
2316 int nq12 = nquad1 * nquad2;
2317
2318 const Array<OneD, const NekDouble> &w0 = m_base[0]->GetW();
2319 const Array<OneD, const NekDouble> &w1 = m_base[1]->GetW();
2320 const Array<OneD, const NekDouble> &w2 = m_base[2]->GetW();
2321
2322 for (int i = 0; i < nq12; ++i)
2323 {
2324 Vmath::Vmul(nquad0, inarray.data() + i * nquad0, 1, w0.data(), 1,
2325 outarray.data() + i * nquad0, 1);
2326 }
2327
2328 for (int i = 0; i < nq12; ++i)
2329 {
2330 Vmath::Smul(nquad0, w1[i % nquad1], outarray.data() + i * nquad0, 1,
2331 outarray.data() + i * nquad0, 1);
2332 }
2333
2334 for (int i = 0; i < nquad2; ++i)
2335 {
2336 Vmath::Smul(nq01, w2[i], outarray.data() + i * nq01, 1,
2337 outarray.data() + i * nq01, 1);
2338 }
2339}
2340
2342 const StdMatrixKey &mkey)
2343{
2344 // Generate an orthonogal expansion
2345 int qa = m_base[0]->GetNumPoints();
2346 int qb = m_base[1]->GetNumPoints();
2347 int qc = m_base[2]->GetNumPoints();
2348 int nmodes_a = m_base[0]->GetNumModes();
2349 int nmodes_b = m_base[1]->GetNumModes();
2350 int nmodes_c = m_base[2]->GetNumModes();
2351 // Declare orthogonal basis.
2355
2359 StdHexExp OrthoExp(Ba, Bb, Bc);
2360
2361 Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2362 int cnt = 0;
2363
2364 // project onto modal space.
2365 OrthoExp.FwdTrans(array, orthocoeffs);
2366
2368 {
2369 // Rodrigo's power kernel
2371 NekDouble SvvDiffCoeff =
2374
2375 for (int i = 0; i < nmodes_a; ++i)
2376 {
2377 for (int j = 0; j < nmodes_b; ++j)
2378 {
2379 NekDouble fac1 = std::max(
2380 pow((1.0 * i) / (nmodes_a - 1), cutoff * nmodes_a),
2381 pow((1.0 * j) / (nmodes_b - 1), cutoff * nmodes_b));
2382
2383 for (int k = 0; k < nmodes_c; ++k)
2384 {
2385 NekDouble fac =
2386 std::max(fac1, pow((1.0 * k) / (nmodes_c - 1),
2387 cutoff * nmodes_c));
2388
2389 orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2390 cnt++;
2391 }
2392 }
2393 }
2394 }
2395 else if (mkey.ConstFactorExists(
2396 eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
2397 {
2400
2401 int max_abc = max(nmodes_a - kSVVDGFiltermodesmin,
2402 nmodes_b - kSVVDGFiltermodesmin);
2403 max_abc = max(max_abc, nmodes_c - kSVVDGFiltermodesmin);
2404 // clamp max_abc
2405 max_abc = max(max_abc, 0);
2406 max_abc = min(max_abc, kSVVDGFiltermodesmax - kSVVDGFiltermodesmin);
2407
2408 for (int i = 0; i < nmodes_a; ++i)
2409 {
2410 for (int j = 0; j < nmodes_b; ++j)
2411 {
2412 int maxij = max(i, j);
2413
2414 for (int k = 0; k < nmodes_c; ++k)
2415 {
2416 int maxijk = max(maxij, k);
2417 maxijk = min(maxijk, kSVVDGFiltermodesmax - 1);
2418
2419 orthocoeffs[cnt] *=
2420 SvvDiffCoeff * kSVVDGFilter[max_abc][maxijk];
2421 cnt++;
2422 }
2423 }
2424 }
2425 }
2426 else
2427 {
2428
2429 int cutoff = (int)(mkey.GetConstFactor(eFactorSVVCutoffRatio) *
2430 min(nmodes_a, nmodes_b));
2431 NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
2432 // Filter just trilinear space
2433 int nmodes = max(nmodes_a, nmodes_b);
2434 nmodes = max(nmodes, nmodes_c);
2435
2436 Array<OneD, NekDouble> fac(nmodes, 1.0);
2437 for (int j = cutoff; j < nmodes; ++j)
2438 {
2439 fac[j] = fabs((j - nmodes) / ((NekDouble)(j - cutoff + 1.0)));
2440 fac[j] *= fac[j]; // added this line to conform with equation
2441 }
2442
2443 for (int i = 0; i < nmodes_a; ++i)
2444 {
2445 for (int j = 0; j < nmodes_b; ++j)
2446 {
2447 for (int k = 0; k < nmodes_c; ++k)
2448 {
2449 if ((i >= cutoff) || (j >= cutoff) || (k >= cutoff))
2450 {
2451 orthocoeffs[i * nmodes_a * nmodes_b + j * nmodes_c +
2452 k] *=
2453 (SvvDiffCoeff * exp(-(fac[i] + fac[j] + fac[k])));
2454 }
2455 else
2456 {
2457 orthocoeffs[i * nmodes_a * nmodes_b + j * nmodes_c +
2458 k] *= 0.0;
2459 }
2460 }
2461 }
2462 }
2463 }
2464
2465 // backward transform to physical space
2466 OrthoExp.BwdTrans(orthocoeffs, array);
2467}
2468
2470 const NekDouble alpha,
2471 const NekDouble exponent,
2472 const NekDouble cutoff)
2473{
2474 // Generate an orthogonal expansion
2475 int qa = m_base[0]->GetNumPoints();
2476 int qb = m_base[1]->GetNumPoints();
2477 int qc = m_base[2]->GetNumPoints();
2478 int nmodesA = m_base[0]->GetNumModes();
2479 int nmodesB = m_base[1]->GetNumModes();
2480 int nmodesC = m_base[2]->GetNumModes();
2481 int P = nmodesA - 1;
2482 int Q = nmodesB - 1;
2483 int R = nmodesC - 1;
2484
2485 // Declare orthogonal basis.
2489
2493 StdHexExp OrthoExp(Ba, Bb, Bc);
2494
2495 // Cutoff
2496 int Pcut = cutoff * P;
2497 int Qcut = cutoff * Q;
2498 int Rcut = cutoff * R;
2499
2500 // Project onto orthogonal space.
2501 Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2502 OrthoExp.FwdTrans(array, orthocoeffs);
2503
2504 //
2505 NekDouble fac, fac1, fac2, fac3;
2506 int index = 0;
2507 for (int i = 0; i < nmodesA; ++i)
2508 {
2509 for (int j = 0; j < nmodesB; ++j)
2510 {
2511 for (int k = 0; k < nmodesC; ++k, ++index)
2512 {
2513 // to filter out only the "high-modes"
2514 if (i > Pcut || j > Qcut || k > Rcut)
2515 {
2516 fac1 = (NekDouble)(i - Pcut) / ((NekDouble)(P - Pcut));
2517 fac2 = (NekDouble)(j - Qcut) / ((NekDouble)(Q - Qcut));
2518 fac3 = (NekDouble)(k - Rcut) / ((NekDouble)(R - Rcut));
2519 fac = max(max(fac1, fac2), fac3);
2520 fac = pow(fac, exponent);
2521 orthocoeffs[index] *= exp(-alpha * fac);
2522 }
2523 }
2524 }
2525 }
2526
2527 // backward transform to physical space
2528 OrthoExp.BwdTrans(orthocoeffs, array);
2529}
2530
2532 Array<OneD, int> &conn, [[maybe_unused]] bool standard)
2533{
2534 int np0 = m_base[0]->GetNumPoints();
2535 int np1 = m_base[1]->GetNumPoints();
2536 int np2 = m_base[2]->GetNumPoints();
2537 int np = max(np0, max(np1, np2));
2538
2539 conn = Array<OneD, int>(6 * (np - 1) * (np - 1) * (np - 1));
2540
2541 int row = 0;
2542 int rowp1 = 0;
2543 int cnt = 0;
2544 int plane = 0;
2545 for (int i = 0; i < np - 1; ++i)
2546 {
2547 for (int j = 0; j < np - 1; ++j)
2548 {
2549 rowp1 += np;
2550 for (int k = 0; k < np - 1; ++k)
2551 {
2552 conn[cnt++] = plane + row + k;
2553 conn[cnt++] = plane + row + k + 1;
2554 conn[cnt++] = plane + rowp1 + k;
2555
2556 conn[cnt++] = plane + rowp1 + k + 1;
2557 conn[cnt++] = plane + rowp1 + k;
2558 conn[cnt++] = plane + row + k + 1;
2559 }
2560 row += np;
2561 }
2562 plane += np * np;
2563 }
2564}
2565} // 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
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)
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
NekDouble BaryTensorDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs)
void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
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
void WeakDerivMatrixOp_MatFree(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
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.
void BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
this function calculates the inner product of a given function f with the different modes of the expa...
Definition: StdExpansion.h:537
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:370
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space.
Definition: StdExpansion.h:433
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:205
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space.
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:218
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:732
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
void MassMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Class representing a hexehedral element in reference space.
Definition: StdHexExp.h:44
int v_NumDGBndryCoeffs() const override
Definition: StdHexExp.cpp:630
void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
Definition: StdHexExp.cpp:2271
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdHexExp.cpp:2308
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
Differentiation Methods.
Definition: StdHexExp.cpp:72
void v_ExponentialFilter(Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff) override
Definition: StdHexExp.cpp:2469
int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false) override
Definition: StdHexExp.cpp:835
void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta) override
Definition: StdHexExp.cpp:505
int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset) override
Definition: StdHexExp.cpp:721
LibUtilities::PointsKey v_GetTracePointsKey(const int i, const int j) const override
Definition: StdHexExp.cpp:701
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) final
Definition: StdHexExp.cpp:569
void v_GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
Definition: StdHexExp.cpp:1836
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdHexExp.cpp:459
DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey) override
Definition: StdHexExp.cpp:2266
int v_GetEdgeNcoeffs(const int i) const override
Definition: StdHexExp.cpp:2160
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multbyweights=true) override
Definition: StdHexExp.cpp:337
StdHexExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc)
Definition: StdHexExp.cpp:45
int v_GetTraceNumPoints(const int i) const override
Definition: StdHexExp.cpp:683
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
Definition: StdHexExp.cpp:1113
void v_GetTraceCoeffMap(const unsigned int fid, Array< OneD, unsigned int > &maparray) override
Definition: StdHexExp.cpp:1124
LibUtilities::ShapeType v_DetShapeType() const override
Definition: StdHexExp.cpp:605
void v_GetTraceNumModes(const int fid, int &numModes0, int &numModes1, Orientation faceOrient=eDir1FwdDir1_Dir2FwdDir2) override
Definition: StdHexExp.cpp:787
void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi) override
Definition: StdHexExp.cpp:513
void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
Definition: StdHexExp.cpp:2278
int v_NumBndryCoeffs() const override
Definition: StdHexExp.cpp:610
int v_GetNedges() const override
Definition: StdHexExp.cpp:595
void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdHexExp.cpp:262
int v_GetNtraces() const override
Definition: StdHexExp.cpp:600
void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
Definition: StdHexExp.cpp:2301
int v_GetNverts() const override
Definition: StdHexExp.cpp:590
bool v_IsBoundaryInteriorExpansion() const override
Definition: StdHexExp.cpp:55
void v_GetInteriorMap(Array< OneD, unsigned int > &outarray) override
Definition: StdHexExp.cpp:963
void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray) override
Definition: StdHexExp.cpp:1024
void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray) override
Definition: StdHexExp.cpp:524
void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true) override
Definition: StdHexExp.cpp:2531
void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
Definition: StdHexExp.cpp:2293
int v_GetTraceIntNcoeffs(const int i) const override
Definition: StdHexExp.cpp:666
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: StdHexExp.cpp:118
void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdHexExp.cpp:320
void v_GetElmtTraceToTraceMap(const unsigned int fid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation faceOrient, int P, int Q) override
Definition: StdHexExp.cpp:1248
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: StdHexExp.cpp:205
void v_GetEdgeInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2) override
Definition: StdHexExp.cpp:1543
void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdHexExp.cpp:153
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey) override
Definition: StdHexExp.cpp:2341
const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int k, bool useGLL=false) const override
Definition: StdHexExp.cpp:731
void v_GetCoords(Array< OneD, NekDouble > &coords_x, Array< OneD, NekDouble > &coords_y, Array< OneD, NekDouble > &coords_z) override
Definition: StdHexExp.cpp:759
DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey) override
Definition: StdHexExp.cpp:2178
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdHexExp.cpp:452
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: StdHexExp.cpp:371
int v_GetTraceNcoeffs(const int i) const override
Definition: StdHexExp.cpp:649
void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdHexExp.cpp:180
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
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 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
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:153
@ 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
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 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.