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