Nektar++
StdQuadExp.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: StdQuadExp.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: Quadrilateral routines built upon StdExpansion2D
32//
33///////////////////////////////////////////////////////////////////////////////
34
38
39using namespace std;
40
41namespace Nektar::StdRegions
42{
43
44/** \brief Constructor using BasisKey class for quadrature
45 * points and order definition
46 */
48 const LibUtilities::BasisKey &Bb)
49 : StdExpansion(Ba.GetNumModes() * Bb.GetNumModes(), 2, Ba, Bb),
50 StdExpansion2D(Ba.GetNumModes() * Bb.GetNumModes(), Ba, Bb)
51{
52}
53
54/////////////////////////
55// Integration Methods //
56/////////////////////////
57
59{
62
63 return StdExpansion2D::Integral(inarray, w0, w1);
64}
65
66/////////////////////////////
67// Differentiation Methods //
68/////////////////////////////
69
70/** \brief Calculate the derivative of the physical points
71 *
72 * For quadrilateral region can use the Tensor_Deriv function
73 * defined under StdExpansion.
74 */
75
79 [[maybe_unused]] Array<OneD, NekDouble> &out_d2)
80{
81 PhysTensorDeriv(inarray, out_d0, out_d1);
82}
83
84void StdQuadExp::v_PhysDeriv(const int dir,
85 const Array<OneD, const NekDouble> &inarray,
86 Array<OneD, NekDouble> &outarray)
87{
88 switch (dir)
89 {
90 case 0:
91 {
92 PhysTensorDeriv(inarray, outarray, NullNekDouble1DArray);
93 }
94 break;
95 case 1:
96 {
97 PhysTensorDeriv(inarray, NullNekDouble1DArray, outarray);
98 }
99 break;
100 default:
101 {
102 ASSERTL1(false, "input dir is out of range");
103 }
104 break;
105 }
106}
107
111 [[maybe_unused]] Array<OneD, NekDouble> &out_d2)
112{
113 StdQuadExp::v_PhysDeriv(inarray, out_d0, out_d1);
114}
115
116void StdQuadExp::v_StdPhysDeriv(const int dir,
117 const Array<OneD, const NekDouble> &inarray,
118 Array<OneD, NekDouble> &outarray)
119{
120 StdQuadExp::v_PhysDeriv(dir, inarray, outarray);
121}
122
123////////////////
124// Transforms //
125////////////////
126
128 Array<OneD, NekDouble> &outarray)
129{
130 if (m_base[0]->Collocation() && m_base[1]->Collocation())
131 {
133 inarray, 1, outarray, 1);
134 }
135 else
136 {
137 StdQuadExp::v_BwdTrans_SumFac(inarray, outarray);
138 }
139}
140
142 Array<OneD, NekDouble> &outarray)
143{
145 m_base[1]->GetNumModes());
146
147 BwdTrans_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(), inarray,
148 outarray, wsp, true, true);
149}
150
151// The arguments doCheckCollDir0 and doCheckCollDir1 allow you to specify
152// whether to check if the basis has collocation properties (i.e. for the
153// classical spectral element basis, In this case the 1D 'B' matrix is equal to
154// the identity matrix which can be exploited to speed up the calculations).
155// However, as this routine also allows to pass the matrix 'DB' (derivative of
156// the basis), the collocation property cannot always be used. Therefor follow
157// this rule: if base0 == m_base[0]->GetBdata() --> set doCheckCollDir0 == true;
158// base1 == m_base[1]->GetBdata() --> set doCheckCollDir1 == true;
159// base0 == m_base[0]->GetDbdata() --> set doCheckCollDir0 == false;
160// base1 == m_base[1]->GetDbdata() --> set doCheckCollDir1 == false;
162 const Array<OneD, const NekDouble> &base0,
163 const Array<OneD, const NekDouble> &base1,
164 const Array<OneD, const NekDouble> &inarray,
166 bool doCheckCollDir0, bool doCheckCollDir1)
167{
168 int nquad0 = m_base[0]->GetNumPoints();
169 int nquad1 = m_base[1]->GetNumPoints();
170 int nmodes0 = m_base[0]->GetNumModes();
171 int nmodes1 = m_base[1]->GetNumModes();
172
173 bool colldir0 = doCheckCollDir0 ? (m_base[0]->Collocation()) : false;
174 bool colldir1 = doCheckCollDir1 ? (m_base[1]->Collocation()) : false;
175
176 if (colldir0 && colldir1)
177 {
178 Vmath::Vcopy(m_ncoeffs, inarray.data(), 1, outarray.data(), 1);
179 }
180 else if (colldir0)
181 {
182 Blas::Dgemm('N', 'T', nquad0, nquad1, nmodes1, 1.0, &inarray[0], nquad0,
183 base1.data(), nquad1, 0.0, &outarray[0], nquad0);
184 }
185 else if (colldir1)
186 {
187 Blas::Dgemm('N', 'N', nquad0, nmodes1, nmodes0, 1.0, base0.data(),
188 nquad0, &inarray[0], nmodes0, 0.0, &outarray[0], nquad0);
189 }
190 else
191 {
192 ASSERTL1(wsp.size() >= nquad0 * nmodes1,
193 "Workspace size is not sufficient");
194
195 // Those two calls correpsond to the operation
196 // out = B0*in*Transpose(B1);
197 Blas::Dgemm('N', 'N', nquad0, nmodes1, nmodes0, 1.0, base0.data(),
198 nquad0, &inarray[0], nmodes0, 0.0, &wsp[0], nquad0);
199 Blas::Dgemm('N', 'T', nquad0, nquad1, nmodes1, 1.0, &wsp[0], nquad0,
200 base1.data(), nquad1, 0.0, &outarray[0], nquad0);
201 }
202}
203
205 Array<OneD, NekDouble> &outarray)
206{
207 if ((m_base[0]->Collocation()) && (m_base[1]->Collocation()))
208 {
209 Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
210 }
211 else
212 {
213 StdQuadExp::v_IProductWRTBase(inarray, outarray);
214
215 // get Mass matrix inverse
216 StdMatrixKey masskey(eInvMass, DetShapeType(), *this);
217 DNekMatSharedPtr matsys = GetStdMatrix(masskey);
218
219 // copy inarray in case inarray == outarray
220 NekVector<NekDouble> in(m_ncoeffs, outarray, eCopy);
222
223 out = (*matsys) * in;
224 }
225}
226
228 const Array<OneD, const NekDouble> &inarray,
229 Array<OneD, NekDouble> &outarray)
230{
231 if ((m_base[0]->Collocation()) && (m_base[1]->Collocation()))
232 {
233 Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
234 }
235 else
236 {
237 int i, j;
238 int npoints[2] = {m_base[0]->GetNumPoints(), m_base[1]->GetNumPoints()};
239 int nmodes[2] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes()};
240
241 fill(outarray.data(), outarray.data() + m_ncoeffs, 0.0);
242
243 Array<OneD, NekDouble> physEdge[4];
244 Array<OneD, NekDouble> coeffEdge[4];
245 for (i = 0; i < 4; i++)
246 {
247 physEdge[i] = Array<OneD, NekDouble>(npoints[i % 2]);
248 coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i % 2]);
249 }
250
251 for (i = 0; i < npoints[0]; i++)
252 {
253 physEdge[0][i] = inarray[i];
254 physEdge[2][i] = inarray[npoints[0] * npoints[1] - 1 - i];
255 }
256
257 for (i = 0; i < npoints[1]; i++)
258 {
259 physEdge[1][i] = inarray[npoints[0] - 1 + i * npoints[0]];
260 physEdge[3][i] =
261 inarray[(npoints[1] - 1) * npoints[0] - i * npoints[0]];
262 }
263
264 StdSegExpSharedPtr segexp[2] = {
266 m_base[0]->GetBasisKey()),
268 m_base[1]->GetBasisKey())};
269
271 Array<OneD, int> signArray;
273
274 for (i = 0; i < 4; i++)
275 {
276 segexp[i % 2]->FwdTransBndConstrained(physEdge[i], coeffEdge[i]);
277
278 GetTraceToElementMap(i, mapArray, signArray);
279 for (j = 0; j < nmodes[i % 2]; j++)
280 {
281 sign = (NekDouble)signArray[j];
282 outarray[mapArray[j]] = sign * coeffEdge[i][j];
283 }
284 }
285
288
289 StdMatrixKey masskey(eMass, DetShapeType(), *this);
290 MassMatrixOp(outarray, tmp0, masskey);
291 IProductWRTBase(inarray, tmp1);
292
293 Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
294
295 // get Mass matrix inverse (only of interior DOF)
296 // use block (1,1) of the static condensed system
297 // note: this block alreay contains the inverse matrix
298 DNekMatSharedPtr matsys =
299 (m_stdStaticCondMatrixManager[masskey])->GetBlock(1, 1);
300
301 int nBoundaryDofs = NumBndryCoeffs();
302 int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
303
304 Array<OneD, NekDouble> rhs(nInteriorDofs);
305 Array<OneD, NekDouble> result(nInteriorDofs);
306
307 GetInteriorMap(mapArray);
308
309 for (i = 0; i < nInteriorDofs; i++)
310 {
311 rhs[i] = tmp1[mapArray[i]];
312 }
313
314 Blas::Dgemv('N', nInteriorDofs, nInteriorDofs, 1.0,
315 &(matsys->GetPtr())[0], nInteriorDofs, rhs.data(), 1, 0.0,
316 result.data(), 1);
317
318 for (i = 0; i < nInteriorDofs; i++)
319 {
320 outarray[mapArray[i]] = result[i];
321 }
322 }
323}
324
325/////////////////////////////
326// Inner Product Functions //
327/////////////////////////////
328
329/** \brief Calculate the inner product of inarray with respect to
330 * the basis B=base0*base1 and put into outarray
331 *
332 * \f$
333 * \begin{array}{rcl}
334 * I_{pq} = (\phi_p \phi_q, u) & = & \sum_{i=0}^{nq_0}
335 * \sum_{j=0}^{nq_1}
336 * \phi_p(\xi_{0,i}) \phi_q(\xi_{1,j}) w^0_i w^1_j u(\xi_{0,i}
337 * \xi_{1,j}) \\
338 * & = & \sum_{i=0}^{nq_0} \phi_p(\xi_{0,i})
339 * \sum_{j=0}^{nq_1} \phi_q(\xi_{1,j}) \tilde{u}_{i,j}
340 * \end{array}
341 * \f$
342 *
343 * where
344 *
345 * \f$ \tilde{u}_{i,j} = w^0_i w^1_j u(\xi_{0,i},\xi_{1,j}) \f$
346 *
347 * which can be implemented as
348 *
349 * \f$ f_{qi} = \sum_{j=0}^{nq_1} \phi_q(\xi_{1,j})
350 * \tilde{u}_{i,j} = {\bf B_1 U} \f$
351 * \f$ I_{pq} = \sum_{i=0}^{nq_0} \phi_p(\xi_{0,i}) f_{qi} =
352 * {\bf B_0 F} \f$
353 */
355 Array<OneD, NekDouble> &outarray)
356{
357 if (m_base[0]->Collocation() && m_base[1]->Collocation())
358 {
359 MultiplyByQuadratureMetric(inarray, outarray);
360 }
361 else
362 {
363 StdQuadExp::v_IProductWRTBase_SumFac(inarray, outarray);
364 }
365}
366
368 const Array<OneD, const NekDouble> &inarray,
369 Array<OneD, NekDouble> &outarray, bool multiplybyweights)
370{
371 int nquad0 = m_base[0]->GetNumPoints();
372 int nquad1 = m_base[1]->GetNumPoints();
373 int order0 = m_base[0]->GetNumModes();
374
375 if (multiplybyweights)
376 {
377 Array<OneD, NekDouble> tmp(nquad0 * nquad1 + nquad1 * order0);
378 Array<OneD, NekDouble> wsp(tmp + nquad0 * nquad1);
379
380 // multiply by integration constants
381 MultiplyByQuadratureMetric(inarray, tmp);
383 m_base[1]->GetBdata(), tmp, outarray, wsp,
384 true, true);
385 }
386 else
387 {
388 Array<OneD, NekDouble> wsp(nquad1 * order0);
390 m_base[1]->GetBdata(), inarray, outarray,
391 wsp, true, true);
392 }
393}
394
396 const int dir, const Array<OneD, const NekDouble> &inarray,
397 Array<OneD, NekDouble> &outarray)
398{
399 v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
400}
401
403 const int dir, const Array<OneD, const NekDouble> &inarray,
404 Array<OneD, NekDouble> &outarray)
405{
406 ASSERTL0((dir == 0) || (dir == 1), "input dir is out of range");
407
408 int nquad0 = m_base[0]->GetNumPoints();
409 int nquad1 = m_base[1]->GetNumPoints();
410 int nqtot = nquad0 * nquad1;
411 int order0 = m_base[0]->GetNumModes();
412
413 Array<OneD, NekDouble> tmp(nqtot + nquad1 * order0);
414 Array<OneD, NekDouble> wsp(tmp + nqtot);
415
416 // multiply by integration constants
417 MultiplyByQuadratureMetric(inarray, tmp);
418
419 if (dir) // dir == 1
420 {
422 m_base[1]->GetDbdata(), tmp, outarray, wsp,
423 true, false);
424 }
425 else // dir == 0
426 {
427 IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
428 m_base[1]->GetBdata(), tmp, outarray, wsp,
429 false, true);
430 }
431}
432
433// the arguments doCheckCollDir0 and doCheckCollDir1 allow you to specify
434// whether to check if the basis has collocation properties (i.e. for the
435// classical spectral element basis, In this case the 1D 'B' matrix is equal to
436// the identity matrix which can be exploited to speed up the calculations).
437// However, as this routine also allows to pass the matrix 'DB' (derivative of
438// the basis), the collocation property cannot always be used. Therefor follow
439// this rule: if base0 == m_base[0]->GetBdata() --> set doCheckCollDir0 == true;
440// base1 == m_base[1]->GetBdata() --> set doCheckCollDir1 == true;
441// base0 == m_base[0]->GetDbdata() --> set doCheckCollDir0 == false;
442// base1 == m_base[1]->GetDbdata() --> set doCheckCollDir1 == false;
444 const Array<OneD, const NekDouble> &base0,
445 const Array<OneD, const NekDouble> &base1,
446 const Array<OneD, const NekDouble> &inarray,
448 bool doCheckCollDir0, bool doCheckCollDir1)
449{
450 int nquad0 = m_base[0]->GetNumPoints();
451 int nquad1 = m_base[1]->GetNumPoints();
452 int nmodes0 = m_base[0]->GetNumModes();
453 int nmodes1 = m_base[1]->GetNumModes();
454
455 bool colldir0 = doCheckCollDir0 ? (m_base[0]->Collocation()) : false;
456 bool colldir1 = doCheckCollDir1 ? (m_base[1]->Collocation()) : false;
457
458 if (colldir0 && colldir1)
459 {
460 Vmath::Vcopy(m_ncoeffs, inarray.data(), 1, outarray.data(), 1);
461 }
462 else if (colldir0)
463 {
464 Blas::Dgemm('N', 'N', nmodes0, nmodes1, nquad1, 1.0, inarray.data(),
465 nmodes0, base1.data(), nquad1, 0.0, outarray.data(),
466 nmodes0);
467 }
468 else if (colldir1)
469 {
470 Blas::Dgemm('T', 'N', nmodes0, nquad1, nquad0, 1.0, base0.data(),
471 nquad0, inarray.data(), nquad0, 0.0, outarray.data(),
472 nmodes0);
473 }
474 else
475 {
476 ASSERTL1(wsp.size() >= nquad1 * nmodes0,
477 "Workspace size is not sufficient");
478
479#if 1
480 Blas::Dgemm('T', 'N', nmodes0, nquad1, nquad0, 1.0, base0.data(),
481 nquad0, inarray.data(), nquad0, 0.0, wsp.data(), nmodes0);
482
483#else
484 for (int i = 0; i < nmodes0; ++i)
485 {
486 for (int j = 0; j < nquad1; ++j)
487 {
488 wsp[j * nmodes0 + i] =
489 Blas::Ddot(nquad0, base0.data() + i * nquad0, 1,
490 inarray.data() + j * nquad0, 1);
491 }
492 }
493#endif
494 Blas::Dgemm('N', 'N', nmodes0, nmodes1, nquad1, 1.0, wsp.data(),
495 nmodes0, base1.data(), nquad1, 0.0, outarray.data(),
496 nmodes0);
497 }
498}
499
500//////////////////////////
501// Evaluation functions //
502//////////////////////////
503
506{
507 eta[0] = xi[0];
508 eta[1] = xi[1];
509}
510
513{
514 xi[0] = eta[0];
515 xi[1] = eta[1];
516}
517
518/** \brief Fill outarray with mode \a mode of expansion
519 *
520 * Note for quadrilateral expansions _base[0] (i.e. p) modes run
521 * fastest
522 */
523
524void StdQuadExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
525{
526 int i;
527 int nquad0 = m_base[0]->GetNumPoints();
528 int nquad1 = m_base[1]->GetNumPoints();
529 Array<OneD, const NekDouble> base0 = m_base[0]->GetBdata();
530 Array<OneD, const NekDouble> base1 = m_base[1]->GetBdata();
531 int btmp0 = m_base[0]->GetNumModes();
532 int mode0 = mode % btmp0;
533 int mode1 = mode / btmp0;
534
535 ASSERTL2(mode1 == (int)floor((1.0 * mode) / btmp0),
536 "Integer Truncation not Equiv to Floor");
537
538 ASSERTL2(m_ncoeffs > mode,
539 "calling argument mode is larger than total expansion order");
540
541 for (i = 0; i < nquad1; ++i)
542 {
543 Vmath::Vcopy(nquad0, (NekDouble *)(base0.data() + mode0 * nquad0), 1,
544 &outarray[0] + i * nquad0, 1);
545 }
546
547 for (i = 0; i < nquad0; ++i)
548 {
549 Vmath::Vmul(nquad1, (NekDouble *)(base1.data() + mode1 * nquad1), 1,
550 &outarray[0] + i, nquad0, &outarray[0] + i, nquad0);
551 }
552}
553
554//////////////////////
555// Helper functions //
556//////////////////////
557
559{
560 return 4;
561}
562
564{
565 return 4;
566}
567
568int StdQuadExp::v_GetTraceNcoeffs(const int i) const
569{
570 ASSERTL2((i >= 0) && (i <= 3), "edge id is out of range");
571
572 if ((i == 0) || (i == 2))
573 {
574 return GetBasisNumModes(0);
575 }
576 else
577 {
578 return GetBasisNumModes(1);
579 }
580}
581
583{
584 ASSERTL2((i >= 0) && (i <= 4), "edge id is out of range");
585 if ((i == 0) || (i == 2))
586 {
587 return GetBasisNumModes(0) - 2;
588 }
589 else
590 {
591 return GetBasisNumModes(1) - 2;
592 }
593}
594
596{
597 ASSERTL2((i >= 0) && (i <= 3), "edge id is out of range");
598
599 if ((i == 0) || (i == 2))
600 {
601 return GetNumPoints(0);
602 }
603 else
604 {
605 return GetNumPoints(1);
606 }
607}
608
610 const int i, [[maybe_unused]] const int j,
611 [[maybe_unused]] bool UseGLL) const
612{
613 ASSERTL2((i >= 0) && (i <= 3), "edge id is out of range");
614
615 if ((i == 0) || (i == 2))
616 {
617 return GetBasis(0)->GetBasisKey();
618 }
619 else
620 {
621 return GetBasis(1)->GetBasisKey();
622 }
623}
624
626{
628}
629
631{
635 "BasisType is not a boundary interior form");
639 "BasisType is not a boundary interior form");
640
641 return 4 + 2 * (GetBasisNumModes(0) - 2) + 2 * (GetBasisNumModes(1) - 2);
642}
643
645{
649 "BasisType is not a boundary interior form");
653 "BasisType is not a boundary interior form");
654
655 return 2 * GetBasisNumModes(0) + 2 * GetBasisNumModes(1);
656}
657
659 const std::vector<unsigned int> &nummodes, int &modes_offset)
660{
661 int nmodes = nummodes[modes_offset] * nummodes[modes_offset + 1];
662 modes_offset += 2;
663
664 return nmodes;
665}
666
668{
669 bool returnval = false;
670
673 {
676 {
677 returnval = true;
678 }
679 }
680
681 return returnval;
682}
683
685 Array<OneD, NekDouble> &coords_1,
686 [[maybe_unused]] Array<OneD, NekDouble> &coords_2)
687{
688 Array<OneD, const NekDouble> z0 = m_base[0]->GetZ();
689 Array<OneD, const NekDouble> z1 = m_base[1]->GetZ();
690 int nq0 = GetNumPoints(0);
691 int nq1 = GetNumPoints(1);
692 int i;
693
694 for (i = 0; i < nq1; ++i)
695 {
696 Blas::Dcopy(nq0, z0.data(), 1, &coords_0[0] + i * nq0, 1);
697 Vmath::Fill(nq0, z1[i], &coords_1[0] + i * nq0, 1);
698 }
699}
700
701/**
702 * @brief This function evaluates the basis function mode @p mode at a
703 * point @p coords of the domain.
704 *
705 * This function uses barycentric interpolation with the tensor
706 * product separation of the basis function to improve performance.
707 *
708 * @param coord The coordinate inside the standard region.
709 * @param mode The mode number to be evaluated.
710 *
711 * @return The value of the basis function @p mode at @p coords.
712 */
714 const Array<OneD, const NekDouble> &coords, int mode)
715{
716 ASSERTL2(coords[0] > -1 - NekConstants::kNekZeroTol, "coord[0] < -1");
717 ASSERTL2(coords[0] < 1 + NekConstants::kNekZeroTol, "coord[0] > 1");
718 ASSERTL2(coords[1] > -1 - NekConstants::kNekZeroTol, "coord[1] < -1");
719 ASSERTL2(coords[1] < 1 + NekConstants::kNekZeroTol, "coord[1] > 1");
720
721 const int nm0 = m_base[0]->GetNumModes();
722 const int nm1 = m_base[1]->GetNumModes();
723
724 return StdExpansion::BaryEvaluateBasis<0>(coords[0], mode % nm1) *
725 StdExpansion::BaryEvaluateBasis<1>(coords[1], mode / nm0);
726}
727
729 const Array<OneD, NekDouble> &coord,
730 const Array<OneD, const NekDouble> &inarray,
731 std::array<NekDouble, 3> &firstOrderDerivs)
732{
733 return BaryTensorDeriv(coord, inarray, firstOrderDerivs);
734}
735
736//////////////
737// Mappings //
738//////////////
739
741{
742 int i;
743 int cnt = 0;
744 int nummodes0, nummodes1;
745 int value1 = 0, value2 = 0;
746 if (outarray.size() != NumBndryCoeffs())
747 {
749 }
750
751 nummodes0 = m_base[0]->GetNumModes();
752 nummodes1 = m_base[1]->GetNumModes();
753
754 const LibUtilities::BasisType Btype0 = GetBasisType(0);
755 const LibUtilities::BasisType Btype1 = GetBasisType(1);
756
757 switch (Btype1)
758 {
761 value1 = nummodes0;
762 break;
764 value1 = 2 * nummodes0;
765 break;
766 default:
767 ASSERTL0(0, "Mapping array is not defined for this expansion");
768 break;
769 }
770
771 for (i = 0; i < value1; i++)
772 {
773 outarray[i] = i;
774 }
775 cnt = value1;
776
777 switch (Btype0)
778 {
781 value2 = value1 + nummodes0 - 1;
782 break;
784 value2 = value1 + 1;
785 break;
786 default:
787 ASSERTL0(0, "Mapping array is not defined for this expansion");
788 break;
789 }
790
791 for (i = 0; i < nummodes1 - 2; i++)
792 {
793 outarray[cnt++] = value1 + i * nummodes0;
794 outarray[cnt++] = value2 + i * nummodes0;
795 }
796
797 if (Btype1 == LibUtilities::eGLL_Lagrange ||
799 {
800 for (i = nummodes0 * (nummodes1 - 1); i < GetNcoeffs(); i++)
801 {
802 outarray[cnt++] = i;
803 }
804 }
805}
806
808{
809 int i, j;
810 int cnt = 0;
811 int nummodes0, nummodes1;
812 int startvalue = 0;
813 if (outarray.size() != GetNcoeffs() - NumBndryCoeffs())
814 {
816 }
817
818 nummodes0 = m_base[0]->GetNumModes();
819 nummodes1 = m_base[1]->GetNumModes();
820
821 const LibUtilities::BasisType Btype0 = GetBasisType(0);
822 const LibUtilities::BasisType Btype1 = GetBasisType(1);
823
824 switch (Btype1)
825 {
827 startvalue = nummodes0;
828 break;
830 startvalue = 2 * nummodes0;
831 break;
832 default:
833 ASSERTL0(0, "Mapping array is not defined for this expansion");
834 break;
835 }
836
837 switch (Btype0)
838 {
840 startvalue++;
841 break;
843 startvalue += 2;
844 break;
845 default:
846 ASSERTL0(0, "Mapping array is not defined for this expansion");
847 break;
848 }
849
850 for (i = 0; i < nummodes1 - 2; i++)
851 {
852 for (j = 0; j < nummodes0 - 2; j++)
853 {
854 outarray[cnt++] = startvalue + j;
855 }
856 startvalue += nummodes0;
857 }
858}
859
860int StdQuadExp::v_GetVertexMap(int localVertexId, bool useCoeffPacking)
861{
862 int localDOF = 0;
863
864 if (useCoeffPacking == true)
865 {
866 switch (localVertexId)
867 {
868 case 0:
869 {
870 localDOF = 0;
871 }
872 break;
873 case 1:
874 {
876 {
877 localDOF = m_base[0]->GetNumModes() - 1;
878 }
879 else
880 {
881 localDOF = 1;
882 }
883 }
884 break;
885 case 2:
886 {
888 {
889 localDOF = m_base[0]->GetNumModes() *
890 (m_base[1]->GetNumModes() - 1);
891 }
892 else
893 {
894 localDOF = m_base[0]->GetNumModes();
895 }
896 }
897 break;
898 case 3:
899 {
901 {
902 localDOF =
903 m_base[0]->GetNumModes() * m_base[1]->GetNumModes() - 1;
904 }
905 else
906 {
907 localDOF = m_base[0]->GetNumModes() + 1;
908 }
909 }
910 break;
911 default:
912 ASSERTL0(false, "eid must be between 0 and 3");
913 break;
914 }
915 }
916 else
917 {
918 switch (localVertexId)
919 {
920 case 0:
921 {
922 localDOF = 0;
923 }
924 break;
925 case 1:
926 {
928 {
929 localDOF = m_base[0]->GetNumModes() - 1;
930 }
931 else
932 {
933 localDOF = 1;
934 }
935 }
936 break;
937 case 2:
938 {
940 {
941 localDOF =
942 m_base[0]->GetNumModes() * m_base[1]->GetNumModes() - 1;
943 }
944 else
945 {
946 localDOF = m_base[0]->GetNumModes() + 1;
947 }
948 }
949 break;
950 case 3:
951 {
953 {
954 localDOF = m_base[0]->GetNumModes() *
955 (m_base[1]->GetNumModes() - 1);
956 }
957 else
958 {
959 localDOF = m_base[0]->GetNumModes();
960 }
961 }
962 break;
963 default:
964 ASSERTL0(false, "eid must be between 0 and 3");
965 break;
966 }
967 }
968 return localDOF;
969}
970
971/** \brief Get the map of the coefficient location to teh
972 * local trace coefficients
973 */
974
975void StdQuadExp::v_GetTraceCoeffMap(const unsigned int traceid,
977{
978 ASSERTL1(traceid < 4, "traceid must be between 0 and 3");
979
980 unsigned int i;
981 unsigned int order0 = m_base[0]->GetNumModes();
982 unsigned int order1 = m_base[1]->GetNumModes();
983 unsigned int numModes = (traceid % 2) ? order1 : order0;
984
985 if (maparray.size() != numModes)
986 {
987 maparray = Array<OneD, unsigned int>(numModes);
988 }
989
990 const LibUtilities::BasisType bType = GetBasisType(traceid % 2);
991
992 if (bType == LibUtilities::eModified_A)
993 {
994 switch (traceid)
995 {
996 case 0:
997 {
998 for (i = 0; i < numModes; i++)
999 {
1000 maparray[i] = i;
1001 }
1002 }
1003 break;
1004 case 1:
1005 {
1006 for (i = 0; i < numModes; i++)
1007 {
1008 maparray[i] = i * order0 + 1;
1009 }
1010 }
1011 break;
1012 case 2:
1013 {
1014 for (i = 0; i < numModes; i++)
1015 {
1016 maparray[i] = order0 + i;
1017 }
1018 }
1019 break;
1020 case 3:
1021 {
1022 for (i = 0; i < numModes; i++)
1023 {
1024 maparray[i] = i * order0;
1025 }
1026 }
1027 break;
1028 default:
1029 break;
1030 }
1031 }
1032 else if (bType == LibUtilities::eGLL_Lagrange ||
1034 {
1035 switch (traceid)
1036 {
1037 case 0:
1038 {
1039 for (i = 0; i < numModes; i++)
1040 {
1041 maparray[i] = i;
1042 }
1043 }
1044 break;
1045 case 1:
1046 {
1047 for (i = 0; i < numModes; i++)
1048 {
1049 maparray[i] = (i + 1) * order0 - 1;
1050 }
1051 }
1052 break;
1053 case 2:
1054 {
1055 for (i = 0; i < numModes; i++)
1056 {
1057 maparray[i] = order0 * (order1 - 1) + i;
1058 }
1059 }
1060 break;
1061 case 3:
1062 {
1063 for (i = 0; i < numModes; i++)
1064 {
1065 maparray[i] = order0 * i;
1066 }
1067 }
1068 break;
1069 default:
1070 break;
1071 }
1072 }
1073 else
1074 {
1075 ASSERTL0(false, "Mapping not defined for this type of basis");
1076 }
1077}
1078
1080 const int eid, Array<OneD, unsigned int> &maparray,
1081 Array<OneD, int> &signarray, const Orientation edgeOrient)
1082{
1083 int i;
1084 const int nummodes0 = m_base[0]->GetNumModes();
1085 const int nummodes1 = m_base[1]->GetNumModes();
1086 const int nEdgeIntCoeffs = GetTraceNcoeffs(eid) - 2;
1087 const LibUtilities::BasisType bType = GetBasisType(eid % 2);
1088
1089 if (maparray.size() != nEdgeIntCoeffs)
1090 {
1091 maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1092 }
1093
1094 if (signarray.size() != nEdgeIntCoeffs)
1095 {
1096 signarray = Array<OneD, int>(nEdgeIntCoeffs, 1);
1097 }
1098 else
1099 {
1100 fill(signarray.data(), signarray.data() + nEdgeIntCoeffs, 1);
1101 }
1102
1103 if (bType == LibUtilities::eModified_A)
1104 {
1105 switch (eid)
1106 {
1107 case 0:
1108 {
1109 for (i = 0; i < nEdgeIntCoeffs; i++)
1110 {
1111 maparray[i] = i + 2;
1112 }
1113 }
1114 break;
1115 case 1:
1116 {
1117 for (i = 0; i < nEdgeIntCoeffs; i++)
1118 {
1119 maparray[i] = (i + 2) * nummodes0 + 1;
1120 }
1121 }
1122 break;
1123 case 2:
1124 {
1125 for (i = 0; i < nEdgeIntCoeffs; i++)
1126 {
1127 maparray[i] = nummodes0 + i + 2;
1128 }
1129 }
1130 break;
1131 case 3:
1132 {
1133 for (i = 0; i < nEdgeIntCoeffs; i++)
1134 {
1135 maparray[i] = (i + 2) * nummodes0;
1136 }
1137 }
1138 break;
1139 default:
1140 ASSERTL0(false, "eid must be between 0 and 3");
1141 break;
1142 }
1143
1144 if (edgeOrient == eBackwards)
1145 {
1146 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1147 {
1148 signarray[i] = -1;
1149 }
1150 }
1151 }
1152 else if (bType == LibUtilities::eGLL_Lagrange)
1153 {
1154 switch (eid)
1155 {
1156 case 0:
1157 {
1158 for (i = 0; i < nEdgeIntCoeffs; i++)
1159 {
1160 maparray[i] = i + 1;
1161 }
1162 }
1163 break;
1164 case 1:
1165 {
1166 for (i = 0; i < nEdgeIntCoeffs; i++)
1167 {
1168 maparray[i] = (i + 2) * nummodes0 - 1;
1169 }
1170 }
1171 break;
1172 case 2:
1173 {
1174 for (i = 0; i < nEdgeIntCoeffs; i++)
1175 {
1176 maparray[i] = nummodes0 * (nummodes1 - 1) + i + 1;
1177 }
1178 }
1179 break;
1180 case 3:
1181 {
1182 for (i = 0; i < nEdgeIntCoeffs; i++)
1183 {
1184 maparray[i] = nummodes0 * (i + 1);
1185 }
1186 }
1187 break;
1188 default:
1189 ASSERTL0(false, "eid must be between 0 and 3");
1190 break;
1191 }
1192 if (edgeOrient == eBackwards)
1193 {
1194 reverse(maparray.data(), maparray.data() + nEdgeIntCoeffs);
1195 }
1196 }
1197 else
1198 {
1199 ASSERTL0(false, "Mapping not defined for this type of basis");
1200 }
1201}
1202
1203///////////////////////
1204// Wrapper Functions //
1205///////////////////////
1206
1208{
1209 int i, j;
1210 int order0 = GetBasisNumModes(0);
1211 int order1 = GetBasisNumModes(1);
1212 MatrixType mtype = mkey.GetMatrixType();
1213
1214 DNekMatSharedPtr Mat;
1215
1216 switch (mtype)
1217 {
1219 {
1220 int nq0 = m_base[0]->GetNumPoints();
1221 int nq1 = m_base[1]->GetNumPoints();
1222 int nq;
1223
1224 // take definition from key
1226 {
1227 nq = (int)mkey.GetConstFactor(eFactorConst);
1228 }
1229 else
1230 {
1231 nq = max(nq0, nq1);
1232 }
1233
1234 int neq =
1237 Array<OneD, NekDouble> coll(2);
1239 Array<OneD, NekDouble> tmp(nq0);
1240
1241 Mat = MemoryManager<DNekMat>::AllocateSharedPtr(neq, nq0 * nq1);
1242 int cnt = 0;
1243
1244 for (i = 0; i < nq; ++i)
1245 {
1246 for (j = 0; j < nq; ++j, ++cnt)
1247 {
1248 coords[cnt] = Array<OneD, NekDouble>(2);
1249 coords[cnt][0] = -1.0 + 2 * j / (NekDouble)(nq - 1);
1250 coords[cnt][1] = -1.0 + 2 * i / (NekDouble)(nq - 1);
1251 }
1252 }
1253
1254 for (i = 0; i < neq; ++i)
1255 {
1256 LocCoordToLocCollapsed(coords[i], coll);
1257
1258 I[0] = m_base[0]->GetI(coll);
1259 I[1] = m_base[1]->GetI(coll + 1);
1260
1261 // interpolate first coordinate direction
1262 for (j = 0; j < nq1; ++j)
1263 {
1264 NekDouble fac = (I[1]->GetPtr())[j];
1265 Vmath::Smul(nq0, fac, I[0]->GetPtr(), 1, tmp, 1);
1266
1267 Vmath::Vcopy(nq0, &tmp[0], 1,
1268 Mat->GetRawPtr() + j * nq0 * neq + i, neq);
1269 }
1270 }
1271 break;
1272 }
1273 case eMass:
1274 {
1276 // For Fourier basis set the imaginary component of mean mode
1277 // to have a unit diagonal component in mass matrix
1279 {
1280 for (i = 0; i < order1; ++i)
1281 {
1282 (*Mat)(order0 *i + 1, i * order0 + 1) = 1.0;
1283 }
1284 }
1285
1287 {
1288 for (i = 0; i < order0; ++i)
1289 {
1290 (*Mat)(order0 + i, order0 + i) = 1.0;
1291 }
1292 }
1293 break;
1294 }
1295 case eFwdTrans:
1296 {
1297 Mat =
1299 StdMatrixKey iprodkey(eIProductWRTBase, DetShapeType(), *this);
1300 DNekMat &Iprod = *GetStdMatrix(iprodkey);
1301 StdMatrixKey imasskey(eInvMass, DetShapeType(), *this);
1302 DNekMat &Imass = *GetStdMatrix(imasskey);
1303
1304 (*Mat) = Imass * Iprod;
1305 break;
1306 }
1307 case eGaussDG:
1308 {
1310
1311 int edge = (int)factors[StdRegions::eFactorGaussEdge];
1312 int dir = (edge + 1) % 2;
1313 int nCoeffs = m_base[dir]->GetNumModes();
1314
1315 const LibUtilities::PointsKey BS_p(
1318 nCoeffs, BS_p);
1319
1320 Array<OneD, NekDouble> coords(1, 0.0);
1321 coords[0] = (edge == 0 || edge == 3) ? -1.0 : 1.0;
1322
1325 DNekMatSharedPtr m_Ix = basis->GetI(coords);
1326
1328 Vmath::Vcopy(nCoeffs, m_Ix->GetPtr(), 1, Mat->GetPtr(), 1);
1329 break;
1330 }
1331 default:
1332 {
1334 break;
1335 }
1336 }
1337
1338 return Mat;
1339}
1340
1342{
1343 return GenMatrix(mkey);
1344}
1345
1346///////////////////////////////////
1347// Operator evaluation functions //
1348///////////////////////////////////
1349
1351 const StdMatrixKey &mkey)
1352{
1353 // Generate an orthonogal expansion
1354 int qa = m_base[0]->GetNumPoints();
1355 int qb = m_base[1]->GetNumPoints();
1356 int nmodes_a = m_base[0]->GetNumModes();
1357 int nmodes_b = m_base[1]->GetNumModes();
1358 int nmodes = min(nmodes_a, nmodes_b);
1359 // Declare orthogonal basis.
1362
1365 StdQuadExp OrthoExp(Ba, Bb);
1366
1367 // SVV parameters loaded from the .xml case file
1368 Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1369
1370 // project onto modal space.
1371 OrthoExp.FwdTrans(array, orthocoeffs);
1372
1373 if (mkey.ConstFactorExists(
1374 eFactorSVVPowerKerDiffCoeff)) // Rodrigo's power kernel
1375 {
1377 NekDouble SvvDiffCoeff =
1380
1381 for (int j = 0; j < nmodes_a; ++j)
1382 {
1383 for (int k = 0; k < nmodes_b; ++k)
1384 {
1385 // linear space but makes high modes very negative
1386 NekDouble fac = std::max(
1387 pow((1.0 * j) / (nmodes_a - 1), cutoff * nmodes_a),
1388 pow((1.0 * k) / (nmodes_b - 1), cutoff * nmodes_b));
1389
1390 orthocoeffs[j * nmodes_b + k] *= SvvDiffCoeff * fac;
1391 }
1392 }
1393 }
1394 else if (mkey.ConstFactorExists(
1395 eFactorSVVDGKerDiffCoeff)) // Rodrigo/mansoor's DG kernel
1396 {
1399 int max_ab = max(nmodes_a - kSVVDGFiltermodesmin,
1400 nmodes_b - kSVVDGFiltermodesmin);
1401 max_ab = max(max_ab, 0);
1402 max_ab = min(max_ab, kSVVDGFiltermodesmax - kSVVDGFiltermodesmin);
1403
1404 for (int j = 0; j < nmodes_a; ++j)
1405 {
1406 for (int k = 0; k < nmodes_b; ++k)
1407 {
1408 int maxjk = max(j, k);
1409 maxjk = min(maxjk, kSVVDGFiltermodesmax - 1);
1410
1411 orthocoeffs[j * nmodes_b + k] *=
1412 SvvDiffCoeff * kSVVDGFilter[max_ab][maxjk];
1413 }
1414 }
1415 }
1416 else
1417 {
1418 NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
1419 // Exponential Kernel implementation
1420 int cutoff = (int)(mkey.GetConstFactor(eFactorSVVCutoffRatio) *
1421 min(nmodes_a, nmodes_b));
1422
1423 //------"New" Version August 22nd '13--------------------
1424 for (int j = 0; j < nmodes_a; ++j)
1425 {
1426 for (int k = 0; k < nmodes_b; ++k)
1427 {
1428 if (j + k >= cutoff) // to filter out only the "high-modes"
1429 {
1430 orthocoeffs[j * nmodes_b + k] *=
1431 (SvvDiffCoeff *
1432 exp(-(j + k - nmodes) * (j + k - nmodes) /
1433 ((NekDouble)((j + k - cutoff + 1) *
1434 (j + k - cutoff + 1)))));
1435 }
1436 else
1437 {
1438 orthocoeffs[j * nmodes_b + k] *= 0.0;
1439 }
1440 }
1441 }
1442 }
1443
1444 // backward transform to physical space
1445 OrthoExp.BwdTrans(orthocoeffs, array);
1446}
1447
1449 const NekDouble alpha,
1450 const NekDouble exponent,
1451 const NekDouble cutoff)
1452{
1453 // Generate an orthogonal expansion
1454 int qa = m_base[0]->GetNumPoints();
1455 int qb = m_base[1]->GetNumPoints();
1456 int nmodesA = m_base[0]->GetNumModes();
1457 int nmodesB = m_base[1]->GetNumModes();
1458 int P = nmodesA - 1;
1459 int Q = nmodesB - 1;
1460
1461 // Declare orthogonal basis.
1464
1467 StdQuadExp OrthoExp(Ba, Bb);
1468
1469 // Cutoff
1470 int Pcut = cutoff * P;
1471 int Qcut = cutoff * Q;
1472
1473 // Project onto orthogonal space.
1474 Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1475 OrthoExp.FwdTrans(array, orthocoeffs);
1476
1477 //
1478 NekDouble fac, fac1, fac2;
1479 for (int i = 0; i < nmodesA; ++i)
1480 {
1481 for (int j = 0; j < nmodesB; ++j)
1482 {
1483 // to filter out only the "high-modes"
1484 if (i > Pcut || j > Qcut)
1485 {
1486 fac1 = (NekDouble)(i - Pcut) / ((NekDouble)(P - Pcut));
1487 fac2 = (NekDouble)(j - Qcut) / ((NekDouble)(Q - Qcut));
1488 fac = max(fac1, fac2);
1489 fac = pow(fac, exponent);
1490 orthocoeffs[i * nmodesB + j] *= exp(-alpha * fac);
1491 }
1492 }
1493 }
1494
1495 // backward transform to physical space
1496 OrthoExp.BwdTrans(orthocoeffs, array);
1497}
1498
1500 int numMin, const Array<OneD, const NekDouble> &inarray,
1501 Array<OneD, NekDouble> &outarray)
1502{
1503 int n_coeffs = inarray.size();
1504
1505 Array<OneD, NekDouble> coeff(n_coeffs);
1506 Array<OneD, NekDouble> coeff_tmp(n_coeffs, 0.0);
1509
1510 int nmodes0 = m_base[0]->GetNumModes();
1511 int nmodes1 = m_base[1]->GetNumModes();
1512 int numMax = nmodes0;
1513
1514 Vmath::Vcopy(n_coeffs, inarray, 1, coeff_tmp, 1);
1515
1516 const LibUtilities::PointsKey Pkey0(nmodes0,
1518 const LibUtilities::PointsKey Pkey1(nmodes1,
1520
1521 LibUtilities::BasisKey b0(m_base[0]->GetBasisType(), nmodes0, Pkey0);
1522 LibUtilities::BasisKey b1(m_base[1]->GetBasisType(), nmodes1, Pkey1);
1523
1524 LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A, nmodes0, Pkey0);
1525 LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_A, nmodes1, Pkey1);
1526
1527 LibUtilities::InterpCoeff2D(b0, b1, coeff_tmp, bortho0, bortho1, coeff);
1528
1529 Vmath::Zero(n_coeffs, coeff_tmp, 1);
1530
1531 int cnt = 0;
1532 for (int i = 0; i < numMin + 1; ++i)
1533 {
1534 Vmath::Vcopy(numMin, tmp = coeff + cnt, 1, tmp2 = coeff_tmp + cnt, 1);
1535
1536 cnt = i * numMax;
1537 }
1538
1539 LibUtilities::InterpCoeff2D(bortho0, bortho1, coeff_tmp, b0, b1, outarray);
1540}
1541
1543 Array<OneD, NekDouble> &outarray,
1544 const StdMatrixKey &mkey)
1545{
1546 StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
1547}
1548
1550 const Array<OneD, const NekDouble> &inarray,
1551 Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
1552{
1553 StdQuadExp::v_LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
1554}
1555
1557 const int k1, const int k2, const Array<OneD, const NekDouble> &inarray,
1558 Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
1559{
1560 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1561}
1562
1564 const int i, const Array<OneD, const NekDouble> &inarray,
1565 Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
1566{
1567 StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
1568}
1569
1571 const Array<OneD, const NekDouble> &inarray,
1572 Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
1573{
1574 StdQuadExp::v_HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
1575}
1576
1577// up to here
1579 const Array<OneD, const NekDouble> &inarray,
1580 Array<OneD, NekDouble> &outarray)
1581{
1582 int i;
1583 int nquad0 = m_base[0]->GetNumPoints();
1584 int nquad1 = m_base[1]->GetNumPoints();
1585
1586 const Array<OneD, const NekDouble> &w0 = m_base[0]->GetW();
1587 const Array<OneD, const NekDouble> &w1 = m_base[1]->GetW();
1588
1589 // multiply by integration constants
1590 for (i = 0; i < nquad1; ++i)
1591 {
1592 Vmath::Vmul(nquad0, inarray.data() + i * nquad0, 1, w0.data(), 1,
1593 outarray.data() + i * nquad0, 1);
1594 }
1595
1596 for (i = 0; i < nquad0; ++i)
1597 {
1598 Vmath::Vmul(nquad1, outarray.data() + i, nquad0, w1.data(), 1,
1599 outarray.data() + i, nquad0);
1600 }
1601}
1602
1604 Array<OneD, int> &conn, [[maybe_unused]] bool standard)
1605{
1606 int np1 = m_base[0]->GetNumPoints();
1607 int np2 = m_base[1]->GetNumPoints();
1608 int np = max(np1, np2);
1609
1610 conn = Array<OneD, int>(6 * (np - 1) * (np - 1));
1611
1612 int row = 0;
1613 int rowp1 = 0;
1614 int cnt = 0;
1615 for (int i = 0; i < np - 1; ++i)
1616 {
1617 rowp1 += np;
1618 for (int j = 0; j < np - 1; ++j)
1619 {
1620 conn[cnt++] = row + j;
1621 conn[cnt++] = row + j + 1;
1622 conn[cnt++] = rowp1 + j;
1623
1624 conn[cnt++] = rowp1 + j + 1;
1625 conn[cnt++] = rowp1 + j;
1626 conn[cnt++] = row + j + 1;
1627 }
1628 row += np;
1629 }
1630}
1631
1632} // 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
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:47
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 PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d0, Array< OneD, NekDouble > &outarray_d1)
Calculate the 2D derivative in the local tensor/collapsed coordinate at the physical points.
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)
NekDouble Integral(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &w0, const Array< OneD, const NekDouble > &w1)
void BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0=true, bool doCheckCollDir1=true)
void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0=true, bool doCheckCollDir1=true)
The base class for all shapes.
Definition: StdExpansion.h:65
const LibUtilities::BasisSharedPtr & GetBasis(int dir) const
This function gets the shared point to basis in the dir direction.
Definition: StdExpansion.h:112
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 MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:761
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)
void GetTraceToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards, int P=-1, int Q=-1)
Definition: StdExpansion.h:693
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:370
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:683
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space.
Definition: StdExpansion.h:433
int GetTraceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th trace.
Definition: StdExpansion.h:261
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:853
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.
LibUtilities::NekManager< StdMatrixKey, DNekBlkMat, StdMatrixKey::opLess > m_stdStaticCondMatrixManager
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
Array< OneD, LibUtilities::BasisSharedPtr > m_base
void MassMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:83
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:138
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:124
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:133
void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Transform a given function from physical quadrature space to coefficient space.
Definition: StdQuadExp.cpp:204
void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdQuadExp.cpp:141
void v_ExponentialFilter(Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff) override
void v_BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1) override
Definition: StdQuadExp.cpp:161
int v_NumDGBndryCoeffs() const final
Definition: StdQuadExp.cpp:644
void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdQuadExp.cpp:227
void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta) override
Definition: StdQuadExp.cpp:504
int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false) override
Definition: StdQuadExp.cpp:860
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdQuadExp.cpp:402
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
bool v_IsBoundaryInteriorExpansion() const override
Definition: StdQuadExp.cpp:667
int v_GetTraceNcoeffs(const int i) const final
Definition: StdQuadExp.cpp:568
int v_GetNtraces() const final
Definition: StdQuadExp.cpp:563
void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
int v_GetTraceNumPoints(const int i) const final
Definition: StdQuadExp.cpp:595
void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
int v_NumBndryCoeffs() const final
Definition: StdQuadExp.cpp:630
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey) override
void v_GetCoords(Array< OneD, NekDouble > &coords_0, Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2) override
Definition: StdQuadExp.cpp:684
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) override
This function evaluates the basis function mode mode at a point coords of the domain.
Definition: StdQuadExp.cpp:713
int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset) override
Definition: StdQuadExp.cpp:658
void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculate the inner product of inarray with respect to the basis B=base0*base1 and put into outarray.
Definition: StdQuadExp.cpp:354
void v_GetInteriorMap(Array< OneD, unsigned int > &outarray) override
Definition: StdQuadExp.cpp:807
int v_GetNverts() const final
Definition: StdQuadExp.cpp:558
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray) override
Integrates the specified function over the domain.
Definition: StdQuadExp.cpp:58
void v_StdPhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray) override
Definition: StdQuadExp.cpp:108
void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdQuadExp.cpp:127
void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi) override
Definition: StdQuadExp.cpp:511
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
Definition: StdQuadExp.cpp:728
void v_GetTraceCoeffMap(const unsigned int traceid, Array< OneD, unsigned int > &maparray) override
Get the map of the coefficient location to teh local trace coefficients.
Definition: StdQuadExp.cpp:975
const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int j, bool UseGLL=false) const final
Definition: StdQuadExp.cpp:609
DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey) override
void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true) override
void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
Definition: StdQuadExp.cpp:367
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdQuadExp.cpp:395
void v_IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1) override
Definition: StdQuadExp.cpp:443
void v_GetTraceInteriorToElementMap(const int eid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation edgeOrient=eForwards) override
void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray) override
Calculate the derivative of the physical points.
Definition: StdQuadExp.cpp:76
void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray) override
Definition: StdQuadExp.cpp:740
DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey) override
int v_GetTraceIntNcoeffs(const int i) const final
Definition: StdQuadExp.cpp:582
StdQuadExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb)
Constructor using BasisKey class for quadrature points and order definition.
Definition: StdQuadExp.cpp:47
void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
LibUtilities::ShapeType v_DetShapeType() const final
Definition: StdQuadExp.cpp:625
void v_FillMode(const int mode, Array< OneD, NekDouble > &array) override
Fill outarray with mode mode of expansion.
Definition: StdQuadExp.cpp:524
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = alpha A x plus beta y where A[m x n].
Definition: Blas.hpp:211
static void Dcopy(const int &n, const double *x, const int &incx, double *y, const int &incy)
BLAS level 1: Copy x to y.
Definition: Blas.hpp:128
static double Ddot(const int &n, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: output = .
Definition: Blas.hpp:163
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition: Blas.hpp:383
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:133
BasisManagerT & BasisManager(void)
std::shared_ptr< Basis > BasisSharedPtr
void InterpCoeff2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
Definition: InterpCoeff.cpp:67
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:46
@ P
Monomial polynomials .
Definition: BasisType.h:62
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
Definition: BasisType.h:57
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:42
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:56
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
@ eFourier
Fourier Expansion .
Definition: BasisType.h:55
static const NekDouble kNekZeroTol
const int kSVVDGFiltermodesmin
Definition: StdRegions.hpp:500
const int kSVVDGFiltermodesmax
Definition: StdRegions.hpp:501
const NekDouble kSVVDGFilter[9][11]
Definition: StdRegions.hpp:503
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:430
std::shared_ptr< StdSegExp > StdSegExpSharedPtr
Definition: StdSegExp.h:217
StdRegions::ConstFactorMap factors
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 Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.hpp:220
STL namespace.