Nektar++
StdTriExp.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: StdTriExp.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: Triangle routines built upon StdExpansion2D
32//
33///////////////////////////////////////////////////////////////////////////////
34
36
37using namespace std;
38
39namespace Nektar::StdRegions
40{
42 const LibUtilities::BasisKey &Bb)
43 : StdExpansion(LibUtilities::StdTriData::getNumberOfCoefficients(
44 Ba.GetNumModes(), Bb.GetNumModes()),
45 2, Ba, Bb),
46 StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(
47 Ba.GetNumModes(), Bb.GetNumModes()),
48 Ba, Bb)
49{
50 ASSERTL0(Ba.GetNumModes() <= Bb.GetNumModes(),
51 "order in 'a' direction is higher than order "
52 "in 'b' direction");
53}
54
55//-------------------------------
56// Integration Methods
57//-------------------------------
59{
60 int i;
61 int nquad1 = m_base[1]->GetNumPoints();
62 Array<OneD, NekDouble> w1_tmp(nquad1);
63
67
68 switch (m_base[1]->GetPointsType())
69 {
70 case LibUtilities::eGaussRadauMAlpha1Beta0: // (0,1) Jacobi Inner
71 // product
72 {
73 Vmath::Smul(nquad1, 0.5, w1, 1, w1_tmp, 1);
74 break;
75 }
76 default:
77 {
78 // include jacobian factor on whatever coordinates are defined.
79 for (i = 0; i < nquad1; ++i)
80 {
81 w1_tmp[i] = 0.5 * (1 - z1[i]) * w1[i];
82 }
83 break;
84 }
85 }
86
87 return StdExpansion2D::Integral(inarray, w0, w1_tmp);
88}
89
90//-----------------------------
91// Differentiation Methods
92//-----------------------------
93
94/**
95 * \brief Calculate the derivative of the physical points.
96 *
97 * \f$ \frac{\partial u}{\partial x_1} = \left .
98 * \frac{2.0}{1-\eta_2} \frac{\partial u}{\partial d\eta_1}
99 * \right |_{\eta_2}\f$
100 *
101 * \f$ \frac{\partial u}{\partial x_2} = \left .
102 * \frac{1+\eta_1}{1-\eta_2} \frac{\partial u}{\partial d\eta_1}
103 * \right |_{\eta_2} + \left . \frac{\partial u}{\partial d\eta_2}
104 * \right |_{\eta_1} \f$
105 */
109 [[maybe_unused]] Array<OneD, NekDouble> &out_d2)
110{
111 int i;
112 int nquad0 = m_base[0]->GetNumPoints();
113 int nquad1 = m_base[1]->GetNumPoints();
114 Array<OneD, NekDouble> wsp(std::max(nquad0, nquad1));
115
116 const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
117 const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
118
119 // set up geometric factor: 2/(1-z1)
120 Vmath::Sadd(nquad1, -1.0, z1, 1, wsp, 1);
121 Vmath::Sdiv(nquad1, -2.0, wsp, 1, wsp, 1);
122
123 if (out_d0.size() > 0)
124 {
125 PhysTensorDeriv(inarray, out_d0, out_d1);
126
127 for (i = 0; i < nquad1; ++i)
128 {
129 Blas::Dscal(nquad0, wsp[i], &out_d0[0] + i * nquad0, 1);
130 }
131
132 // if no d1 required do not need to calculate both deriv
133 if (out_d1.size() > 0)
134 {
135 // set up geometric factor: (1+z0)/2
136 Vmath::Sadd(nquad0, 1.0, z0, 1, wsp, 1);
137 Vmath::Smul(nquad0, 0.5, wsp, 1, wsp, 1);
138
139 for (i = 0; i < nquad1; ++i)
140 {
141 Vmath::Vvtvp(nquad0, &wsp[0], 1, &out_d0[0] + i * nquad0, 1,
142 &out_d1[0] + i * nquad0, 1,
143 &out_d1[0] + i * nquad0, 1);
144 }
145 }
146 }
147 else if (out_d1.size() > 0)
148 {
149 Array<OneD, NekDouble> diff0(nquad0 * nquad1);
150 PhysTensorDeriv(inarray, diff0, out_d1);
151
152 for (i = 0; i < nquad1; ++i)
153 {
154 Blas::Dscal(nquad0, wsp[i], &diff0[0] + i * nquad0, 1);
155 }
156
157 Vmath::Sadd(nquad0, 1.0, z0, 1, wsp, 1);
158 Vmath::Smul(nquad0, 0.5, wsp, 1, wsp, 1);
159
160 for (i = 0; i < nquad1; ++i)
161 {
162 Vmath::Vvtvp(nquad0, &wsp[0], 1, &diff0[0] + i * nquad0, 1,
163 &out_d1[0] + i * nquad0, 1, &out_d1[0] + i * nquad0,
164 1);
165 }
166 }
167}
168
169void StdTriExp::v_PhysDeriv(const int dir,
170 const Array<OneD, const NekDouble> &inarray,
171 Array<OneD, NekDouble> &outarray)
172{
173 switch (dir)
174 {
175 case 0:
176 {
177 v_PhysDeriv(inarray, outarray, NullNekDouble1DArray);
178 break;
179 }
180 case 1:
181 {
182 v_PhysDeriv(inarray, NullNekDouble1DArray, outarray);
183 break;
184 }
185 default:
186 {
187 ASSERTL1(false, "input dir is out of range");
188 break;
189 }
190 }
191}
192
196 [[maybe_unused]] Array<OneD, NekDouble> &out_d2)
197{
198 StdTriExp::v_PhysDeriv(inarray, out_d0, out_d1);
199}
200
201void StdTriExp::v_StdPhysDeriv(const int dir,
202 const Array<OneD, const NekDouble> &inarray,
203 Array<OneD, NekDouble> &outarray)
204{
205 StdTriExp::v_PhysDeriv(dir, inarray, outarray);
206}
207
208//---------------------------------------
209// Transforms
210//---------------------------------------
211
212/**
213 * \brief Backward tranform for triangular elements
214 *
215 * @note 'q' (base[1]) runs fastest in this element.
216 */
218 Array<OneD, NekDouble> &outarray)
219{
220 v_BwdTrans_SumFac(inarray, outarray);
221}
222
224 Array<OneD, NekDouble> &outarray)
225{
227 m_base[0]->GetNumModes());
228
229 BwdTrans_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(), inarray,
230 outarray, wsp);
231}
232
234 const Array<OneD, const NekDouble> &base0,
235 const Array<OneD, const NekDouble> &base1,
236 const Array<OneD, const NekDouble> &inarray,
238 [[maybe_unused]] bool doCheckCollDir0,
239 [[maybe_unused]] bool doCheckCollDir1)
240{
241 int i;
242 int mode;
243 int nquad0 = m_base[0]->GetNumPoints();
244 int nquad1 = m_base[1]->GetNumPoints();
245 int nmodes0 = m_base[0]->GetNumModes();
246 int nmodes1 = m_base[1]->GetNumModes();
247
248 ASSERTL1(wsp.size() >= nquad1 * nmodes0,
249 "Workspace size is not sufficient");
252 "Basis[1] is not of general tensor type");
253
254 for (i = mode = 0; i < nmodes0; ++i)
255 {
256 Blas::Dgemv('N', nquad1, nmodes1 - i, 1.0, base1.data() + mode * nquad1,
257 nquad1, &inarray[0] + mode, 1, 0.0, &wsp[0] + i * nquad1,
258 1);
259 mode += nmodes1 - i;
260 }
261
262 // fix for modified basis by splitting top vertex mode
264 {
265 Blas::Daxpy(nquad1, inarray[1], base1.data() + nquad1, 1,
266 &wsp[0] + nquad1, 1);
267 }
268
269 Blas::Dgemm('N', 'T', nquad0, nquad1, nmodes0, 1.0, base0.data(), nquad0,
270 &wsp[0], nquad1, 0.0, &outarray[0], nquad0);
271}
272
274 Array<OneD, NekDouble> &outarray)
275{
276 v_IProductWRTBase(inarray, outarray);
277
278 // get Mass matrix inverse
279 StdMatrixKey masskey(eInvMass, DetShapeType(), *this);
280 DNekMatSharedPtr matsys = GetStdMatrix(masskey);
281
282 // copy inarray in case inarray == outarray
283 NekVector<NekDouble> in(m_ncoeffs, outarray, eCopy);
285
286 out = (*matsys) * in;
287}
288
290 const Array<OneD, const NekDouble> &inarray,
291 Array<OneD, NekDouble> &outarray)
292{
293 int i, j;
294 int npoints[2] = {m_base[0]->GetNumPoints(), m_base[1]->GetNumPoints()};
295 int nmodes[2] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes()};
296
297 fill(outarray.data(), outarray.data() + m_ncoeffs, 0.0);
298
299 Array<OneD, NekDouble> physEdge[3];
300 Array<OneD, NekDouble> coeffEdge[3];
301 for (i = 0; i < 3; i++)
302 {
303 physEdge[i] = Array<OneD, NekDouble>(npoints[i != 0]);
304 coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i != 0]);
305 }
306
307 for (i = 0; i < npoints[0]; i++)
308 {
309 physEdge[0][i] = inarray[i];
310 }
311
312 for (i = 0; i < npoints[1]; i++)
313 {
314 physEdge[1][i] = inarray[npoints[0] - 1 + i * npoints[0]];
315 physEdge[2][i] =
316 inarray[(npoints[1] - 1) * npoints[0] - i * npoints[0]];
317 }
318
319 StdSegExpSharedPtr segexp[2] = {
321 m_base[0]->GetBasisKey()),
323 m_base[1]->GetBasisKey())};
324
326 Array<OneD, int> signArray;
328
329 for (i = 0; i < 3; i++)
330 {
331 segexp[i != 0]->FwdTransBndConstrained(physEdge[i], coeffEdge[i]);
332
333 GetTraceToElementMap(i, mapArray, signArray);
334 for (j = 0; j < nmodes[i != 0]; j++)
335 {
336 sign = (NekDouble)signArray[j];
337 outarray[mapArray[j]] = sign * coeffEdge[i][j];
338 }
339 }
340
343
344 StdMatrixKey masskey(eMass, DetShapeType(), *this);
345 MassMatrixOp(outarray, tmp0, masskey);
346 v_IProductWRTBase(inarray, tmp1);
347
348 Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
349
350 // get Mass matrix inverse (only of interior DOF)
351 // use block (1,1) of the static condensed system
352 // note: this block alreay contains the inverse matrix
353 DNekMatSharedPtr matsys =
354 (m_stdStaticCondMatrixManager[masskey])->GetBlock(1, 1);
355
356 int nBoundaryDofs = v_NumBndryCoeffs();
357 int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
358
359 Array<OneD, NekDouble> rhs(nInteriorDofs);
360 Array<OneD, NekDouble> result(nInteriorDofs);
361
362 v_GetInteriorMap(mapArray);
363
364 for (i = 0; i < nInteriorDofs; i++)
365 {
366 rhs[i] = tmp1[mapArray[i]];
367 }
368
369 Blas::Dgemv('N', nInteriorDofs, nInteriorDofs, 1.0, &(matsys->GetPtr())[0],
370 nInteriorDofs, rhs.data(), 1, 0.0, result.data(), 1);
371
372 for (i = 0; i < nInteriorDofs; i++)
373 {
374 outarray[mapArray[i]] = result[i];
375 }
376}
377
378//---------------------------------------
379// Inner product functions
380//---------------------------------------
381
382/**
383 * \brief Calculate the inner product of inarray with respect to the
384 * basis B=base0[p]*base1[pq] and put into outarray.
385 *
386 * \f$
387 * \begin{array}{rcl}
388 * I_{pq} = (\phi^A_q \phi^B_{pq}, u) &=&
389 * \sum_{i=0}^{nq_0}\sum_{j=0}^{nq_1}
390 * \phi^A_p(\eta_{0,i})\phi^B_{pq}(\eta_{1,j}) w^0_i w^1_j
391 * u(\xi_{0,i} \xi_{1,j}) \\
392 * & = & \sum_{i=0}^{nq_0} \phi^A_p(\eta_{0,i})
393 * \sum_{j=0}^{nq_1} \phi^B_{pq}(\eta_{1,j}) \tilde{u}_{i,j}
394 * \end{array}
395 * \f$
396 *
397 * where
398 *
399 * \f$ \tilde{u}_{i,j} = w^0_i w^1_j u(\xi_{0,i},\xi_{1,j}) \f$
400 *
401 * which can be implemented as
402 *
403 * \f$ f_{pj} = \sum_{i=0}^{nq_0} \phi^A_p(\eta_{0,i})
404 * \tilde{u}_{i,j}
405 * \rightarrow {\bf B_1 U} \f$
406 * \f$ I_{pq} = \sum_{j=0}^{nq_1} \phi^B_{pq}(\eta_{1,j}) f_{pj}
407 * \rightarrow {\bf B_2[p*skip] f[skip]} \f$
408 *
409 * \b Recall: \f$ \eta_{1} = \frac{2(1+\xi_1)}{(1-\xi_2)}-1, \,
410 * \eta_2 = \xi_2\f$
411 *
412 * \b Note: For the orthgonality of this expansion to be realised the
413 * 'q' ordering must run fastest in contrast to the Quad and Hex
414 * ordering where 'p' index runs fastest to be consistent with the
415 * quadrature ordering.
416 *
417 * In the triangular space the i (i.e. \f$\eta_1\f$ direction)
418 * ordering still runs fastest by convention.
419 */
421 Array<OneD, NekDouble> &outarray)
422{
423 StdTriExp::v_IProductWRTBase_SumFac(inarray, outarray);
424}
425
427 const Array<OneD, const NekDouble> &inarray,
428 Array<OneD, NekDouble> &outarray, bool multiplybyweights)
429{
430 int nquad0 = m_base[0]->GetNumPoints();
431 int nquad1 = m_base[1]->GetNumPoints();
432 int order0 = m_base[0]->GetNumModes();
433
434 if (multiplybyweights)
435 {
436 Array<OneD, NekDouble> tmp(nquad0 * nquad1 + nquad1 * order0);
437 Array<OneD, NekDouble> wsp(tmp + nquad0 * nquad1);
438
439 // multiply by integration constants
440 MultiplyByQuadratureMetric(inarray, tmp);
442 m_base[1]->GetBdata(), tmp, outarray, wsp);
443 }
444 else
445 {
446 Array<OneD, NekDouble> wsp(nquad1 * order0);
448 m_base[1]->GetBdata(), inarray, outarray,
449 wsp);
450 }
451}
452
454 const Array<OneD, const NekDouble> &base0,
455 const Array<OneD, const NekDouble> &base1,
456 const Array<OneD, const NekDouble> &inarray,
458 [[maybe_unused]] bool doCheckCollDir0,
459 [[maybe_unused]] bool doCheckCollDir1)
460{
461 int i;
462 int mode;
463 int nquad0 = m_base[0]->GetNumPoints();
464 int nquad1 = m_base[1]->GetNumPoints();
465 int nmodes0 = m_base[0]->GetNumModes();
466 int nmodes1 = m_base[1]->GetNumModes();
467
468 ASSERTL1(wsp.size() >= nquad1 * nmodes0,
469 "Workspace size is not sufficient");
470
471 Blas::Dgemm('T', 'N', nquad1, nmodes0, nquad0, 1.0, inarray.data(), nquad0,
472 base0.data(), nquad0, 0.0, wsp.data(), nquad1);
473
474 // Inner product with respect to 'b' direction
475 for (mode = i = 0; i < nmodes0; ++i)
476 {
477 Blas::Dgemv('T', nquad1, nmodes1 - i, 1.0, base1.data() + mode * nquad1,
478 nquad1, wsp.data() + i * nquad1, 1, 0.0,
479 outarray.data() + mode, 1);
480 mode += nmodes1 - i;
481 }
482
483 // fix for modified basis by splitting top vertex mode
485 {
486 outarray[1] += Blas::Ddot(nquad1, base1.data() + nquad1, 1,
487 wsp.data() + nquad1, 1);
488 }
489}
490
492 const int dir, const Array<OneD, const NekDouble> &inarray,
493 Array<OneD, NekDouble> &outarray)
494{
495 StdTriExp::v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
496}
497
499 const int dir, const Array<OneD, const NekDouble> &inarray,
500 Array<OneD, NekDouble> &outarray)
501{
502 int i;
503 int nquad0 = m_base[0]->GetNumPoints();
504 int nquad1 = m_base[1]->GetNumPoints();
505 int nqtot = nquad0 * nquad1;
506 int nmodes0 = m_base[0]->GetNumModes();
507 int wspsize = max(max(nqtot, m_ncoeffs), nquad1 * nmodes0);
508
509 Array<OneD, NekDouble> gfac0(2 * wspsize);
510 Array<OneD, NekDouble> tmp0(gfac0 + wspsize);
511
512 const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
513
514 // set up geometric factor: 2/(1-z1)
515 for (i = 0; i < nquad1; ++i)
516 {
517 gfac0[i] = 2.0 / (1 - z1[i]);
518 }
519
520 for (i = 0; i < nquad1; ++i)
521 {
522 Vmath::Smul(nquad0, gfac0[i], &inarray[0] + i * nquad0, 1,
523 &tmp0[0] + i * nquad0, 1);
524 }
525
526 MultiplyByQuadratureMetric(tmp0, tmp0);
527
528 switch (dir)
529 {
530 case 0:
531 {
532 IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
533 m_base[1]->GetBdata(), tmp0, outarray,
534 gfac0);
535 break;
536 }
537 case 1:
538 {
540 const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
541
542 for (i = 0; i < nquad0; ++i)
543 {
544 gfac0[i] = 0.5 * (1 + z0[i]);
545 }
546
547 for (i = 0; i < nquad1; ++i)
548 {
549 Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp0[0] + i * nquad0, 1,
550 &tmp0[0] + i * nquad0, 1);
551 }
552
553 IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
554 m_base[1]->GetBdata(), tmp0, tmp3,
555 gfac0);
556
557 MultiplyByQuadratureMetric(inarray, tmp0);
559 m_base[1]->GetDbdata(), tmp0, outarray,
560 gfac0);
561 Vmath::Vadd(m_ncoeffs, &tmp3[0], 1, &outarray[0], 1, &outarray[0],
562 1);
563 break;
564 }
565 default:
566 {
567 ASSERTL1(false, "input dir is out of range");
568 break;
569 }
570 }
571}
572
573//---------------------------------------
574// Evaluation functions
575//---------------------------------------
576
579{
580 NekDouble d1 = 1. - xi[1];
581 if (fabs(d1) < NekConstants::kNekZeroTol)
582 {
583 if (d1 >= 0.)
584 {
586 }
587 else
588 {
590 }
591 }
592 eta[0] = 2. * (1. + xi[0]) / d1 - 1.0;
593 eta[1] = xi[1];
594}
595
598{
599 xi[0] = (1.0 + eta[0]) * (1.0 - eta[1]) * 0.5 - 1.0;
600 xi[1] = eta[1];
601}
602
603void StdTriExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
604{
605 int i, m;
606 int nquad0 = m_base[0]->GetNumPoints();
607 int nquad1 = m_base[1]->GetNumPoints();
608 int order0 = m_base[0]->GetNumModes();
609 int order1 = m_base[1]->GetNumModes();
610 int mode0 = 0;
611 Array<OneD, const NekDouble> base0 = m_base[0]->GetBdata();
612 Array<OneD, const NekDouble> base1 = m_base[1]->GetBdata();
613
614 ASSERTL2(mode <= m_ncoeffs, "calling argument mode is larger than "
615 "total expansion order");
616
617 m = order1;
618 for (i = 0; i < order0; ++i, m += order1 - i)
619 {
620 if (m > mode)
621 {
622 mode0 = i;
623 break;
624 }
625 }
626
627 // deal with top vertex mode in modified basis
628 if (mode == 1 && m_base[0]->GetBasisType() == LibUtilities::eModified_A)
629 {
630 Vmath::Fill(nquad0 * nquad1, 1.0, outarray, 1);
631 }
632 else
633 {
634 for (i = 0; i < nquad1; ++i)
635 {
636 Vmath::Vcopy(nquad0, (NekDouble *)(base0.data() + mode0 * nquad0),
637 1, &outarray[0] + i * nquad0, 1);
638 }
639 }
640
641 for (i = 0; i < nquad0; ++i)
642 {
643 Vmath::Vmul(nquad1, (NekDouble *)(base1.data() + mode * nquad1), 1,
644 &outarray[0] + i, nquad0, &outarray[0] + i, nquad0);
645 }
646}
647
649 const Array<OneD, const NekDouble> &coords, int mode)
650{
652 LocCoordToLocCollapsed(coords, coll);
653
654 // From mode we need to determine mode0 and mode1 in the (p,q)
655 // direction. mode1 can be directly inferred from mode.
656 const int nm1 = m_base[1]->GetNumModes();
657 const double c = 1 + 2 * nm1;
658 const int mode0 = floor(0.5 * (c - sqrt(c * c - 8 * mode)));
659
660 if (mode == 1 && m_base[0]->GetBasisType() == LibUtilities::eModified_A)
661 {
662 // Account for collapsed vertex.
663 return StdExpansion::BaryEvaluateBasis<1>(coll[1], 1);
664 }
665 else
666 {
667 return StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
668 StdExpansion::BaryEvaluateBasis<1>(coll[1], mode);
669 }
670}
671
673 const Array<OneD, NekDouble> &coord,
674 const Array<OneD, const NekDouble> &inarray,
675 std::array<NekDouble, 3> &firstOrderDerivs)
676{
677 // Collapse coordinates
678 Array<OneD, NekDouble> coll(2, 0.0);
679 LocCoordToLocCollapsed(coord, coll);
680
681 // If near singularity do the old interpolation matrix method
682 if ((1 - coll[1]) < 1e-5)
683 {
684 int totPoints = GetTotPoints();
685 Array<OneD, NekDouble> EphysDeriv0(totPoints), EphysDeriv1(totPoints);
686 PhysDeriv(inarray, EphysDeriv0, EphysDeriv1);
687
689 I[0] = GetBase()[0]->GetI(coll);
690 I[1] = GetBase()[1]->GetI(coll + 1);
691
692 firstOrderDerivs[0] = PhysEvaluate(I, EphysDeriv0);
693 firstOrderDerivs[1] = PhysEvaluate(I, EphysDeriv1);
694 return PhysEvaluate(I, inarray);
695 }
696
697 // set up geometric factor: 2.0/(1.0-z1)
698 NekDouble fac0 = 2 / (1 - coll[1]);
699
700 NekDouble val = BaryTensorDeriv(coll, inarray, firstOrderDerivs);
701
702 // Copy d0 into temp for d1
703 NekDouble temp;
704 temp = firstOrderDerivs[0];
705
706 // Multiply by geometric factor
707 firstOrderDerivs[0] = firstOrderDerivs[0] * fac0;
708
709 // set up geometric factor: (1+z0)/(1-z1)
710 NekDouble fac1 = fac0 * (coll[0] + 1) / 2;
711
712 // Multiply out_d0 by geometric factor and add to out_d1
713 firstOrderDerivs[1] += fac1 * temp;
714
715 return val;
716}
717
719{
720 return 3;
721}
722
724{
725 return 3;
726}
727
729{
731}
732
734{
736 "BasisType is not a boundary interior form");
738 "BasisType is not a boundary interior form");
739
740 return 3 + (GetBasisNumModes(0) - 2) + 2 * (GetBasisNumModes(1) - 2);
741}
742
744{
746 "BasisType is not a boundary interior form");
748 "BasisType is not a boundary interior form");
749
750 return GetBasisNumModes(0) + 2 * GetBasisNumModes(1);
751}
752
753int StdTriExp::v_GetTraceNcoeffs(const int i) const
754{
755 ASSERTL2(i >= 0 && i <= 2, "edge id is out of range");
756
757 if (i == 0)
758 {
759 return GetBasisNumModes(0);
760 }
761 else
762 {
763 return GetBasisNumModes(1);
764 }
765}
766
768{
769 ASSERTL2(i >= 0 && i <= 2, "edge id is out of range");
770
771 if (i == 0)
772 {
773 return GetBasisNumModes(0) - 2;
774 }
775 else
776 {
777 return GetBasisNumModes(1) - 2;
778 }
779}
780
781int StdTriExp::v_GetTraceNumPoints(const int i) const
782{
783 ASSERTL2((i >= 0) && (i <= 2), "edge id is out of range");
784
785 if (i == 0)
786 {
787 return GetNumPoints(0);
788 }
789 else
790 {
791 return GetNumPoints(1);
792 }
793}
794
796 const std::vector<unsigned int> &nummodes, int &modes_offset)
797{
799 nummodes[modes_offset], nummodes[modes_offset + 1]);
800 modes_offset += 2;
801
802 return nmodes;
803}
804
806 Array<OneD, NekDouble> &coords_1,
807 [[maybe_unused]] Array<OneD, NekDouble> &coords_2)
808{
809 Array<OneD, const NekDouble> z0 = m_base[0]->GetZ();
810 Array<OneD, const NekDouble> z1 = m_base[1]->GetZ();
811 int nq0 = GetNumPoints(0);
812 int nq1 = GetNumPoints(1);
813 int i, j;
814
815 for (i = 0; i < nq1; ++i)
816 {
817 for (j = 0; j < nq0; ++j)
818 {
819 coords_0[i * nq0 + j] = (1 + z0[j]) * (1 - z1[i]) / 2.0 - 1.0;
820 }
821 Vmath::Fill(nq0, z1[i], &coords_1[0] + i * nq0, 1);
822 }
823}
824
826{
827 return m_base[0]->GetBasisType() == LibUtilities::eModified_A &&
828 m_base[1]->GetBasisType() == LibUtilities::eModified_B;
829}
830
832 const int i, [[maybe_unused]] const int j,
833 [[maybe_unused]] bool UseGLL) const
834{
835 ASSERTL2(i >= 0 && i <= 2, "edge id is out of range");
836
837 // Get basiskey (0 or 1) according to edge id i
838 int dir = (i != 0);
839
840 switch (m_base[dir]->GetBasisType())
841 {
844 {
845 switch (m_base[dir]->GetPointsType())
846 {
848 {
849 return m_base[dir]->GetBasisKey();
850 }
851 break;
852 default:
853 {
855 "Unexpected points distribution " +
857 [m_base[dir]->GetPointsType()] +
858 " in StdTriExp::v_GetTraceBasisKey");
859 }
860 }
861 }
862 break;
864 {
865 switch (m_base[dir]->GetPointsType())
866 {
867 case LibUtilities::eGaussRadauMAlpha1Beta0:
868 {
870 m_base[dir]
871 ->GetBasisKey()
872 .GetPointsKey()
873 .GetNumPoints() +
874 1,
877 m_base[dir]->GetNumModes(),
878 pkey);
879 }
880 break;
881 // Currently this does not increase the points by
882 // 1 since when using this quadrature we are
883 // presuming it is already been increased by one
884 // when comopared to
885 // GaussRadauMAlpha1Beta0. Currently used in the
886 // GJP option
887 //
888 // Note have put down it back to numpoints +1 to
889 // test for use on tri faces and GJP.
891 {
893 m_base[dir]
894 ->GetBasisKey()
895 .GetPointsKey()
896 .GetNumPoints() +
897 1,
900 m_base[dir]->GetNumModes(),
901 pkey);
902 }
903 break;
905 {
907 m_base[dir]
908 ->GetBasisKey()
909 .GetPointsKey()
910 .GetNumPoints() +
911 1,
914 m_base[dir]->GetNumModes(),
915 pkey);
916 }
917 break;
918 default:
919 {
921 "Unexpected points distribution " +
923 [m_base[dir]->GetPointsType()] +
924 " in StdTriExp::v_GetTraceBasisKey");
925 }
926 }
927 }
928 break;
930 {
931 switch (m_base[dir]->GetPointsType())
932 {
933 case LibUtilities::eGaussRadauMAlpha1Beta0:
934 {
936 m_base[dir]
937 ->GetBasisKey()
938 .GetPointsKey()
939 .GetNumPoints() +
940 1,
943 m_base[dir]->GetNumModes(),
944 pkey);
945 }
946 break;
947 default:
948 {
950 "Unexpected points distribution " +
952 [m_base[dir]->GetPointsType()] +
953 " in StdTriExp::v_GetTraceBasisKey");
954 }
955 }
956 }
957 break;
958 default:
959 {
961 "Information not available to set edge key");
962 }
963 }
965}
966
967//--------------------------
968// Mappings
969//--------------------------
970
971int StdTriExp::v_GetVertexMap(const int localVertexId, bool useCoeffPacking)
972{
975 "Mapping not defined for this type of basis");
976
977 int localDOF = 0;
978 if (useCoeffPacking == true)
979 {
980 switch (localVertexId)
981 {
982 case 0:
983 {
984 localDOF = 0;
985 break;
986 }
987 case 1:
988 {
989 localDOF = 1;
990 break;
991 }
992 case 2:
993 {
994 localDOF = m_base[1]->GetNumModes();
995 break;
996 }
997 default:
998 {
999 ASSERTL0(false, "eid must be between 0 and 2");
1000 break;
1001 }
1002 }
1003 }
1004 else // follow book format for vertex indexing.
1005 {
1006 switch (localVertexId)
1007 {
1008 case 0:
1009 {
1010 localDOF = 0;
1011 break;
1012 }
1013 case 1:
1014 {
1015 localDOF = m_base[1]->GetNumModes();
1016 break;
1017 }
1018 case 2:
1019 {
1020 localDOF = 1;
1021 break;
1022 }
1023 default:
1024 {
1025 ASSERTL0(false, "eid must be between 0 and 2");
1026 break;
1027 }
1028 }
1029 }
1030
1031 return localDOF;
1032}
1033
1035{
1038 "Expansion not of a proper type");
1039
1040 int i, j;
1041 int cnt = 0;
1042 int nummodes0, nummodes1;
1043 int startvalue;
1044 if (outarray.size() != GetNcoeffs() - NumBndryCoeffs())
1045 {
1047 }
1048
1049 nummodes0 = m_base[0]->GetNumModes();
1050 nummodes1 = m_base[1]->GetNumModes();
1051
1052 startvalue = 2 * nummodes1;
1053
1054 for (i = 0; i < nummodes0 - 2; i++)
1055 {
1056 for (j = 0; j < nummodes1 - 3 - i; j++)
1057 {
1058 outarray[cnt++] = startvalue + j;
1059 }
1060 startvalue += nummodes1 - 2 - i;
1061 }
1062}
1063
1065{
1068 "Expansion not of expected type");
1069 int i;
1070 int cnt;
1071 int nummodes0, nummodes1;
1072 int value;
1073
1074 if (outarray.size() != NumBndryCoeffs())
1075 {
1077 }
1078
1079 nummodes0 = m_base[0]->GetNumModes();
1080 nummodes1 = m_base[1]->GetNumModes();
1081
1082 value = 2 * nummodes1 - 1;
1083 for (i = 0; i < value; i++)
1084 {
1085 outarray[i] = i;
1086 }
1087 cnt = value;
1088
1089 for (i = 0; i < nummodes0 - 2; i++)
1090 {
1091 outarray[cnt++] = value;
1092 value += nummodes1 - 2 - i;
1093 }
1094}
1095
1096void StdTriExp::v_GetTraceCoeffMap(const unsigned int eid,
1097 Array<OneD, unsigned int> &maparray)
1098{
1099
1102 "Mapping not defined for this type of basis");
1103
1104 ASSERTL1(eid < 3, "eid must be between 0 and 2");
1105
1106 int i;
1107 int order0 = m_base[0]->GetNumModes();
1108 int order1 = m_base[1]->GetNumModes();
1109 int numModes = (eid == 0) ? order0 : order1;
1110
1111 if (maparray.size() != numModes)
1112 {
1113 maparray = Array<OneD, unsigned int>(numModes);
1114 }
1115
1116 switch (eid)
1117 {
1118 case 0:
1119 {
1120 int cnt = 0;
1121 for (i = 0; i < numModes; cnt += order1 - i, ++i)
1122 {
1123 maparray[i] = cnt;
1124 }
1125 break;
1126 }
1127 case 1:
1128 {
1129 maparray[0] = order1;
1130 maparray[1] = 1;
1131 for (i = 2; i < numModes; i++)
1132 {
1133 maparray[i] = order1 - 1 + i;
1134 }
1135 break;
1136 }
1137 case 2:
1138 {
1139 for (i = 0; i < numModes; i++)
1140 {
1141 maparray[i] = i;
1142 }
1143 break;
1144 }
1145 default:
1146 ASSERTL0(false, "eid must be between 0 and 2");
1147 break;
1148 }
1149}
1150
1152 const int eid, Array<OneD, unsigned int> &maparray,
1153 Array<OneD, int> &signarray, const Orientation edgeOrient)
1154{
1157 "Mapping not defined for this type of basis");
1158 int i;
1159 const int nummodes1 = m_base[1]->GetNumModes();
1160 const int nEdgeIntCoeffs = GetTraceNcoeffs(eid) - 2;
1161
1162 if (maparray.size() != nEdgeIntCoeffs)
1163 {
1164 maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1165 }
1166
1167 if (signarray.size() != nEdgeIntCoeffs)
1168 {
1169 signarray = Array<OneD, int>(nEdgeIntCoeffs, 1);
1170 }
1171 else
1172 {
1173 fill(signarray.data(), signarray.data() + nEdgeIntCoeffs, 1);
1174 }
1175
1176 switch (eid)
1177 {
1178 case 0:
1179 {
1180 int cnt = 2 * nummodes1 - 1;
1181 for (i = 0; i < nEdgeIntCoeffs; cnt += nummodes1 - 2 - i, ++i)
1182 {
1183 maparray[i] = cnt;
1184 }
1185 break;
1186 }
1187 case 1:
1188 {
1189 for (i = 0; i < nEdgeIntCoeffs; i++)
1190 {
1191 maparray[i] = nummodes1 + 1 + i;
1192 }
1193 break;
1194 }
1195 case 2:
1196 {
1197 for (i = 0; i < nEdgeIntCoeffs; i++)
1198 {
1199 maparray[i] = 2 + i;
1200 }
1201 break;
1202 }
1203 default:
1204 {
1205 ASSERTL0(false, "eid must be between 0 and 2");
1206 break;
1207 }
1208 }
1209
1210 if (edgeOrient == eBackwards)
1211 {
1212 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1213 {
1214 signarray[i] = -1;
1215 }
1216 }
1217}
1218
1219//---------------------------------------
1220// Wrapper functions
1221//---------------------------------------
1222
1224{
1225
1226 MatrixType mtype = mkey.GetMatrixType();
1227
1228 DNekMatSharedPtr Mat;
1229
1230 switch (mtype)
1231 {
1233 {
1234 int nq0, nq1, nq;
1235
1236 nq0 = m_base[0]->GetNumPoints();
1237 nq1 = m_base[1]->GetNumPoints();
1238
1239 // take definition from key
1241 {
1242 nq = (int)mkey.GetConstFactor(eFactorConst);
1243 }
1244 else
1245 {
1246 nq = max(nq0, nq1);
1247 }
1248
1251 Array<OneD, NekDouble> coll(2);
1253 Array<OneD, NekDouble> tmp(nq0);
1254
1255 Mat = MemoryManager<DNekMat>::AllocateSharedPtr(neq, nq0 * nq1);
1256 int cnt = 0;
1257
1258 for (int i = 0; i < nq; ++i)
1259 {
1260 for (int j = 0; j < nq - i; ++j, ++cnt)
1261 {
1262 coords[cnt] = Array<OneD, NekDouble>(2);
1263 coords[cnt][0] = -1.0 + 2 * j / (NekDouble)(nq - 1);
1264 coords[cnt][1] = -1.0 + 2 * i / (NekDouble)(nq - 1);
1265 }
1266 }
1267
1268 for (int i = 0; i < neq; ++i)
1269 {
1270 LocCoordToLocCollapsed(coords[i], coll);
1271
1272 I[0] = m_base[0]->GetI(coll);
1273 I[1] = m_base[1]->GetI(coll + 1);
1274
1275 // interpolate first coordinate direction
1276 for (int j = 0; j < nq1; ++j)
1277 {
1278 NekDouble fac = (I[1]->GetPtr())[j];
1279 Vmath::Smul(nq0, fac, I[0]->GetPtr(), 1, tmp, 1);
1280
1281 Vmath::Vcopy(nq0, &tmp[0], 1,
1282 Mat->GetRawPtr() + j * nq0 * neq + i, neq);
1283 }
1284 }
1285 break;
1286 }
1287 default:
1288 {
1290 break;
1291 }
1292 }
1293
1294 return Mat;
1295}
1296
1298{
1299 return v_GenMatrix(mkey);
1300}
1301
1302//---------------------------------------
1303// Operator evaluation functions
1304//---------------------------------------
1305
1307 Array<OneD, NekDouble> &outarray,
1308 const StdMatrixKey &mkey)
1309{
1310 StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
1311}
1312
1314 Array<OneD, NekDouble> &outarray,
1315 const StdMatrixKey &mkey)
1316{
1317 StdTriExp::v_LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
1318}
1319
1320void StdTriExp::v_LaplacianMatrixOp(const int k1, const int k2,
1321 const Array<OneD, const NekDouble> &inarray,
1322 Array<OneD, NekDouble> &outarray,
1323 const StdMatrixKey &mkey)
1324{
1325 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1326}
1327
1329 const Array<OneD, const NekDouble> &inarray,
1330 Array<OneD, NekDouble> &outarray,
1331 const StdMatrixKey &mkey)
1332{
1333 StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
1334}
1335
1337 Array<OneD, NekDouble> &outarray,
1338 const StdMatrixKey &mkey)
1339{
1340 StdTriExp::v_HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
1341}
1342
1344 const StdMatrixKey &mkey)
1345{
1346 int qa = m_base[0]->GetNumPoints();
1347 int qb = m_base[1]->GetNumPoints();
1348 int nmodes_a = m_base[0]->GetNumModes();
1349 int nmodes_b = m_base[1]->GetNumModes();
1350
1351 // Declare orthogonal basis.
1354
1357 StdTriExp OrthoExp(Ba, Bb);
1358
1359 Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1360
1361 // project onto physical space.
1362 OrthoExp.FwdTrans(array, orthocoeffs);
1363
1364 if (mkey.ConstFactorExists(
1365 eFactorSVVPowerKerDiffCoeff)) // Rodrigo's power kern
1366 {
1368 NekDouble SvvDiffCoeff =
1371
1372 int cnt = 0;
1373 for (int j = 0; j < nmodes_a; ++j)
1374 {
1375 for (int k = 0; k < nmodes_b - j; ++k, ++cnt)
1376 {
1377 NekDouble fac = std::max(
1378 pow((1.0 * j) / (nmodes_a - 1), cutoff * nmodes_a),
1379 pow((1.0 * k) / (nmodes_b - 1), cutoff * nmodes_b));
1380
1381 orthocoeffs[cnt] *= (SvvDiffCoeff * fac);
1382 }
1383 }
1384 }
1385 else if (mkey.ConstFactorExists(
1386 eFactorSVVDGKerDiffCoeff)) // Rodrigo/mansoor's DG kernel
1387 {
1390 int max_ab = max(nmodes_a - kSVVDGFiltermodesmin,
1391 nmodes_b - kSVVDGFiltermodesmin);
1392 max_ab = max(max_ab, 0);
1393 max_ab = min(max_ab, kSVVDGFiltermodesmax - kSVVDGFiltermodesmin);
1394
1395 int cnt = 0;
1396 for (int j = 0; j < nmodes_a; ++j)
1397 {
1398 for (int k = 0; k < nmodes_b - j; ++k, ++cnt)
1399 {
1400 int maxjk = max(j, k);
1401 maxjk = min(maxjk, kSVVDGFiltermodesmax - 1);
1402
1403 orthocoeffs[cnt] *= SvvDiffCoeff * kSVVDGFilter[max_ab][maxjk];
1404 }
1405 }
1406 }
1407 else
1408 {
1409 NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
1410
1411 int cutoff = (int)(mkey.GetConstFactor(eFactorSVVCutoffRatio) *
1412 min(nmodes_a, nmodes_b));
1413
1414 NekDouble epsilon = 1.0;
1415 int nmodes = min(nmodes_a, nmodes_b);
1416
1417 int cnt = 0;
1418
1419 // apply SVV filter (JEL)
1420 for (int j = 0; j < nmodes_a; ++j)
1421 {
1422 for (int k = 0; k < nmodes_b - j; ++k)
1423 {
1424 if (j + k >= cutoff)
1425 {
1426 orthocoeffs[cnt] *=
1427 (SvvDiffCoeff *
1428 exp(-(j + k - nmodes) * (j + k - nmodes) /
1429 ((NekDouble)((j + k - cutoff + epsilon) *
1430 (j + k - cutoff + epsilon)))));
1431 }
1432 else
1433 {
1434 orthocoeffs[cnt] *= 0.0;
1435 }
1436 cnt++;
1437 }
1438 }
1439 }
1440
1441 // backward transform to physical space
1442 OrthoExp.BwdTrans(orthocoeffs, array);
1443}
1444
1446 const Array<OneD, const NekDouble> &inarray,
1447 Array<OneD, NekDouble> &outarray)
1448{
1449 int n_coeffs = inarray.size();
1450 int nquad0 = m_base[0]->GetNumPoints();
1451 int nquad1 = m_base[1]->GetNumPoints();
1452 Array<OneD, NekDouble> coeff(n_coeffs);
1453 Array<OneD, NekDouble> coeff_tmp(n_coeffs, 0.0);
1456 int nqtot = nquad0 * nquad1;
1457 Array<OneD, NekDouble> phys_tmp(nqtot, 0.0);
1458
1459 int nmodes0 = m_base[0]->GetNumModes();
1460 int nmodes1 = m_base[1]->GetNumModes();
1461 int numMin2 = nmodes0;
1462 int i;
1463
1464 const LibUtilities::PointsKey Pkey0(nmodes0,
1466 const LibUtilities::PointsKey Pkey1(nmodes1,
1468
1469 LibUtilities::BasisKey b0(m_base[0]->GetBasisType(), nmodes0, Pkey0);
1470 LibUtilities::BasisKey b1(m_base[1]->GetBasisType(), nmodes1, Pkey1);
1471
1472 LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A, nmodes0, Pkey0);
1473 LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_B, nmodes1, Pkey1);
1474
1475 StdRegions::StdTriExpSharedPtr m_OrthoTriExp;
1477
1480 bortho0, bortho1);
1481
1482 m_TriExp->BwdTrans(inarray, phys_tmp);
1483 m_OrthoTriExp->FwdTrans(phys_tmp, coeff);
1484
1485 for (i = 0; i < n_coeffs; i++)
1486 {
1487 if (i == numMin)
1488 {
1489 coeff[i] = 0.0;
1490 numMin += numMin2 - 1;
1491 numMin2 -= 1.0;
1492 }
1493 }
1494
1495 m_OrthoTriExp->BwdTrans(coeff, phys_tmp);
1496 m_TriExp->FwdTrans(phys_tmp, outarray);
1497}
1498
1499//---------------------------------------
1500// Private helper functions
1501//---------------------------------------
1502
1504 const Array<OneD, const NekDouble> &inarray,
1505 Array<OneD, NekDouble> &outarray)
1506{
1507 int i;
1508 int nquad0 = m_base[0]->GetNumPoints();
1509 int nquad1 = m_base[1]->GetNumPoints();
1510
1511 const Array<OneD, const NekDouble> &w0 = m_base[0]->GetW();
1512 const Array<OneD, const NekDouble> &w1 = m_base[1]->GetW();
1513 const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
1514
1515 // multiply by integration constants
1516 for (i = 0; i < nquad1; ++i)
1517 {
1518 Vmath::Vmul(nquad0, inarray.data() + i * nquad0, 1, w0.data(), 1,
1519 outarray.data() + i * nquad0, 1);
1520 }
1521
1522 switch (m_base[1]->GetPointsType())
1523 {
1524 // (1,0) Jacobi Inner product
1525 case LibUtilities::eGaussRadauMAlpha1Beta0:
1526 for (i = 0; i < nquad1; ++i)
1527 {
1528 Blas::Dscal(nquad0, 0.5 * w1[i], outarray.data() + i * nquad0,
1529 1);
1530 }
1531 break;
1532 // Legendre inner product
1533 default:
1534 for (i = 0; i < nquad1; ++i)
1535 {
1536 Blas::Dscal(nquad0, 0.5 * (1 - z1[i]) * w1[i],
1537 outarray.data() + i * nquad0, 1);
1538 }
1539 break;
1540 }
1541}
1542
1544 Array<OneD, int> &conn, [[maybe_unused]] bool standard)
1545{
1546 int np1 = m_base[0]->GetNumPoints();
1547 int np2 = m_base[1]->GetNumPoints();
1548 int np = max(np1, np2);
1549
1550 conn = Array<OneD, int>(3 * (np - 1) * (np - 1));
1551
1552 int row = 0;
1553 int rowp1 = 0;
1554 int cnt = 0;
1555 for (int i = 0; i < np - 1; ++i)
1556 {
1557 rowp1 += np - i;
1558 for (int j = 0; j < np - i - 2; ++j)
1559 {
1560 conn[cnt++] = row + j;
1561 conn[cnt++] = row + j + 1;
1562 conn[cnt++] = rowp1 + j;
1563
1564 conn[cnt++] = rowp1 + j + 1;
1565 conn[cnt++] = rowp1 + j;
1566 conn[cnt++] = row + j + 1;
1567 }
1568
1569 conn[cnt++] = row + np - i - 2;
1570 conn[cnt++] = row + np - i - 1;
1571 conn[cnt++] = rowp1 + np - i - 2;
1572
1573 row += np - i;
1574 }
1575}
1576} // namespace Nektar::StdRegions
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
#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
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:74
Defines a specification for a set of points.
Definition: Points.h:50
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
void 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
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:124
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:134
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
const Array< OneD, const LibUtilities::BasisSharedPtr > & GetBase() const
This function gets the shared point to basis.
Definition: StdExpansion.h:100
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
NekDouble PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals)
This function evaluates the expansion at a single (arbitrary) point of the domain.
Definition: StdExpansion.h:928
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 BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space.
Definition: StdExpansion.h:433
int GetTraceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th trace.
Definition: StdExpansion.h:261
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:205
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space.
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
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Definition: StdExpansion.h:858
Array< OneD, LibUtilities::BasisSharedPtr > m_base
void MassMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:83
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:124
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:133
void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
Definition: StdTriExp.cpp:1328
LibUtilities::ShapeType v_DetShapeType() const final
Definition: StdTriExp.cpp:728
int v_GetTraceNumPoints(const int i) const override
Definition: StdTriExp.cpp:781
void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Backward tranform for triangular elements.
Definition: StdTriExp.cpp:217
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdTriExp.cpp:491
void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi) override
Definition: StdTriExp.cpp:596
void v_GetCoords(Array< OneD, NekDouble > &coords_x, Array< OneD, NekDouble > &coords_y, Array< OneD, NekDouble > &coords_z) override
Definition: StdTriExp.cpp:805
int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset) override
Definition: StdTriExp.cpp:795
const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int j, bool UseGLL=false) const override
Definition: StdTriExp.cpp:831
void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true) override
Definition: StdTriExp.cpp:1543
void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
Definition: StdTriExp.cpp:1306
void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray) override
Definition: StdTriExp.cpp:1064
int v_GetNtraces() const final
Definition: StdTriExp.cpp:723
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: StdTriExp.cpp:193
void v_GetTraceInteriorToElementMap(const int eid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation edgeOrient=eForwards) override
Definition: StdTriExp.cpp:1151
int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false) override
Definition: StdTriExp.cpp:971
DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey) override
Definition: StdTriExp.cpp:1297
void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdTriExp.cpp:223
DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey) override
Definition: StdTriExp.cpp:1223
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[p]*base1[pq] and put into ou...
Definition: StdTriExp.cpp:420
int v_GetTraceNcoeffs(const int i) const override
Definition: StdTriExp.cpp:753
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: StdTriExp.cpp:106
void v_GetTraceCoeffMap(const unsigned int traceid, Array< OneD, unsigned int > &maparray) override
Definition: StdTriExp.cpp:1096
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: StdTriExp.cpp:233
bool v_IsBoundaryInteriorExpansion() const override
Definition: StdTriExp.cpp:825
void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray) override
Definition: StdTriExp.cpp:603
void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdTriExp.cpp:1445
void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdTriExp.cpp:289
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdTriExp.cpp:498
int v_GetTraceIntNcoeffs(const int i) const override
Definition: StdTriExp.cpp:767
int v_NumDGBndryCoeffs() const override
Definition: StdTriExp.cpp:743
void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta) override
Definition: StdTriExp.cpp:577
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) final
Definition: StdTriExp.cpp:648
void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
Definition: StdTriExp.cpp:1336
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray) override
Integrates the specified function over the domain.
Definition: StdTriExp.cpp:58
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
Definition: StdTriExp.cpp:672
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: StdTriExp.cpp:273
void v_GetInteriorMap(Array< OneD, unsigned int > &outarray) override
Definition: StdTriExp.cpp:1034
int v_GetNverts() const final
Definition: StdTriExp.cpp:718
void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
Definition: StdTriExp.cpp:1313
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey) override
Definition: StdTriExp.cpp:1343
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdTriExp.cpp:1503
int v_NumBndryCoeffs() const override
Definition: StdTriExp.cpp:733
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
Definition: StdTriExp.cpp:426
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: StdTriExp.cpp:453
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = alpha A x plus beta y where A[m x n].
Definition: Blas.hpp:211
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:149
static double Ddot(const int &n, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: output = .
Definition: Blas.hpp:163
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition: Blas.hpp:383
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
Definition: Blas.hpp:135
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:109
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
const std::string kPointsTypeStr[]
Definition: Foundations.hpp:52
@ eGaussRadauMLegendre
1D Gauss-Radau-Legendre quadrature points, pinned at x=-1
Definition: PointsType.h:47
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
@ ePolyEvenlySpaced
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:73
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:42
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:56
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:44
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
static const NekDouble kNekZeroTol
const int kSVVDGFiltermodesmin
Definition: StdRegions.hpp:500
const int kSVVDGFiltermodesmax
Definition: StdRegions.hpp:501
const NekDouble kSVVDGFilter[9][11]
Definition: StdRegions.hpp:503
std::shared_ptr< StdTriExp > StdTriExpSharedPtr
Definition: StdTriExp.h:218
std::shared_ptr< StdSegExp > StdSegExpSharedPtr
Definition: StdSegExp.h:217
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 Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.hpp:366
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/x.
Definition: Vmath.hpp:154
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.hpp:194
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.
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285