Nektar++
TriExp.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: TriExp.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// The above copyright notice and this permission notice shall be included
20// in all copies or substantial portions of the Software.
21//
22// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28// DEALINGS IN THE SOFTWARE.
29//
30// Description: Expasion for triangular elements.
31//
32///////////////////////////////////////////////////////////////////////////////
33
37#include <LocalRegions/SegExp.h>
38#include <LocalRegions/TriExp.h>
40
41using namespace std;
42
44{
46 const LibUtilities::BasisKey &Bb,
48 : StdExpansion(LibUtilities::StdTriData::getNumberOfCoefficients(
49 Ba.GetNumModes(), (Bb.GetNumModes())),
50 2, Ba, Bb),
51 StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(
52 Ba.GetNumModes(), (Bb.GetNumModes())),
53 Ba, Bb),
54 StdTriExp(Ba, Bb), Expansion(geom), Expansion2D(geom),
55 m_matrixManager(
56 std::bind(&Expansion2D::CreateMatrix, this, std::placeholders::_1),
57 std::string("TriExpMatrix")),
58 m_staticCondMatrixManager(std::bind(&Expansion::CreateStaticCondMatrix,
59 this, std::placeholders::_1),
60 std::string("TriExpStaticCondMatrix"))
61{
62}
63
65 : StdExpansion(T), StdExpansion2D(T), StdTriExp(T), Expansion(T),
66 Expansion2D(T), m_matrixManager(T.m_matrixManager),
67 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
68{
69}
70
72{
73 int nquad0 = m_base[0]->GetNumPoints();
74 int nquad1 = m_base[1]->GetNumPoints();
76 NekDouble ival;
77 Array<OneD, NekDouble> tmp(nquad0 * nquad1);
78
79 // multiply inarray with Jacobian
80 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
81 {
82 Vmath::Vmul(nquad0 * nquad1, jac, 1, inarray, 1, tmp, 1);
83 }
84 else
85 {
86 Vmath::Smul(nquad0 * nquad1, jac[0], inarray, 1, tmp, 1);
87 }
88
89 // call StdQuadExp version;
90 ival = StdTriExp::v_Integral(tmp);
91 return ival;
92}
93
98{
99 int nquad0 = m_base[0]->GetNumPoints();
100 int nquad1 = m_base[1]->GetNumPoints();
101 int nqtot = nquad0 * nquad1;
103 m_metricinfo->GetDerivFactors(GetPointsKeys());
104
105 Array<OneD, NekDouble> diff0(2 * nqtot);
106 Array<OneD, NekDouble> diff1(diff0 + nqtot);
107
108 StdTriExp::v_PhysDeriv(inarray, diff0, diff1);
109
110 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
111 {
112 if (out_d0.size())
113 {
114 Vmath::Vmul(nqtot, df[0], 1, diff0, 1, out_d0, 1);
115 Vmath::Vvtvp(nqtot, df[1], 1, diff1, 1, out_d0, 1, out_d0, 1);
116 }
117
118 if (out_d1.size())
119 {
120 Vmath::Vmul(nqtot, df[2], 1, diff0, 1, out_d1, 1);
121 Vmath::Vvtvp(nqtot, df[3], 1, diff1, 1, out_d1, 1, out_d1, 1);
122 }
123
124 if (out_d2.size())
125 {
126 Vmath::Vmul(nqtot, df[4], 1, diff0, 1, out_d2, 1);
127 Vmath::Vvtvp(nqtot, df[5], 1, diff1, 1, out_d2, 1, out_d2, 1);
128 }
129 }
130 else // regular geometry
131 {
132 if (out_d0.size())
133 {
134 Vmath::Smul(nqtot, df[0][0], diff0, 1, out_d0, 1);
135 Blas::Daxpy(nqtot, df[1][0], diff1, 1, out_d0, 1);
136 }
137
138 if (out_d1.size())
139 {
140 Vmath::Smul(nqtot, df[2][0], diff0, 1, out_d1, 1);
141 Blas::Daxpy(nqtot, df[3][0], diff1, 1, out_d1, 1);
142 }
143
144 if (out_d2.size())
145 {
146 Vmath::Smul(nqtot, df[4][0], diff0, 1, out_d2, 1);
147 Blas::Daxpy(nqtot, df[5][0], diff1, 1, out_d2, 1);
148 }
149 }
150}
151
152void TriExp::v_PhysDeriv(const int dir,
153 const Array<OneD, const NekDouble> &inarray,
154 Array<OneD, NekDouble> &outarray)
155{
156 switch (dir)
157 {
158 case 0:
159 {
160 PhysDeriv(inarray, outarray, NullNekDouble1DArray,
162 }
163 break;
164 case 1:
165 {
166 PhysDeriv(inarray, NullNekDouble1DArray, outarray,
168 }
169 break;
170 case 2:
171 {
173 outarray);
174 }
175 break;
176 default:
177 {
178 ASSERTL1(false, "input dir is out of range");
179 }
180 break;
181 }
182}
183
185 const Array<OneD, const NekDouble> &inarray,
187{
188 if (!out.size())
189 {
190 return;
191 }
192
193 int nquad0 = m_base[0]->GetNumPoints();
194 int nquad1 = m_base[1]->GetNumPoints();
195 int nqtot = nquad0 * nquad1;
196
198 m_metricinfo->GetDerivFactors(GetPointsKeys());
199
200 Array<OneD, NekDouble> diff0(2 * nqtot);
201 Array<OneD, NekDouble> diff1(diff0 + nqtot);
202
203 // diff0 = du/d_xi, diff1 = du/d_eta
204 StdTriExp::v_PhysDeriv(inarray, diff0, diff1);
205
206 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
207 {
209
210 // D^v_xi = v_x*d_xi/dx + v_y*d_xi/dy + v_z*d_xi/dz
211 // D^v_eta = v_x*d_eta/dx + v_y*d_eta/dy + v_z*d_eta/dz
212 for (int i = 0; i < 2; ++i)
213 {
214 tangmat[i] = Array<OneD, NekDouble>(nqtot, 0.0);
215 for (int k = 0; k < (m_geom->GetCoordim()); ++k)
216 {
217 Vmath::Vvtvp(nqtot, &df[2 * k + i][0], 1, &direction[k * nqtot],
218 1, &tangmat[i][0], 1, &tangmat[i][0], 1);
219 }
220 }
221
222 /// D_v = D^v_xi * du/d_xi + D^v_eta * du/d_eta
223 Vmath::Vmul(nqtot, &tangmat[0][0], 1, &diff0[0], 1, &out[0], 1);
224 Vmath::Vvtvp(nqtot, &tangmat[1][0], 1, &diff1[0], 1, &out[0], 1,
225 &out[0], 1);
226 }
227 else
228 {
230
231 for (int i = 0; i < 2; ++i)
232 {
233 tangmat[i] = Array<OneD, NekDouble>(nqtot, 0.0);
234 for (int k = 0; k < (m_geom->GetCoordim()); ++k)
235 {
236 Vmath::Svtvp(nqtot, df[2 * k + i][0], &direction[k * nqtot], 1,
237 &tangmat[i][0], 1, &tangmat[i][0], 1);
238 }
239 }
240
241 /// D_v = D^v_xi * du/d_xi + D^v_eta * du/d_eta
242 Vmath::Vmul(nqtot, &tangmat[0][0], 1, &diff0[0], 1, &out[0], 1);
243
244 Vmath::Vvtvp(nqtot, &tangmat[1][0], 1, &diff1[0], 1, &out[0], 1,
245 &out[0], 1);
246 }
247}
248
250 Array<OneD, NekDouble> &outarray)
251{
252 IProductWRTBase(inarray, outarray);
253
254 // get Mass matrix inverse
256 DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
257
258 // copy inarray in case inarray == outarray
259 NekVector<NekDouble> in(m_ncoeffs, outarray, eCopy);
261
262 out = (*matsys) * in;
263}
264
266 const Array<OneD, const NekDouble> &inarray,
267 Array<OneD, NekDouble> &outarray)
268{
269 int i, j;
270 int npoints[2] = {m_base[0]->GetNumPoints(), m_base[1]->GetNumPoints()};
271 int nmodes[2] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes()};
272
273 fill(outarray.data(), outarray.data() + m_ncoeffs, 0.0);
274
275 if (nmodes[0] == 1 && nmodes[1] == 1)
276 {
277 outarray[0] = inarray[0];
278 return;
279 }
280
281 Array<OneD, NekDouble> physEdge[3];
282 Array<OneD, NekDouble> coeffEdge[3];
283 for (i = 0; i < 3; i++)
284 {
285 // define physEdge and add 1 so can interpolate grl10 points if
286 // necessary
287 physEdge[i] = Array<OneD, NekDouble>(max(npoints[i != 0], npoints[0]));
288 coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i != 0]);
289 }
290
291 for (i = 0; i < npoints[0]; i++)
292 {
293 physEdge[0][i] = inarray[i];
294 }
295
296 // extract data in cartesian directions
297 for (i = 0; i < npoints[1]; i++)
298 {
299 physEdge[1][i] = inarray[npoints[0] - 1 + i * npoints[0]];
300 physEdge[2][i] = inarray[i * npoints[0]];
301 }
302
303 SegExpSharedPtr segexp[3];
305 m_base[0]->GetBasisKey(), GetGeom2D()->GetEdge(0));
306
308 {
309 for (i = 1; i < 3; i++)
310 {
312 m_base[i != 0]->GetBasisKey(), GetGeom2D()->GetEdge(i));
313 }
314 }
315 else // interploate using edge 0 GLL distribution
316 {
317 for (i = 1; i < 3; i++)
318 {
320 m_base[0]->GetBasisKey(), GetGeom2D()->GetEdge(i));
321
322 LibUtilities::Interp1D(m_base[1]->GetPointsKey(), physEdge[i],
323 m_base[0]->GetPointsKey(), physEdge[i]);
324 }
325 npoints[1] = npoints[0];
326 }
327
329 Array<OneD, int> signArray;
331 // define an orientation to get EdgeToElmtMapping from Cartesian data
332 StdRegions::Orientation orient[3] = {
334
335 for (i = 0; i < 3; i++)
336 {
337 segexp[i]->FwdTransBndConstrained(physEdge[i], coeffEdge[i]);
338
339 // this orient goes with the one above and so could
340 // probably set both to eForwards
341 GetTraceToElementMap(i, mapArray, signArray, orient[i]);
342 for (j = 0; j < nmodes[i != 0]; j++)
343 {
344 sign = (NekDouble)signArray[j];
345 outarray[mapArray[j]] = sign * coeffEdge[i][j];
346 }
347 }
348
349 int nBoundaryDofs = NumBndryCoeffs();
350 int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
351
352 if (nInteriorDofs > 0)
353 {
356
358 *this);
359 MassMatrixOp(outarray, tmp0, stdmasskey);
360 IProductWRTBase(inarray, tmp1);
361
362 Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
363
364 // get Mass matrix inverse (only of interior DOF)
365 // use block (1,1) of the static condensed system
366 // note: this block alreay contains the inverse matrix
367 MatrixKey masskey(StdRegions::eMass, DetShapeType(), *this);
368 DNekScalMatSharedPtr matsys =
369 (m_staticCondMatrixManager[masskey])->GetBlock(1, 1);
370
371 Array<OneD, NekDouble> rhs(nInteriorDofs);
372 Array<OneD, NekDouble> result(nInteriorDofs);
373
374 GetInteriorMap(mapArray);
375
376 for (i = 0; i < nInteriorDofs; i++)
377 {
378 rhs[i] = tmp1[mapArray[i]];
379 }
380
381 Blas::Dgemv('N', nInteriorDofs, nInteriorDofs, matsys->Scale(),
382 &((matsys->GetOwnedMatrix())->GetPtr())[0], nInteriorDofs,
383 rhs.data(), 1, 0.0, result.data(), 1);
384
385 for (i = 0; i < nInteriorDofs; i++)
386 {
387 outarray[mapArray[i]] = result[i];
388 }
389 }
390}
391
393 Array<OneD, NekDouble> &outarray)
394{
395 IProductWRTBase_SumFac(inarray, outarray);
396}
397
399 const Array<OneD, const NekDouble> &inarray,
400 Array<OneD, NekDouble> &outarray)
401{
402 IProductWRTDerivBase_SumFac(dir, inarray, outarray);
403}
404
406 const Array<OneD, const NekDouble> &inarray,
407 Array<OneD, NekDouble> &outarray, bool multiplybyweights)
408{
409 int nquad0 = m_base[0]->GetNumPoints();
410 int nquad1 = m_base[1]->GetNumPoints();
411 int order0 = m_base[0]->GetNumModes();
412
413 if (multiplybyweights)
414 {
415 Array<OneD, NekDouble> tmp(nquad0 * nquad1 + nquad1 * order0);
416 Array<OneD, NekDouble> wsp(tmp + nquad0 * nquad1);
417
418 MultiplyByQuadratureMetric(inarray, tmp);
420 m_base[1]->GetBdata(), tmp, outarray, wsp);
421 }
422 else
423 {
424 Array<OneD, NekDouble> wsp(+nquad1 * order0);
425
427 m_base[1]->GetBdata(), inarray, outarray,
428 wsp);
429 }
430}
431
433 const int dir, const Array<OneD, const NekDouble> &inarray,
434 Array<OneD, NekDouble> &outarray)
435{
436 int nquad0 = m_base[0]->GetNumPoints();
437 int nquad1 = m_base[1]->GetNumPoints();
438 int nqtot = nquad0 * nquad1;
439 int nmodes0 = m_base[0]->GetNumModes();
440 int wspsize = max(max(nqtot, m_ncoeffs), nquad1 * nmodes0);
441
442 Array<OneD, NekDouble> tmp0(4 * wspsize);
443 Array<OneD, NekDouble> tmp1(tmp0 + wspsize);
444 Array<OneD, NekDouble> tmp2(tmp0 + 2 * wspsize);
445 Array<OneD, NekDouble> tmp3(tmp0 + 3 * wspsize);
446
448 tmp2D[0] = tmp1;
449 tmp2D[1] = tmp2;
450
451 TriExp::v_AlignVectorToCollapsedDir(dir, inarray, tmp2D);
452
453 MultiplyByQuadratureMetric(tmp1, tmp1);
454 MultiplyByQuadratureMetric(tmp2, tmp2);
455
456 IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
457 tmp1, tmp3, tmp0);
458 IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
459 tmp2, outarray, tmp0);
460 Vmath::Vadd(m_ncoeffs, tmp3, 1, outarray, 1, outarray, 1);
461}
462
464 const int dir, const Array<OneD, const NekDouble> &inarray,
466{
467 ASSERTL1((dir == 0) || (dir == 1) || (dir == 2), "Invalid direction.");
468 ASSERTL1((dir == 2) ? (m_geom->GetCoordim() == 3) : true,
469 "Invalid direction.");
470
471 int nquad0 = m_base[0]->GetNumPoints();
472 int nquad1 = m_base[1]->GetNumPoints();
473 int nqtot = nquad0 * nquad1;
474 int nmodes0 = m_base[0]->GetNumModes();
475 int wspsize = max(max(nqtot, m_ncoeffs), nquad1 * nmodes0);
476
478 m_metricinfo->GetDerivFactors(GetPointsKeys());
479
480 Array<OneD, NekDouble> tmp0(wspsize);
481 Array<OneD, NekDouble> tmp3(wspsize);
482 Array<OneD, NekDouble> gfac0(wspsize);
483 Array<OneD, NekDouble> gfac1(wspsize);
484
485 Array<OneD, NekDouble> tmp1 = outarray[0];
486 Array<OneD, NekDouble> tmp2 = outarray[1];
487
488 const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
489 const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
490
491 // set up geometric factor: 2/(1-z1)
492 for (int i = 0; i < nquad1; ++i)
493 {
494 gfac0[i] = 2.0 / (1 - z1[i]);
495 }
496 for (int i = 0; i < nquad0; ++i)
497 {
498 gfac1[i] = 0.5 * (1 + z0[i]);
499 }
500
501 for (int i = 0; i < nquad1; ++i)
502 {
503 Vmath::Smul(nquad0, gfac0[i], &inarray[0] + i * nquad0, 1,
504 &tmp0[0] + i * nquad0, 1);
505 }
506
507 for (int i = 0; i < nquad1; ++i)
508 {
509 Vmath::Vmul(nquad0, &gfac1[0], 1, &tmp0[0] + i * nquad0, 1,
510 &tmp1[0] + i * nquad0, 1);
511 }
512
513 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
514 {
515 Vmath::Vmul(nqtot, &df[2 * dir][0], 1, &tmp0[0], 1, &tmp0[0], 1);
516 Vmath::Vmul(nqtot, &df[2 * dir + 1][0], 1, &tmp1[0], 1, &tmp1[0], 1);
517 Vmath::Vmul(nqtot, &df[2 * dir + 1][0], 1, &inarray[0], 1, &tmp2[0], 1);
518 }
519 else
520 {
521 Vmath::Smul(nqtot, df[2 * dir][0], tmp0, 1, tmp0, 1);
522 Vmath::Smul(nqtot, df[2 * dir + 1][0], tmp1, 1, tmp1, 1);
523 Vmath::Smul(nqtot, df[2 * dir + 1][0], inarray, 1, tmp2, 1);
524 }
525 Vmath::Vadd(nqtot, tmp0, 1, tmp1, 1, tmp1, 1);
526}
527
529 const Array<OneD, const NekDouble> &direction,
530 const Array<OneD, const NekDouble> &inarray,
531 Array<OneD, NekDouble> &outarray)
532{
533 IProductWRTDirectionalDerivBase_SumFac(direction, inarray, outarray);
534}
535
536/**
537 * @brief Directinoal Derivative in the modal space in the dir
538 * direction of varcoeffs.
539 */
541 const Array<OneD, const NekDouble> &direction,
542 const Array<OneD, const NekDouble> &inarray,
543 Array<OneD, NekDouble> &outarray)
544{
545 int i;
546 int shapedim = 2;
547 int nquad0 = m_base[0]->GetNumPoints();
548 int nquad1 = m_base[1]->GetNumPoints();
549 int nqtot = nquad0 * nquad1;
550 int nmodes0 = m_base[0]->GetNumModes();
551 int wspsize = max(max(nqtot, m_ncoeffs), nquad1 * nmodes0);
552
554 m_metricinfo->GetDerivFactors(GetPointsKeys());
555
556 Array<OneD, NekDouble> tmp0(6 * wspsize);
557 Array<OneD, NekDouble> tmp1(tmp0 + wspsize);
558 Array<OneD, NekDouble> tmp2(tmp0 + 2 * wspsize);
559 Array<OneD, NekDouble> tmp3(tmp0 + 3 * wspsize);
560 Array<OneD, NekDouble> gfac0(tmp0 + 4 * wspsize);
561 Array<OneD, NekDouble> gfac1(tmp0 + 5 * wspsize);
562
563 const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
564 const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
565
566 // set up geometric factor: 2/(1-z1)
567 for (i = 0; i < nquad1; ++i)
568 {
569 gfac0[i] = 2.0 / (1 - z1[i]);
570 }
571 for (i = 0; i < nquad0; ++i)
572 {
573 gfac1[i] = 0.5 * (1 + z0[i]);
574 }
575 for (i = 0; i < nquad1; ++i)
576 {
577 Vmath::Smul(nquad0, gfac0[i], &inarray[0] + i * nquad0, 1,
578 &tmp0[0] + i * nquad0, 1);
579 }
580 for (i = 0; i < nquad1; ++i)
581 {
582 Vmath::Vmul(nquad0, &gfac1[0], 1, &tmp0[0] + i * nquad0, 1,
583 &tmp1[0] + i * nquad0, 1);
584 }
585
586 // Compute gmat \cdot e^j
587 Array<OneD, Array<OneD, NekDouble>> dfdir(shapedim);
588 Expansion::ComputeGmatcdotMF(df, direction, dfdir);
589
590 Vmath::Vmul(nqtot, &dfdir[0][0], 1, &tmp0[0], 1, &tmp0[0], 1);
591 Vmath::Vmul(nqtot, &dfdir[1][0], 1, &tmp1[0], 1, &tmp1[0], 1);
592 Vmath::Vmul(nqtot, &dfdir[1][0], 1, &inarray[0], 1, &tmp2[0], 1);
593
594 Vmath::Vadd(nqtot, &tmp0[0], 1, &tmp1[0], 1, &tmp1[0], 1);
595
596 MultiplyByQuadratureMetric(tmp1, tmp1);
597 MultiplyByQuadratureMetric(tmp2, tmp2);
598
599 IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
600 tmp1, tmp3, tmp0);
601 IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
602 tmp2, outarray, tmp0);
603 Vmath::Vadd(m_ncoeffs, tmp3, 1, outarray, 1, outarray, 1);
604}
605
609 Array<OneD, NekDouble> &outarray)
610{
611 int nq = m_base[0]->GetNumPoints() * m_base[1]->GetNumPoints();
613
615 GetLeftAdjacentElementExp()->GetTraceNormal(
617
618 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
619 {
620 Vmath::Vvtvvtp(nq, &normals[0][0], 1, &Fx[0], 1, &normals[1][0], 1,
621 &Fy[0], 1, &Fn[0], 1);
622 Vmath::Vvtvp(nq, &normals[2][0], 1, &Fz[0], 1, &Fn[0], 1, &Fn[0], 1);
623 }
624 else
625 {
626 Vmath::Svtsvtp(nq, normals[0][0], &Fx[0], 1, normals[1][0], &Fy[0], 1,
627 &Fn[0], 1);
628 Vmath::Svtvp(nq, normals[2][0], &Fz[0], 1, &Fn[0], 1, &Fn[0], 1);
629 }
630
631 IProductWRTBase(Fn, outarray);
632}
633
635 const Array<OneD, const Array<OneD, NekDouble>> &Fvec,
636 Array<OneD, NekDouble> &outarray)
637{
638 NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
639}
640
642{
643
645 m_base[0]->GetBasisKey(), m_base[1]->GetBasisKey());
646}
647
649{
651 m_base[0]->GetPointsKey());
653 m_base[1]->GetPointsKey());
654
656 bkey1);
657}
658
661{
662 int i;
663
664 ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 && Lcoords[1] >= -1.0 &&
665 Lcoords[1] <= 1.0,
666 "Local coordinates are not in region [-1,1]");
667
668 m_geom->FillGeom();
669
670 for (i = 0; i < m_geom->GetCoordim(); ++i)
671 {
672 coords[i] = m_geom->GetCoord(i, Lcoords);
673 }
674}
675
677 Array<OneD, NekDouble> &coords_1,
678 Array<OneD, NekDouble> &coords_2)
679{
680 Expansion::v_GetCoords(coords_0, coords_1, coords_2);
681}
682
683/**
684 * Given the local cartesian coordinate \a Lcoord evaluate the
685 * value of physvals at this point by calling through to the
686 * StdExpansion method
687 */
689 const Array<OneD, const NekDouble> &Lcoord,
690 const Array<OneD, const NekDouble> &physvals)
691{
692 // Evaluate point in local (eta) coordinates.
693 return StdExpansion2D::v_PhysEvaluate(Lcoord, physvals);
694}
695
697 const Array<OneD, const NekDouble> &physvals)
698{
700
701 ASSERTL0(m_geom, "m_geom not defined");
702 m_geom->GetLocCoords(coord, Lcoord);
703
704 return StdExpansion2D::v_PhysEvaluate(Lcoord, physvals);
705}
706
708 const Array<OneD, NekDouble> &coord,
709 const Array<OneD, const NekDouble> &inarray,
710 std::array<NekDouble, 3> &firstOrderDerivs)
711{
712 Array<OneD, NekDouble> Lcoord(2);
713 ASSERTL0(m_geom, "m_geom not defined");
714 m_geom->GetLocCoords(coord, Lcoord);
715 return StdTriExp::v_PhysEvalFirstDeriv(Lcoord, inarray, firstOrderDerivs);
716}
717
719 const int edge, const StdRegions::StdExpansionSharedPtr &EdgeExp,
720 const Array<OneD, const NekDouble> &inarray,
722{
723 int nquad0 = m_base[0]->GetNumPoints();
724 int nquad1 = m_base[1]->GetNumPoints();
725 int nt = 0;
726 // Extract in Cartesian direction because we have to deal with
727 // e.g. Gauss-Radau points.
728 switch (edge)
729 {
730 case 0:
731 Vmath::Vcopy(nquad0, &(inarray[0]), 1, &(outarray[0]), 1);
732 nt = nquad0;
733 break;
734 case 1:
735 Vmath::Vcopy(nquad1, &(inarray[0]) + (nquad0 - 1), nquad0,
736 &(outarray[0]), 1);
737 nt = nquad1;
738 break;
739 case 2:
740 Vmath::Vcopy(nquad1, &(inarray[0]), nquad0, &(outarray[0]), 1);
741 nt = nquad1;
742 break;
743 default:
744 ASSERTL0(false, "edge value (< 3) is out of range");
745 break;
746 }
747
748 ASSERTL1(EdgeExp->GetBasis(0)->GetPointsType() ==
750 "Edge expansion should be GLL");
751
752 // Interpolate if required
753 if (m_base[edge ? 1 : 0]->GetPointsKey() !=
754 EdgeExp->GetBasis(0)->GetPointsKey())
755 {
756 Array<OneD, NekDouble> outtmp(max(nquad0, nquad1));
757
758 Vmath::Vcopy(nt, outarray, 1, outtmp, 1);
759
760 LibUtilities::Interp1D(m_base[edge ? 1 : 0]->GetPointsKey(), outtmp,
761 EdgeExp->GetBasis(0)->GetPointsKey(), outarray);
762 }
763
764 if (orient == StdRegions::eNoOrientation)
765 {
766 orient = GetTraceOrient(edge);
767 }
768
769 // Reverse data if necessary
770 if (orient == StdRegions::eBackwards)
771 {
772 Vmath::Reverse(EdgeExp->GetNumPoints(0), &outarray[0], 1, &outarray[0],
773 1);
774 }
775}
776
778 [[maybe_unused]] const int edge,
779 [[maybe_unused]] Array<OneD, NekDouble> &outarray)
780{
781 ASSERTL0(false, "Routine not implemented for triangular elements");
782}
783
784void TriExp::v_GetTracePhysMap(const int edge, Array<OneD, int> &outarray)
785{
786 int nquad0 = m_base[0]->GetNumPoints();
787 int nquad1 = m_base[1]->GetNumPoints();
788
789 // Get points in Cartesian orientation
790 switch (edge)
791 {
792 case 0:
793 outarray = Array<OneD, int>(nquad0);
794 for (int i = 0; i < nquad0; ++i)
795 {
796 outarray[i] = i;
797 }
798 break;
799 case 1:
800 outarray = Array<OneD, int>(nquad1);
801 for (int i = 0; i < nquad1; ++i)
802 {
803 outarray[i] = (nquad0 - 1) + i * nquad0;
804 }
805 break;
806 case 2:
807 outarray = Array<OneD, int>(nquad1);
808 for (int i = 0; i < nquad1; ++i)
809 {
810 outarray[i] = i * nquad0;
811 }
812 break;
813 default:
814 ASSERTL0(false, "edge value (< 3) is out of range");
815 break;
816 }
817}
818
820{
821 int i;
822 const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
823 GetGeom()->GetMetricInfo();
824
826 for (i = 0; i < ptsKeys.size(); ++i)
827 {
828 // Need at least 2 points for computing normals
829 if (ptsKeys[i].GetNumPoints() == 1)
830 {
831 LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
832 ptsKeys[i] = pKey;
833 }
834 }
835
836 const SpatialDomains::GeomType type = geomFactors->GetGtype();
838 geomFactors->GetDerivFactors(ptsKeys);
839 const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
840
841 // The points of normals should follow trace basis, not local basis.
843
844 int nqe = tobasis.GetNumPoints();
845 int dim = GetCoordim();
846
849 for (i = 0; i < dim; ++i)
850 {
851 normal[i] = Array<OneD, NekDouble>(nqe);
852 }
853
854 size_t nqb = nqe;
855 size_t nbnd = edge;
858
859 // Regular geometry case
860 if ((type == SpatialDomains::eRegular) ||
862 {
863 NekDouble fac;
864 // Set up normals
865 switch (edge)
866 {
867 case 0:
868 for (i = 0; i < GetCoordim(); ++i)
869 {
870 Vmath::Fill(nqe, -df[2 * i + 1][0], normal[i], 1);
871 }
872 break;
873 case 1:
874 for (i = 0; i < GetCoordim(); ++i)
875 {
876 Vmath::Fill(nqe, df[2 * i + 1][0] + df[2 * i][0], normal[i],
877 1);
878 }
879 break;
880 case 2:
881 for (i = 0; i < GetCoordim(); ++i)
882 {
883 Vmath::Fill(nqe, -df[2 * i][0], normal[i], 1);
884 }
885 break;
886 default:
887 ASSERTL0(false, "Edge is out of range (edge < 3)");
888 }
889
890 // normalise
891 fac = 0.0;
892 for (i = 0; i < GetCoordim(); ++i)
893 {
894 fac += normal[i][0] * normal[i][0];
895 }
896 fac = 1.0 / sqrt(fac);
897
898 Vmath::Fill(nqb, fac, length, 1);
899
900 for (i = 0; i < GetCoordim(); ++i)
901 {
902 Vmath::Smul(nqe, fac, normal[i], 1, normal[i], 1);
903 }
904 }
905 else // Set up deformed normals
906 {
907 int j;
908
909 int nquad0 = ptsKeys[0].GetNumPoints();
910 int nquad1 = ptsKeys[1].GetNumPoints();
911
913
914 Array<OneD, NekDouble> normals(GetCoordim() * max(nquad0, nquad1), 0.0);
915 Array<OneD, NekDouble> edgejac(GetCoordim() * max(nquad0, nquad1), 0.0);
916
917 // Extract Jacobian along edges and recover local
918 // derivates (dx/dr) for polynomial interpolation by
919 // multiplying m_gmat by jacobian
920 switch (edge)
921 {
922 case 0:
923 for (j = 0; j < nquad0; ++j)
924 {
925 edgejac[j] = jac[j];
926 for (i = 0; i < GetCoordim(); ++i)
927 {
928 normals[i * nquad0 + j] =
929 -df[2 * i + 1][j] * edgejac[j];
930 }
931 }
932 from_key = ptsKeys[0];
933 break;
934 case 1:
935 for (j = 0; j < nquad1; ++j)
936 {
937 edgejac[j] = jac[nquad0 * j + nquad0 - 1];
938 for (i = 0; i < GetCoordim(); ++i)
939 {
940 normals[i * nquad1 + j] =
941 (df[2 * i][nquad0 * j + nquad0 - 1] +
942 df[2 * i + 1][nquad0 * j + nquad0 - 1]) *
943 edgejac[j];
944 }
945 }
946 from_key = ptsKeys[1];
947 break;
948 case 2:
949 for (j = 0; j < nquad1; ++j)
950 {
951 edgejac[j] = jac[nquad0 * j];
952 for (i = 0; i < GetCoordim(); ++i)
953 {
954 normals[i * nquad1 + j] =
955 -df[2 * i][nquad0 * j] * edgejac[j];
956 }
957 }
958 from_key = ptsKeys[1];
959 break;
960 default:
961 ASSERTL0(false, "edge is out of range (edge < 3)");
962 }
963
964 int nq = from_key.GetNumPoints();
965 Array<OneD, NekDouble> work(nqe, 0.0);
966
967 // interpolate Jacobian and invert
968 LibUtilities::Interp1D(from_key, jac, tobasis.GetPointsKey(), work);
969 Vmath::Sdiv(nqe, 1.0, &work[0], 1, &work[0], 1);
970
971 // interpolate
972 for (i = 0; i < GetCoordim(); ++i)
973 {
974 LibUtilities::Interp1D(from_key, &normals[i * nq],
975 tobasis.GetPointsKey(), &normal[i][0]);
976 Vmath::Vmul(nqe, work, 1, normal[i], 1, normal[i], 1);
977 }
978
979 // normalise normal vectors
980 Vmath::Zero(nqe, work, 1);
981 for (i = 0; i < GetCoordim(); ++i)
982 {
983 Vmath::Vvtvp(nqe, normal[i], 1, normal[i], 1, work, 1, work, 1);
984 }
985
986 Vmath::Vsqrt(nqe, work, 1, work, 1);
987 Vmath::Sdiv(nqe, 1.0, work, 1, work, 1);
988
989 Vmath::Vcopy(nqb, work, 1, length, 1);
990
991 for (i = 0; i < GetCoordim(); ++i)
992 {
993 Vmath::Vmul(nqe, normal[i], 1, work, 1, normal[i], 1);
994 }
995 }
996
997 if (GetGeom()->GetEorient(edge) == StdRegions::eBackwards)
998 {
999 for (i = 0; i < GetCoordim(); ++i)
1000 {
1001 if (geomFactors->GetGtype() == SpatialDomains::eDeformed)
1002 {
1003 Vmath::Reverse(nqe, normal[i], 1, normal[i], 1);
1004 }
1005 }
1006 }
1007}
1008
1010 const NekDouble *data, const std::vector<unsigned int> &nummodes,
1011 const int mode_offset, NekDouble *coeffs,
1012 [[maybe_unused]] std::vector<LibUtilities::BasisType> &fromType)
1013{
1014 int data_order0 = nummodes[mode_offset];
1015 int fillorder0 = min(m_base[0]->GetNumModes(), data_order0);
1016 int data_order1 = nummodes[mode_offset + 1];
1017 int order1 = m_base[1]->GetNumModes();
1018 int fillorder1 = min(order1, data_order1);
1019
1020 switch (m_base[0]->GetBasisType())
1021 {
1024 {
1025 int i;
1026 int cnt = 0;
1027 int cnt1 = 0;
1028
1031 "Extraction routine not set up for this basis");
1032
1033 Vmath::Zero(m_ncoeffs, coeffs, 1);
1034 for (i = 0; i < fillorder0; ++i)
1035 {
1036 Vmath::Vcopy(fillorder1 - i, &data[cnt], 1, &coeffs[cnt1], 1);
1037 cnt += data_order1 - i;
1038 cnt1 += order1 - i;
1039 }
1040 }
1041 break;
1042 default:
1043 ASSERTL0(false, "basis is either not set up or not hierarchicial");
1044 }
1045}
1046
1048{
1049 return GetGeom2D()->GetEorient(edge);
1050}
1051
1053{
1054 DNekMatSharedPtr returnval;
1055 switch (mkey.GetMatrixType())
1056 {
1064 returnval = Expansion2D::v_GenMatrix(mkey);
1065 break;
1066 default:
1067 returnval = StdTriExp::v_GenMatrix(mkey);
1068 break;
1069 }
1070
1071 return returnval;
1072}
1073
1075{
1076 LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1077 LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1080
1081 return tmp->GetStdMatrix(mkey);
1082}
1083
1085{
1086 return m_matrixManager[mkey];
1087}
1088
1090{
1091 m_matrixManager.DeleteObject(mkey);
1092}
1093
1095{
1096 return m_staticCondMatrixManager[mkey];
1097}
1098
1100{
1101 m_staticCondMatrixManager.DeleteObject(mkey);
1102}
1103
1105 Array<OneD, NekDouble> &outarray,
1106 const StdRegions::StdMatrixKey &mkey)
1107{
1108 StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
1109}
1110
1112 Array<OneD, NekDouble> &outarray,
1113 const StdRegions::StdMatrixKey &mkey)
1114{
1115 TriExp::LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
1116}
1117
1118void TriExp::v_LaplacianMatrixOp(const int k1, const int k2,
1119 const Array<OneD, const NekDouble> &inarray,
1120 Array<OneD, NekDouble> &outarray,
1121 const StdRegions::StdMatrixKey &mkey)
1122{
1123 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1124}
1125
1127 const Array<OneD, const NekDouble> &inarray,
1128 Array<OneD, NekDouble> &outarray,
1129 const StdRegions::StdMatrixKey &mkey)
1130{
1131 StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
1132}
1133
1135 const Array<OneD, const NekDouble> &inarray,
1136 Array<OneD, NekDouble> &outarray, const StdRegions::StdMatrixKey &mkey)
1137{
1138 StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray, outarray, mkey);
1139}
1140
1142 const Array<OneD, const NekDouble> &inarray,
1143 Array<OneD, NekDouble> &outarray, const StdRegions::StdMatrixKey &mkey)
1144{
1145 StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray, outarray, mkey);
1146}
1147
1149 Array<OneD, NekDouble> &outarray,
1150 const StdRegions::StdMatrixKey &mkey)
1151{
1152 TriExp::HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
1153}
1154
1156 const Array<OneD, const NekDouble> &inarray,
1158{
1159 if (m_metrics.count(eMetricLaplacian00) == 0)
1160 {
1162 }
1163
1164 int nquad0 = m_base[0]->GetNumPoints();
1165 int nquad1 = m_base[1]->GetNumPoints();
1166 int nqtot = nquad0 * nquad1;
1167 int nmodes0 = m_base[0]->GetNumModes();
1168 int nmodes1 = m_base[1]->GetNumModes();
1169 int wspsize =
1170 max(max(max(nqtot, m_ncoeffs), nquad1 * nmodes0), nquad0 * nmodes1);
1171
1172 ASSERTL1(wsp.size() >= 3 * wspsize, "Workspace is of insufficient size.");
1173
1174 const Array<OneD, const NekDouble> &base0 = m_base[0]->GetBdata();
1175 const Array<OneD, const NekDouble> &base1 = m_base[1]->GetBdata();
1176 const Array<OneD, const NekDouble> &dbase0 = m_base[0]->GetDbdata();
1177 const Array<OneD, const NekDouble> &dbase1 = m_base[1]->GetDbdata();
1178 const Array<OneD, const NekDouble> &metric00 =
1180 const Array<OneD, const NekDouble> &metric01 =
1182 const Array<OneD, const NekDouble> &metric11 =
1184
1185 // Allocate temporary storage
1186 Array<OneD, NekDouble> wsp0(wsp);
1187 Array<OneD, NekDouble> wsp1(wsp + wspsize);
1188 Array<OneD, NekDouble> wsp2(wsp + 2 * wspsize);
1189
1190 StdExpansion2D::PhysTensorDeriv(inarray, wsp1, wsp2);
1191
1192 // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1193 // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1194 // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1195 // especially for this purpose
1196 Vmath::Vvtvvtp(nqtot, &metric00[0], 1, &wsp1[0], 1, &metric01[0], 1,
1197 &wsp2[0], 1, &wsp0[0], 1);
1198 Vmath::Vvtvvtp(nqtot, &metric01[0], 1, &wsp1[0], 1, &metric11[0], 1,
1199 &wsp2[0], 1, &wsp2[0], 1);
1200
1201 // outarray = m = (D_xi1 * B)^T * k
1202 // wsp1 = n = (D_xi2 * B)^T * l
1203 IProductWRTBase_SumFacKernel(dbase0, base1, wsp0, outarray, wsp1);
1204 IProductWRTBase_SumFacKernel(base0, dbase1, wsp2, wsp1, wsp0);
1205
1206 // outarray = outarray + wsp1
1207 // = L * u_hat
1208 Vmath::Vadd(m_ncoeffs, wsp1.data(), 1, outarray.data(), 1, outarray.data(),
1209 1);
1210}
1211
1213{
1214 if (m_metrics.count(eMetricQuadrature) == 0)
1215 {
1217 }
1218
1219 unsigned int i, j;
1220 const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
1221 const unsigned int nqtot = GetTotPoints();
1222 const unsigned int dim = 2;
1223 const MetricType m[3][3] = {
1227
1228 Array<OneD, NekDouble> dEta_dXi[2] = {Array<OneD, NekDouble>(nqtot, 1.0),
1229 Array<OneD, NekDouble>(nqtot, 1.0)};
1230
1231 for (i = 0; i < dim; ++i)
1232 {
1233 for (j = i; j < dim; ++j)
1234 {
1235 m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1236 }
1237 }
1238
1239 const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
1240 const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
1241 const unsigned int nquad0 = m_base[0]->GetNumPoints();
1242 const unsigned int nquad1 = m_base[1]->GetNumPoints();
1244 m_metricinfo->GetDerivFactors(GetPointsKeys());
1245
1246 for (i = 0; i < nquad1; i++)
1247 {
1248 Blas::Dscal(nquad0, 2.0 / (1 - z1[i]), &dEta_dXi[0][0] + i * nquad0, 1);
1249 Blas::Dscal(nquad0, 2.0 / (1 - z1[i]), &dEta_dXi[1][0] + i * nquad0, 1);
1250 }
1251 for (i = 0; i < nquad0; i++)
1252 {
1253 Blas::Dscal(nquad1, 0.5 * (1 + z0[i]), &dEta_dXi[1][0] + i, nquad0);
1254 }
1255
1256 Array<OneD, NekDouble> tmp(nqtot);
1257 if ((type == SpatialDomains::eRegular ||
1259 {
1260 Vmath::Smul(nqtot, df[0][0], &dEta_dXi[0][0], 1, &tmp[0], 1);
1261 Vmath::Svtvp(nqtot, df[1][0], &dEta_dXi[1][0], 1, &tmp[0], 1, &tmp[0],
1262 1);
1263
1264 Vmath::Vmul(nqtot, &tmp[0], 1, &tmp[0], 1,
1266 Vmath::Smul(nqtot, df[1][0], &tmp[0], 1,
1268
1269 Vmath::Smul(nqtot, df[2][0], &dEta_dXi[0][0], 1, &tmp[0], 1);
1270 Vmath::Svtvp(nqtot, df[3][0], &dEta_dXi[1][0], 1, &tmp[0], 1, &tmp[0],
1271 1);
1272
1273 Vmath::Vvtvp(nqtot, &tmp[0], 1, &tmp[0], 1,
1276 Vmath::Svtvp(nqtot, df[3][0], &tmp[0], 1,
1279
1280 if (GetCoordim() == 3)
1281 {
1282 Vmath::Smul(nqtot, df[4][0], &dEta_dXi[0][0], 1, &tmp[0], 1);
1283 Vmath::Svtvp(nqtot, df[5][0], &dEta_dXi[1][0], 1, &tmp[0], 1,
1284 &tmp[0], 1);
1285
1286 Vmath::Vvtvp(nqtot, &tmp[0], 1, &tmp[0], 1,
1289 Vmath::Svtvp(nqtot, df[5][0], &tmp[0], 1,
1292 }
1293
1294 NekDouble g2 = df[1][0] * df[1][0] + df[3][0] * df[3][0];
1295 if (GetCoordim() == 3)
1296 {
1297 g2 += df[5][0] * df[5][0];
1298 }
1299 Vmath::Fill(nqtot, g2, &m_metrics[eMetricLaplacian11][0], 1);
1300 }
1301 else
1302 {
1303
1304 Vmath::Vmul(nqtot, &df[0][0], 1, &dEta_dXi[0][0], 1, &tmp[0], 1);
1305 Vmath::Vvtvp(nqtot, &df[1][0], 1, &dEta_dXi[1][0], 1, &tmp[0], 1,
1306 &tmp[0], 1);
1307
1308 Vmath::Vmul(nqtot, &tmp[0], 1, &tmp[0], 1,
1310 Vmath::Vmul(nqtot, &df[1][0], 1, &tmp[0], 1,
1312 Vmath::Vmul(nqtot, &df[1][0], 1, &df[1][0], 1,
1314
1315 Vmath::Vmul(nqtot, &df[2][0], 1, &dEta_dXi[0][0], 1, &tmp[0], 1);
1316 Vmath::Vvtvp(nqtot, &df[3][0], 1, &dEta_dXi[1][0], 1, &tmp[0], 1,
1317 &tmp[0], 1);
1318
1319 Vmath::Vvtvp(nqtot, &tmp[0], 1, &tmp[0], 1,
1322 Vmath::Vvtvp(nqtot, &df[3][0], 1, &tmp[0], 1,
1325 Vmath::Vvtvp(nqtot, &df[3][0], 1, &df[3][0], 1,
1328
1329 if (GetCoordim() == 3)
1330 {
1331 Vmath::Vmul(nqtot, &df[4][0], 1, &dEta_dXi[0][0], 1, &tmp[0], 1);
1332 Vmath::Vvtvp(nqtot, &df[5][0], 1, &dEta_dXi[1][0], 1, &tmp[0], 1,
1333 &tmp[0], 1);
1334
1335 Vmath::Vvtvp(nqtot, &tmp[0], 1, &tmp[0], 1,
1338 Vmath::Vvtvp(nqtot, &df[5][0], 1, &tmp[0], 1,
1341 Vmath::Vvtvp(nqtot, &df[5][0], 1, &df[5][0], 1,
1344 }
1345 }
1346
1347 for (unsigned int i = 0; i < dim; ++i)
1348 {
1349 for (unsigned int j = i; j < dim; ++j)
1350 {
1352 }
1353 }
1354}
1355
1356/**
1357 * Function is used to compute exactly the advective numerical flux on
1358 * theinterface of two elements with different expansions, hence an
1359 * appropriate number of Gauss points has to be used. The number of
1360 * Gauss points has to be equal to the number used by the highest
1361 * polynomial degree of the two adjacent elements. Furthermore, this
1362 * function is used to compute the sensor value in each element.
1363 *
1364 * @param numMin Is the reduced polynomial order
1365 * @param inarray Input array of coefficients
1366 * @param dumpVar Output array of reduced coefficients.
1367 */
1369 const Array<OneD, const NekDouble> &inarray,
1370 Array<OneD, NekDouble> &outarray)
1371{
1372 int n_coeffs = inarray.size();
1373 int nquad0 = m_base[0]->GetNumPoints();
1374 int nquad1 = m_base[1]->GetNumPoints();
1375 int nqtot = nquad0 * nquad1;
1376 int nmodes0 = m_base[0]->GetNumModes();
1377 int nmodes1 = m_base[1]->GetNumModes();
1378 int numMin2 = nmodes0, i;
1379
1380 Array<OneD, NekDouble> coeff(n_coeffs, 0.0);
1381 Array<OneD, NekDouble> phys_tmp(nqtot, 0.0);
1382 Array<OneD, NekDouble> tmp, tmp2;
1383
1384 const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
1385 const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
1386
1387 LibUtilities::BasisKey b0(m_base[0]->GetBasisType(), nmodes0, Pkey0);
1388 LibUtilities::BasisKey b1(m_base[1]->GetBasisType(), nmodes1, Pkey1);
1389 LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A, nmodes0, Pkey0);
1390 LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_B, nmodes1, Pkey1);
1391
1392 // Check if it is also possible to use the same InterCoeff routine
1393 // which is also used for Quadrilateral and Hexagonal shaped
1394 // elements
1395
1396 // For now, set up the used basis on the standard element to
1397 // calculate the phys values, set up the orthogonal basis to do a
1398 // forward transform, to obtain the coefficients in orthogonal
1399 // coefficient space
1400 StdRegions::StdTriExpSharedPtr m_OrthoTriExp;
1402
1405 bortho0, bortho1);
1406
1407 m_TriExp->BwdTrans(inarray, phys_tmp);
1408 m_OrthoTriExp->FwdTrans(phys_tmp, coeff);
1409
1410 for (i = 0; i < n_coeffs; i++)
1411 {
1412 if (i == numMin)
1413 {
1414 coeff[i] = 0.0;
1415 numMin += numMin2 - 1;
1416 numMin2 -= 1.0;
1417 }
1418 }
1419
1420 m_OrthoTriExp->BwdTrans(coeff, phys_tmp);
1421 m_TriExp->FwdTrans(phys_tmp, outarray);
1422}
1423
1425 const StdRegions::StdMatrixKey &mkey)
1426{
1427 int nq = GetTotPoints();
1428
1429 // Calculate sqrt of the Jacobian
1431 Array<OneD, NekDouble> sqrt_jac(nq);
1432 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1433 {
1434 Vmath::Vsqrt(nq, jac, 1, sqrt_jac, 1);
1435 }
1436 else
1437 {
1438 Vmath::Fill(nq, sqrt(jac[0]), sqrt_jac, 1);
1439 }
1440
1441 // Multiply array by sqrt(Jac)
1442 Vmath::Vmul(nq, sqrt_jac, 1, array, 1, array, 1);
1443
1444 // Apply std region filter
1445 StdTriExp::v_SVVLaplacianFilter(array, mkey);
1446
1447 // Divide by sqrt(Jac)
1448 Vmath::Vdiv(nq, array, 1, sqrt_jac, 1, array, 1);
1449}
1450
1451/** @brief: This method gets all of the factors which are
1452 required as part of the Gradient Jump Penalty
1453 stabilisation and involves the product of the normal and
1454 geometric factors along the element trace.
1455*/
1457 Array<OneD, Array<OneD, NekDouble>> &d0factors,
1458 Array<OneD, Array<OneD, NekDouble>> &d1factors,
1459 [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &d2factors)
1460{
1461 int nquad0 = GetNumPoints(0);
1462 int nquad1 = GetNumPoints(1);
1463
1465 m_metricinfo->GetDerivFactors(GetPointsKeys());
1466
1467 if (d0factors.size() != 3)
1468 {
1469 d0factors = Array<OneD, Array<OneD, NekDouble>>(3);
1470 d1factors = Array<OneD, Array<OneD, NekDouble>>(3);
1471 }
1472
1473 if (d0factors[0].size() != nquad0)
1474 {
1475 d0factors[0] = Array<OneD, NekDouble>(nquad0);
1476 d1factors[0] = Array<OneD, NekDouble>(nquad0);
1477 }
1478
1479 if (d0factors[1].size() != nquad1)
1480 {
1481 d0factors[1] = Array<OneD, NekDouble>(nquad1);
1482 d0factors[2] = Array<OneD, NekDouble>(nquad1);
1483 d1factors[1] = Array<OneD, NekDouble>(nquad1);
1484 d1factors[2] = Array<OneD, NekDouble>(nquad1);
1485 }
1486
1487 // Outwards normals
1489 GetTraceNormal(0);
1491 GetTraceNormal(1);
1493 GetTraceNormal(2);
1494
1495 int ncoords = normal_0.size();
1496
1497 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1498 {
1499
1500 // d xi_2/dx n_x
1501 for (int i = 0; i < nquad0; ++i)
1502 {
1503 d1factors[0][i] = df[1][i] * normal_0[0][i];
1504 }
1505
1506 // d xi_1/dx n_x
1507 for (int i = 0; i < nquad1; ++i)
1508 {
1509 d0factors[1][i] = df[0][(i + 1) * nquad0 - 1] * normal_1[0][i];
1510 d0factors[2][i] = df[0][i * nquad0] * normal_2[0][i];
1511 }
1512
1513 for (int n = 1; n < ncoords; ++n)
1514 {
1515 // d xi_2/dy n_y
1516 // needs checking for 3D coords
1517 for (int i = 0; i < nquad0; ++i)
1518 {
1519 d1factors[0][i] += df[2 * n + 1][i] * normal_0[n][i];
1520 }
1521
1522 // d xi_1/dy n_y
1523 // needs checking for 3D coords
1524 for (int i = 0; i < nquad1; ++i)
1525 {
1526 d0factors[1][i] +=
1527 df[2 * n][(i + 1) * nquad0 - 1] * normal_1[n][i];
1528 d0factors[2][i] += df[2 * n][i * nquad0] * normal_2[n][i];
1529 }
1530 }
1531
1532 // d0 factors
1533 // d xi_1/dx n_x
1534 for (int i = 0; i < nquad0; ++i)
1535 {
1536 d0factors[0][i] = df[0][i] * normal_0[0][i];
1537 }
1538
1539 // d xi_2/dx n_x
1540 for (int i = 0; i < nquad1; ++i)
1541 {
1542 d1factors[1][i] = df[1][(i + 1) * nquad0 - 1] * normal_1[0][i];
1543 d1factors[2][i] = df[1][i * nquad0] * normal_2[0][i];
1544 }
1545
1546 for (int n = 1; n < ncoords; ++n)
1547 {
1548 // d xi_1/dy n_y
1549 // needs checking for 3D coords
1550 for (int i = 0; i < nquad0; ++i)
1551 {
1552 d0factors[0][i] += df[2 * n][i] * normal_0[n][i];
1553 }
1554
1555 // d xi_2/dy n_y
1556 // needs checking for 3D coords
1557 for (int i = 0; i < nquad1; ++i)
1558 {
1559 d1factors[1][i] +=
1560 df[2 * n + 1][(i + 1) * nquad0 - 1] * normal_1[n][i];
1561 d1factors[2][i] += df[2 * n + 1][i * nquad0] * normal_2[n][i];
1562 }
1563 }
1564 }
1565 else
1566 {
1567 // d xi_2/dx n_x
1568 for (int i = 0; i < nquad0; ++i)
1569 {
1570 d1factors[0][i] = df[1][0] * normal_0[0][i];
1571 }
1572
1573 // d xi_1/dx n_x
1574 for (int i = 0; i < nquad1; ++i)
1575 {
1576 d0factors[1][i] = df[0][0] * normal_1[0][i];
1577 d0factors[2][i] = df[0][0] * normal_2[0][i];
1578 }
1579
1580 for (int n = 1; n < ncoords; ++n)
1581 {
1582 // d xi_2/dy n_y
1583 // needs checking for 3D coords
1584 for (int i = 0; i < nquad0; ++i)
1585 {
1586 d1factors[0][i] += df[2 * n + 1][0] * normal_0[n][i];
1587 }
1588
1589 // d xi_1/dy n_y
1590 // needs checking for 3D coords
1591 for (int i = 0; i < nquad1; ++i)
1592 {
1593 d0factors[1][i] += df[2 * n][0] * normal_1[n][i];
1594 d0factors[2][i] += df[2 * n][0] * normal_2[n][i];
1595 }
1596 }
1597
1598 // d1factors
1599 // d xi_1/dx n_x
1600 for (int i = 0; i < nquad0; ++i)
1601 {
1602 d0factors[0][i] = df[0][0] * normal_0[0][i];
1603 }
1604
1605 // d xi_2/dx n_x
1606 for (int i = 0; i < nquad1; ++i)
1607 {
1608 d1factors[1][i] = df[1][0] * normal_1[0][i];
1609 d1factors[2][i] = df[1][0] * normal_2[0][i];
1610 }
1611
1612 for (int n = 1; n < ncoords; ++n)
1613 {
1614 // d xi_1/dy n_y
1615 // needs checking for 3D coords
1616 for (int i = 0; i < nquad0; ++i)
1617 {
1618 d0factors[0][i] += df[2 * n][0] * normal_0[n][i];
1619 }
1620
1621 // d xi_2/dy n_y
1622 // needs checking for 3D coords
1623 for (int i = 0; i < nquad1; ++i)
1624 {
1625 d1factors[1][i] += df[2 * n + 1][0] * normal_1[n][i];
1626 d1factors[2][i] += df[2 * n + 1][0] * normal_2[n][i];
1627 }
1628 }
1629 }
1630}
1631} // namespace Nektar::LocalRegions
#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 sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:47
Describes the specification for a Basis.
Definition: Basis.h:45
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:120
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:137
Defines a specification for a set of points.
Definition: Points.h:50
size_t GetNumPoints() const
Definition: Points.h:85
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
Definition: Expansion2D.h:164
std::map< int, NormalVector > m_traceNormals
Definition: Expansion.h:276
std::map< int, Array< OneD, NekDouble > > m_elmtBndNormDirElmtLen
the element length in each element boundary(Vertex, edge or face) normal direction calculated based o...
Definition: Expansion.h:286
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:167
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:273
ExpansionSharedPtr GetLeftAdjacentElementExp() const
Definition: Expansion.h:441
void ComputeGmatcdotMF(const Array< TwoD, const NekDouble > &df, const Array< OneD, const NekDouble > &direction, Array< OneD, Array< OneD, NekDouble > > &dfdir)
Definition: Expansion.cpp:607
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
Definition: Expansion.cpp:534
int GetLeftAdjacentElementTrace() const
Definition: Expansion.h:454
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:274
StdRegions::Orientation GetTraceOrient(int trace)
Definition: Expansion.h:168
const NormalVector & GetTraceNormal(const int id)
Definition: Expansion.cpp:251
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: TriExp.cpp:432
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: TriExp.h:245
StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const override
Definition: TriExp.cpp:641
void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: TriExp.cpp:1148
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
Definition: TriExp.cpp:707
void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: TriExp.cpp:1126
NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coord, const Array< OneD, const NekDouble > &physvals) override
This function evaluates the expansion at a single (arbitrary) point of the domain.
Definition: TriExp.cpp:696
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: TriExp.cpp:398
void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray) override
Definition: TriExp.cpp:463
void v_DropLocMatrix(const MatrixKey &mkey) override
Definition: TriExp.cpp:1089
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
Definition: TriExp.cpp:405
void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp) override
Definition: TriExp.cpp:1155
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
Definition: TriExp.cpp:676
DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey) override
Definition: TriExp.cpp:1094
void v_GetTracePhysMap(const int edge, Array< OneD, int > &outarray) override
Definition: TriExp.cpp:784
DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey) override
Definition: TriExp.cpp:1074
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: TriExp.h:247
void v_IProductWRTDirectionalDerivBase(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: TriExp.cpp:528
void v_PhysDirectionalDeriv(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &direction, Array< OneD, NekDouble > &out) override
Physical derivative along a direction vector.
Definition: TriExp.cpp:184
StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const override
Definition: TriExp.cpp:648
void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: TriExp.cpp:1111
void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords) override
Definition: TriExp.cpp:659
DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey) override
Definition: TriExp.cpp:1084
StdRegions::Orientation v_GetTraceOrient(int edge) override
Definition: TriExp.cpp:1047
NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals) override
Definition: TriExp.cpp:688
void v_ComputeTraceNormal(const int edge) override
Definition: TriExp.cpp:819
void v_GetTracePhysVals(const int edge, const StdRegions::StdExpansionSharedPtr &EdgeExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient) override
Definition: TriExp.cpp:718
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
Definition: TriExp.cpp:1052
void v_DropLocStaticCondMatrix(const MatrixKey &mkey) override
Definition: TriExp.cpp:1099
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: TriExp.cpp:392
void v_GetTraceQFactors(const int edge, Array< OneD, NekDouble > &outarray) override
Definition: TriExp.cpp:777
void v_WeakDirectionalDerivMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: TriExp.cpp:1134
void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: TriExp.cpp:265
void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: TriExp.cpp:1368
void v_ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int mode_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType) override
Definition: TriExp.cpp:1009
void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: TriExp.cpp:1104
void v_NormalTraceDerivFactors(Array< OneD, Array< OneD, NekDouble > > &factors, Array< OneD, Array< OneD, NekDouble > > &d0factors, Array< OneD, Array< OneD, NekDouble > > &d1factors) override
: This method gets all of the factors which are required as part of the Gradient Jump Penalty stabili...
Definition: TriExp.cpp:1456
void v_IProductWRTDirectionalDerivBase_SumFac(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Directinoal Derivative in the modal space in the dir direction of varcoeffs.
Definition: TriExp.cpp:540
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey) override
Definition: TriExp.cpp:1424
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: TriExp.cpp:94
void v_NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, const Array< OneD, const NekDouble > &Fz, Array< OneD, NekDouble > &outarray) override
Definition: TriExp.cpp:606
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: TriExp.cpp:249
void v_MassLevelCurvatureMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: TriExp.cpp:1141
TriExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const SpatialDomains::Geometry2DSharedPtr &geom)
Constructor using BasisKey class for quadrature points and order definition.
Definition: TriExp.cpp:45
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray) override
Integrates the specified function over the domain.
Definition: TriExp.cpp:71
void v_ComputeLaplacianMetric() override
Definition: TriExp.cpp:1212
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
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)
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:134
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:156
void HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
const LibUtilities::PointsKeyVector GetPointsKeys() const
void MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:761
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:622
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
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:205
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:218
void IProductWRTDirectionalDerivBase_SumFac(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:732
const LibUtilities::BasisKey GetTraceBasisKey(const int i, int k=-1, bool UseGLL=false) const
This function returns the basis key belonging to the i-th trace.
Definition: StdExpansion.h:301
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Definition: StdExpansion.h:858
Array< OneD, LibUtilities::BasisSharedPtr > m_base
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:83
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 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
void Interp1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
this function interpolates a 1D function evaluated at the quadrature points of the basis fbasis0 to ...
Definition: Interp.cpp:47
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:231
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:42
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:44
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
std::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:253
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:60
GeomType
Indicates the type of element geometry.
@ eRegular
Geometry is straight-sided with constant geometric factors.
@ eMovingRegular
Currently unused.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:62
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::shared_ptr< StdTriExp > StdTriExpSharedPtr
Definition: StdTriExp.h:218
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
Svtsvtp (scalar times vector plus scalar times vector):
Definition: Vmath.hpp:473
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvp (scalar times vector plus vector): z = alpha*x + y.
Definition: Vmath.hpp:396
void 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 Vdiv(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:126
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 Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:844
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.hpp:439
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