Nektar++
QuadExp.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: QuadExp.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: Expansion for quadrilateral elements.
32//
33///////////////////////////////////////////////////////////////////////////////
34
41#include <LocalRegions/SegExp.h>
42
43using namespace std;
44
46{
48 const LibUtilities::BasisKey &Bb,
50 : StdExpansion(Ba.GetNumModes() * Bb.GetNumModes(), 2, Ba, Bb),
51 StdExpansion2D(Ba.GetNumModes() * Bb.GetNumModes(), Ba, Bb),
52 StdQuadExp(Ba, Bb), Expansion(geom), Expansion2D(geom),
53 m_matrixManager(
54 std::bind(&Expansion2D::CreateMatrix, this, std::placeholders::_1),
55 std::string("QuadExpMatrix")),
56 m_staticCondMatrixManager(std::bind(&Expansion::CreateStaticCondMatrix,
57 this, std::placeholders::_1),
58 std::string("QuadExpStaticCondMatrix"))
59{
60}
61
63 : StdExpansion(T), StdExpansion2D(T), StdQuadExp(T), Expansion(T),
64 Expansion2D(T), m_matrixManager(T.m_matrixManager),
65 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
66{
67}
68
70{
71 int nquad0 = m_base[0]->GetNumPoints();
72 int nquad1 = m_base[1]->GetNumPoints();
74 NekDouble ival;
75 Array<OneD, NekDouble> tmp(nquad0 * nquad1);
76
77 // multiply inarray with Jacobian
78 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
79 {
80 Vmath::Vmul(nquad0 * nquad1, jac, 1, inarray, 1, tmp, 1);
81 }
82 else
83 {
84 Vmath::Smul(nquad0 * nquad1, jac[0], inarray, 1, tmp, 1);
85 }
86
87 // call StdQuadExp version;
88 ival = StdQuadExp::v_Integral(tmp);
89 return ival;
90}
91
96{
97 int nquad0 = m_base[0]->GetNumPoints();
98 int nquad1 = m_base[1]->GetNumPoints();
99 int nqtot = nquad0 * nquad1;
101 m_metricinfo->GetDerivFactors(GetPointsKeys());
102 Array<OneD, NekDouble> diff0(2 * nqtot);
103 Array<OneD, NekDouble> diff1(diff0 + nqtot);
104
105 StdQuadExp::v_PhysDeriv(inarray, diff0, diff1);
106
107 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
108 {
109 if (out_d0.size())
110 {
111 Vmath::Vmul(nqtot, df[0], 1, diff0, 1, out_d0, 1);
112 Vmath::Vvtvp(nqtot, df[1], 1, diff1, 1, out_d0, 1, out_d0, 1);
113 }
114
115 if (out_d1.size())
116 {
117 Vmath::Vmul(nqtot, df[2], 1, diff0, 1, out_d1, 1);
118 Vmath::Vvtvp(nqtot, df[3], 1, diff1, 1, out_d1, 1, out_d1, 1);
119 }
120
121 if (out_d2.size())
122 {
123 Vmath::Vmul(nqtot, df[4], 1, diff0, 1, out_d2, 1);
124 Vmath::Vvtvp(nqtot, df[5], 1, diff1, 1, out_d2, 1, out_d2, 1);
125 }
126 }
127 else // regular geometry
128 {
129 if (out_d0.size())
130 {
131 Vmath::Smul(nqtot, df[0][0], diff0, 1, out_d0, 1);
132 Blas::Daxpy(nqtot, df[1][0], diff1, 1, out_d0, 1);
133 }
134
135 if (out_d1.size())
136 {
137 Vmath::Smul(nqtot, df[2][0], diff0, 1, out_d1, 1);
138 Blas::Daxpy(nqtot, df[3][0], diff1, 1, out_d1, 1);
139 }
140
141 if (out_d2.size())
142 {
143 Vmath::Smul(nqtot, df[4][0], diff0, 1, out_d2, 1);
144 Blas::Daxpy(nqtot, df[5][0], diff1, 1, out_d2, 1);
145 }
146 }
147}
148
149void QuadExp::v_PhysDeriv(const int dir,
150 const Array<OneD, const NekDouble> &inarray,
151 Array<OneD, NekDouble> &outarray)
152{
153 switch (dir)
154 {
155 case 0:
156 {
157 v_PhysDeriv(inarray, outarray, NullNekDouble1DArray,
159 }
160 break;
161 case 1:
162 {
163 v_PhysDeriv(inarray, NullNekDouble1DArray, outarray,
165 }
166 break;
167 case 2:
168 {
170 outarray);
171 }
172 break;
173 default:
174 {
175 ASSERTL1(false, "input dir is out of range");
176 }
177 break;
178 }
179}
180
182 const Array<OneD, const NekDouble> &inarray,
184{
185 int nquad0 = m_base[0]->GetNumPoints();
186 int nquad1 = m_base[1]->GetNumPoints();
187 int nqtot = nquad0 * nquad1;
188
190 m_metricinfo->GetDerivFactors(GetPointsKeys());
191
192 Array<OneD, NekDouble> diff0(2 * nqtot);
193 Array<OneD, NekDouble> diff1(diff0 + nqtot);
194
195 StdQuadExp::v_PhysDeriv(inarray, diff0, diff1);
196
197 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
198 {
200
201 // d/dx_v^s = v_x*ds/dx + v_y*ds/dy + v_z*dx/dz
202 for (int i = 0; i < 2; ++i)
203 {
204 tangmat[i] = Array<OneD, NekDouble>(nqtot, 0.0);
205 for (int k = 0; k < (m_geom->GetCoordim()); ++k)
206 {
207 Vmath::Vvtvp(nqtot, &df[2 * k + i][0], 1, &direction[k * nqtot],
208 1, &tangmat[i][0], 1, &tangmat[i][0], 1);
209 }
210 }
211
212 /// D_v = d/dx_v^s + d/dx_v^r
213 if (out.size())
214 {
215 Vmath::Vmul(nqtot, &tangmat[0][0], 1, &diff0[0], 1, &out[0], 1);
216 Vmath::Vvtvp(nqtot, &tangmat[1][0], 1, &diff1[0], 1, &out[0], 1,
217 &out[0], 1);
218 }
219 }
220 else
221 {
223 "Wrong route");
224 }
225}
226
228 Array<OneD, NekDouble> &outarray)
229{
230 if ((m_base[0]->Collocation()) && (m_base[1]->Collocation()))
231 {
232 Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
233 }
234 else
235 {
236 IProductWRTBase(inarray, outarray);
237
238 // get Mass matrix inverse
240 DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
241
242 // copy inarray in case inarray == outarray
243 NekVector<NekDouble> in(m_ncoeffs, outarray, eCopy);
245
246 out = (*matsys) * in;
247 }
248}
249
251 const Array<OneD, const NekDouble> &inarray,
252 Array<OneD, NekDouble> &outarray)
253{
254 if ((m_base[0]->Collocation()) && (m_base[1]->Collocation()))
255 {
256 Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
257 }
258 else
259 {
260 int i, j;
261 int npoints[2] = {m_base[0]->GetNumPoints(), m_base[1]->GetNumPoints()};
262 int nmodes[2] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes()};
263
264 fill(outarray.get(), outarray.get() + m_ncoeffs, 0.0);
265
266 if (nmodes[0] == 1 && nmodes[1] == 1)
267 {
268 outarray[0] = inarray[0];
269 return;
270 }
271
272 Array<OneD, NekDouble> physEdge[4];
273 Array<OneD, NekDouble> coeffEdge[4];
274 StdRegions::Orientation orient[4];
275 for (i = 0; i < 4; i++)
276 {
277 physEdge[i] = Array<OneD, NekDouble>(npoints[i % 2]);
278 coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i % 2]);
279 orient[i] = GetTraceOrient(i);
280 }
281
282 for (i = 0; i < npoints[0]; i++)
283 {
284 physEdge[0][i] = inarray[i];
285 physEdge[2][i] = inarray[npoints[0] * (npoints[1] - 1) + i];
286 }
287
288 for (i = 0; i < npoints[1]; i++)
289 {
290 physEdge[1][i] = inarray[npoints[0] - 1 + i * npoints[0]];
291 physEdge[3][i] = inarray[i * npoints[0]];
292 }
293
294 for (i = 0; i < 4; i++)
295 {
296 if (orient[i] == StdRegions::eBackwards)
297 {
298 reverse((physEdge[i]).get(),
299 (physEdge[i]).get() + npoints[i % 2]);
300 }
301 }
302
303 SegExpSharedPtr segexp[4];
304 for (i = 0; i < 4; i++)
305 {
307 m_base[i % 2]->GetBasisKey(), GetGeom2D()->GetEdge(i));
308 }
309
311 Array<OneD, int> signArray;
313
314 for (i = 0; i < 4; i++)
315 {
316 segexp[i % 2]->FwdTransBndConstrained(physEdge[i], coeffEdge[i]);
317
318 GetTraceToElementMap(i, mapArray, signArray, orient[i]);
319 for (j = 0; j < nmodes[i % 2]; j++)
320 {
321 sign = (NekDouble)signArray[j];
322 outarray[mapArray[j]] = sign * coeffEdge[i][j];
323 }
324 }
325
326 int nBoundaryDofs = NumBndryCoeffs();
327 int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
328
329 if (nInteriorDofs > 0)
330 {
333
335 DetShapeType(), *this);
336 MassMatrixOp(outarray, tmp0, stdmasskey);
337 IProductWRTBase(inarray, tmp1);
338
339 Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
340
341 // get Mass matrix inverse (only of interior DOF)
342 // use block (1,1) of the static condensed system
343 // note: this block alreay contains the inverse matrix
344 MatrixKey masskey(StdRegions::eMass, DetShapeType(), *this);
345 DNekScalMatSharedPtr matsys =
346 (m_staticCondMatrixManager[masskey])->GetBlock(1, 1);
347
348 Array<OneD, NekDouble> rhs(nInteriorDofs);
349 Array<OneD, NekDouble> result(nInteriorDofs);
350
351 GetInteriorMap(mapArray);
352
353 for (i = 0; i < nInteriorDofs; i++)
354 {
355 rhs[i] = tmp1[mapArray[i]];
356 }
357
358 Blas::Dgemv('N', nInteriorDofs, nInteriorDofs, matsys->Scale(),
359 &((matsys->GetOwnedMatrix())->GetPtr())[0],
360 nInteriorDofs, rhs.get(), 1, 0.0, result.get(), 1);
361
362 for (i = 0; i < nInteriorDofs; i++)
363 {
364 outarray[mapArray[i]] = result[i];
365 }
366 }
367 }
368}
369
371 Array<OneD, NekDouble> &outarray)
372{
373 if (m_base[0]->Collocation() && m_base[1]->Collocation())
374 {
375 MultiplyByQuadratureMetric(inarray, outarray);
376 }
377 else
378 {
379 IProductWRTBase_SumFac(inarray, outarray);
380 }
381}
382
384 const int dir, const Array<OneD, const NekDouble> &inarray,
385 Array<OneD, NekDouble> &outarray)
386{
387 IProductWRTDerivBase_SumFac(dir, inarray, outarray);
388}
389
391 const Array<OneD, const NekDouble> &inarray,
392 Array<OneD, NekDouble> &outarray, bool multiplybyweights)
393{
394 int nquad0 = m_base[0]->GetNumPoints();
395 int nquad1 = m_base[1]->GetNumPoints();
396 int order0 = m_base[0]->GetNumModes();
397
398 if (multiplybyweights)
399 {
400 Array<OneD, NekDouble> tmp(nquad0 * nquad1 + nquad1 * order0);
401 Array<OneD, NekDouble> wsp(tmp + nquad0 * nquad1);
402
403 MultiplyByQuadratureMetric(inarray, tmp);
404 StdQuadExp::IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
405 m_base[1]->GetBdata(), tmp,
406 outarray, wsp, true, true);
407 }
408 else
409 {
410 Array<OneD, NekDouble> wsp(nquad1 * order0);
411
412 StdQuadExp::IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
413 m_base[1]->GetBdata(), inarray,
414 outarray, wsp, true, true);
415 }
416}
417
419 const int dir, const Array<OneD, const NekDouble> &inarray,
420 Array<OneD, NekDouble> &outarray)
421{
422 ASSERTL1((dir == 0) || (dir == 1) || (dir == 2), "Invalid direction.");
423 ASSERTL1((dir == 2) ? (m_geom->GetCoordim() == 3) : true,
424 "Invalid direction.");
425
426 int nquad0 = m_base[0]->GetNumPoints();
427 int nquad1 = m_base[1]->GetNumPoints();
428 int nqtot = nquad0 * nquad1;
429 int nmodes0 = m_base[0]->GetNumModes();
430
431 Array<OneD, NekDouble> tmp1(2 * nqtot + m_ncoeffs + nmodes0 * nquad1);
432 Array<OneD, NekDouble> tmp2(tmp1 + nqtot);
433 Array<OneD, NekDouble> tmp3(tmp1 + 2 * nqtot);
434 Array<OneD, NekDouble> tmp4(tmp1 + 2 * nqtot + m_ncoeffs);
435
437 tmp2D[0] = tmp1;
438 tmp2D[1] = tmp2;
439
440 QuadExp::v_AlignVectorToCollapsedDir(dir, inarray, tmp2D);
441
442 MultiplyByQuadratureMetric(tmp1, tmp1);
443 MultiplyByQuadratureMetric(tmp2, tmp2);
444
445 IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
446 tmp1, tmp3, tmp4, false, true);
447 IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
448 tmp2, outarray, tmp4, true, false);
449 Vmath::Vadd(m_ncoeffs, tmp3, 1, outarray, 1, outarray, 1);
450}
451
453 const int dir, const Array<OneD, const NekDouble> &inarray,
455{
456 ASSERTL1((dir == 0) || (dir == 1) || (dir == 2), "Invalid direction.");
457 ASSERTL1((dir == 2) ? (m_geom->GetCoordim() == 3) : true,
458 "Invalid direction.");
459
460 int nquad0 = m_base[0]->GetNumPoints();
461 int nquad1 = m_base[1]->GetNumPoints();
462 int nqtot = nquad0 * nquad1;
463 int nmodes0 = m_base[0]->GetNumModes();
464
466 m_metricinfo->GetDerivFactors(GetPointsKeys());
467
468 Array<OneD, NekDouble> tmp1 = outarray[0];
469 Array<OneD, NekDouble> tmp2 = outarray[1];
471 Array<OneD, NekDouble> tmp4(nmodes0 * nquad1);
472
473 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
474 {
475 Vmath::Vmul(nqtot, &df[2 * dir][0], 1, inarray.get(), 1, tmp1.get(), 1);
476 Vmath::Vmul(nqtot, &df[2 * dir + 1][0], 1, inarray.get(), 1, tmp2.get(),
477 1);
478 }
479 else
480 {
481 Vmath::Smul(nqtot, df[2 * dir][0], inarray.get(), 1, tmp1.get(), 1);
482 Vmath::Smul(nqtot, df[2 * dir + 1][0], inarray.get(), 1, tmp2.get(), 1);
483 }
484}
485
490{
491 int nq = m_base[0]->GetNumPoints() * m_base[1]->GetNumPoints();
493
495 GetLeftAdjacentElementExp()->GetTraceNormal(
497
498 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
499 {
500 Vmath::Vvtvvtp(nq, &normals[0][0], 1, &Fx[0], 1, &normals[1][0], 1,
501 &Fy[0], 1, &Fn[0], 1);
502 Vmath::Vvtvp(nq, &normals[2][0], 1, &Fz[0], 1, &Fn[0], 1, &Fn[0], 1);
503 }
504 else
505 {
506 Vmath::Svtsvtp(nq, normals[0][0], &Fx[0], 1, normals[1][0], &Fy[0], 1,
507 &Fn[0], 1);
508 Vmath::Svtvp(nq, normals[2][0], &Fz[0], 1, &Fn[0], 1, &Fn[0], 1);
509 }
510
511 IProductWRTBase(Fn, outarray);
512}
513
515 const Array<OneD, const Array<OneD, NekDouble>> &Fvec,
516 Array<OneD, NekDouble> &outarray)
517{
518 NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
519}
520
522{
524 m_base[0]->GetBasisKey(), m_base[1]->GetBasisKey());
525}
526
528{
530 m_base[0]->GetPointsKey());
532 m_base[1]->GetPointsKey());
533
535 bkey1);
536}
537
539 Array<OneD, NekDouble> &coords_1,
540 Array<OneD, NekDouble> &coords_2)
541{
542 Expansion::v_GetCoords(coords_0, coords_1, coords_2);
543}
544
547{
548 int i;
549
550 ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 && Lcoords[1] >= -1.0 &&
551 Lcoords[1] <= 1.0,
552 "Local coordinates are not in region [-1,1]");
553
554 m_geom->FillGeom();
555 for (i = 0; i < m_geom->GetCoordim(); ++i)
556 {
557 coords[i] = m_geom->GetCoord(i, Lcoords);
558 }
559}
560
561/**
562 * Given the local cartesian coordinate \a Lcoord evaluate the
563 * value of physvals at this point by calling through to the
564 * StdExpansion method
565 */
567 const Array<OneD, const NekDouble> &Lcoord,
568 const Array<OneD, const NekDouble> &physvals)
569{
570 // Evaluate point in local (eta) coordinates.
571 return StdExpansion2D::v_PhysEvaluate(Lcoord, physvals);
572}
573
575 const Array<OneD, const NekDouble> &physvals)
576{
578
579 ASSERTL0(m_geom, "m_geom not defined");
580 m_geom->GetLocCoords(coord, Lcoord);
581
582 return StdExpansion2D::v_PhysEvaluate(Lcoord, physvals);
583}
584
586 const Array<OneD, const NekDouble> &inarray,
587 std::array<NekDouble, 3> &firstOrderDerivs)
588{
589 Array<OneD, NekDouble> Lcoord(2);
590 ASSERTL0(m_geom, "m_geom not defined");
591 m_geom->GetLocCoords(coord, Lcoord);
592 return StdQuadExp::v_PhysEvaluate(Lcoord, inarray, firstOrderDerivs);
593}
594
596 const int edge, const StdRegions::StdExpansionSharedPtr &EdgeExp,
597 const Array<OneD, const NekDouble> &inarray,
599{
600 int nquad0 = m_base[0]->GetNumPoints();
601 int nquad1 = m_base[1]->GetNumPoints();
602
603 // Implementation for all the basis except Gauss points
606 {
607 switch (edge)
608 {
609 case 0:
610 Vmath::Vcopy(nquad0, &(inarray[0]), 1, &(outarray[0]), 1);
611 break;
612 case 1:
613 Vmath::Vcopy(nquad1, &(inarray[0]) + (nquad0 - 1), nquad0,
614 &(outarray[0]), 1);
615 break;
616 case 2:
617 Vmath::Vcopy(nquad0, &(inarray[0]) + nquad0 * (nquad1 - 1), 1,
618 &(outarray[0]), 1);
619 break;
620 case 3:
621 Vmath::Vcopy(nquad1, &(inarray[0]), nquad0, &(outarray[0]), 1);
622 break;
623 default:
624 ASSERTL0(false, "edge value (< 3) is out of range");
625 break;
626 }
627 }
628 else
629 {
630 QuadExp::GetEdgeInterpVals(edge, inarray, outarray);
631 }
632
633 // Interpolate if required
634 if (m_base[edge % 2]->GetPointsKey() !=
635 EdgeExp->GetBasis(0)->GetPointsKey())
636 {
637 Array<OneD, NekDouble> outtmp(max(nquad0, nquad1));
638
639 outtmp = outarray;
640
641 LibUtilities::Interp1D(m_base[edge % 2]->GetPointsKey(), outtmp,
642 EdgeExp->GetBasis(0)->GetPointsKey(), outarray);
643 }
644
645 if (orient == StdRegions::eNoOrientation)
646 {
647 orient = GetTraceOrient(edge);
648 }
649 // Reverse data if necessary
650 if (orient == StdRegions::eBackwards)
651 {
652 Vmath::Reverse(EdgeExp->GetNumPoints(0), &outarray[0], 1, &outarray[0],
653 1);
654 }
655}
656
657void QuadExp::GetEdgeInterpVals(const int edge,
658 const Array<OneD, const NekDouble> &inarray,
659 Array<OneD, NekDouble> &outarray)
660{
661 int i;
662 int nq0 = m_base[0]->GetNumPoints();
663 int nq1 = m_base[1]->GetNumPoints();
664
667
669 *this, factors);
670
671 DNekScalMatSharedPtr mat_gauss = m_matrixManager[key];
672
673 switch (edge)
674 {
675 case 0:
676 {
677 for (i = 0; i < nq0; i++)
678 {
679 outarray[i] =
680 Blas::Ddot(nq1, mat_gauss->GetOwnedMatrix()->GetPtr().get(),
681 1, &inarray[i], nq0);
682 }
683 break;
684 }
685 case 1:
686 {
687 for (i = 0; i < nq1; i++)
688 {
689 outarray[i] =
690 Blas::Ddot(nq0, mat_gauss->GetOwnedMatrix()->GetPtr().get(),
691 1, &inarray[i * nq0], 1);
692 }
693 break;
694 }
695 case 2:
696 {
697 for (i = 0; i < nq0; i++)
698 {
699 outarray[i] =
700 Blas::Ddot(nq1, mat_gauss->GetOwnedMatrix()->GetPtr().get(),
701 1, &inarray[i], nq0);
702 }
703 break;
704 }
705 case 3:
706 {
707 for (i = 0; i < nq1; i++)
708 {
709 outarray[i] =
710 Blas::Ddot(nq0, mat_gauss->GetOwnedMatrix()->GetPtr().get(),
711 1, &inarray[i * nq0], 1);
712 }
713 break;
714 }
715 default:
716 ASSERTL0(false, "edge value (< 3) is out of range");
717 break;
718 }
719}
720
721void QuadExp::v_GetTracePhysMap(const int edge, Array<OneD, int> &outarray)
722{
723 int nquad0 = m_base[0]->GetNumPoints();
724 int nquad1 = m_base[1]->GetNumPoints();
725
726 // Get points in Cartesian orientation
727 switch (edge)
728 {
729 case 0:
730 outarray = Array<OneD, int>(nquad0);
731 for (int i = 0; i < nquad0; ++i)
732 {
733 outarray[i] = i;
734 }
735 break;
736 case 1:
737 outarray = Array<OneD, int>(nquad1);
738 for (int i = 0; i < nquad1; ++i)
739 {
740 outarray[i] = (nquad0 - 1) + i * nquad0;
741 }
742 break;
743 case 2:
744 outarray = Array<OneD, int>(nquad0);
745 for (int i = 0; i < nquad0; ++i)
746 {
747 outarray[i] = i + nquad0 * (nquad1 - 1);
748 }
749 break;
750 case 3:
751 outarray = Array<OneD, int>(nquad1);
752 for (int i = 0; i < nquad1; ++i)
753 {
754 outarray[i] = i * nquad0;
755 }
756 break;
757 default:
758 ASSERTL0(false, "edge value (< 3) is out of range");
759 break;
760 }
761}
762
763void QuadExp::v_GetTraceQFactors(const int edge,
764 Array<OneD, NekDouble> &outarray)
765{
766 int i;
767 int nquad0 = m_base[0]->GetNumPoints();
768 int nquad1 = m_base[1]->GetNumPoints();
769
771 const Array<OneD, const NekDouble> &jac = m_metricinfo->GetJac(ptsKeys);
773 m_metricinfo->GetDerivFactors(ptsKeys);
774
775 Array<OneD, NekDouble> j(max(nquad0, nquad1), 0.0);
776 Array<OneD, NekDouble> g0(max(nquad0, nquad1), 0.0);
777 Array<OneD, NekDouble> g1(max(nquad0, nquad1), 0.0);
778 Array<OneD, NekDouble> g2(max(nquad0, nquad1), 0.0);
779 Array<OneD, NekDouble> g3(max(nquad0, nquad1), 0.0);
780
781 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
782 {
783 // Implementation for all the basis except Gauss points
786 {
787 switch (edge)
788 {
789 case 0:
790 Vmath::Vcopy(nquad0, &(df[1][0]), 1, &(g1[0]), 1);
791 Vmath::Vcopy(nquad0, &(df[3][0]), 1, &(g3[0]), 1);
792 Vmath::Vcopy(nquad0, &(jac[0]), 1, &(j[0]), 1);
793
794 for (i = 0; i < nquad0; ++i)
795 {
796 outarray[i] =
797 j[i] * sqrt(g1[i] * g1[i] + g3[i] * g3[i]);
798 }
799 break;
800 case 1:
801 Vmath::Vcopy(nquad1, &(df[0][0]) + (nquad0 - 1), nquad0,
802 &(g0[0]), 1);
803
804 Vmath::Vcopy(nquad1, &(df[2][0]) + (nquad0 - 1), nquad0,
805 &(g2[0]), 1);
806
807 Vmath::Vcopy(nquad1, &(jac[0]) + (nquad0 - 1), nquad0,
808 &(j[0]), 1);
809
810 for (i = 0; i < nquad1; ++i)
811 {
812 outarray[i] =
813 j[i] * sqrt(g0[i] * g0[i] + g2[i] * g2[i]);
814 }
815 break;
816 case 2:
817
818 Vmath::Vcopy(nquad0, &(df[1][0]) + (nquad0 * (nquad1 - 1)),
819 1, &(g1[0]), 1);
820
821 Vmath::Vcopy(nquad0, &(df[3][0]) + (nquad0 * (nquad1 - 1)),
822 1, &(g3[0]), 1);
823
824 Vmath::Vcopy(nquad0, &(jac[0]) + (nquad0 * (nquad1 - 1)), 1,
825 &(j[0]), 1);
826
827 for (i = 0; i < nquad0; ++i)
828 {
829 outarray[i] =
830 j[i] * sqrt(g1[i] * g1[i] + g3[i] * g3[i]);
831 }
832 break;
833 case 3:
834
835 Vmath::Vcopy(nquad1, &(df[0][0]), nquad0, &(g0[0]), 1);
836 Vmath::Vcopy(nquad1, &(df[2][0]), nquad0, &(g2[0]), 1);
837 Vmath::Vcopy(nquad1, &(jac[0]), nquad0, &(j[0]), 1);
838
839 for (i = 0; i < nquad1; ++i)
840 {
841 outarray[i] =
842 j[i] * sqrt(g0[i] * g0[i] + g2[i] * g2[i]);
843 }
844 break;
845 default:
846 ASSERTL0(false, "edge value (< 3) is out of range");
847 break;
848 }
849 }
850 else
851 {
852 int nqtot = nquad0 * nquad1;
853 Array<OneD, NekDouble> tmp_gmat0(nqtot, 0.0);
854 Array<OneD, NekDouble> tmp_gmat1(nqtot, 0.0);
855 Array<OneD, NekDouble> tmp_gmat2(nqtot, 0.0);
856 Array<OneD, NekDouble> tmp_gmat3(nqtot, 0.0);
857 Array<OneD, NekDouble> g0_edge(max(nquad0, nquad1), 0.0);
858 Array<OneD, NekDouble> g1_edge(max(nquad0, nquad1), 0.0);
859 Array<OneD, NekDouble> g2_edge(max(nquad0, nquad1), 0.0);
860 Array<OneD, NekDouble> g3_edge(max(nquad0, nquad1), 0.0);
861 Array<OneD, NekDouble> jac_edge(max(nquad0, nquad1), 0.0);
862
863 switch (edge)
864 {
865 case 0:
866 Vmath::Vmul(nqtot, &(df[1][0]), 1, &jac[0], 1,
867 &(tmp_gmat1[0]), 1);
868 Vmath::Vmul(nqtot, &(df[3][0]), 1, &jac[0], 1,
869 &(tmp_gmat3[0]), 1);
870 QuadExp::GetEdgeInterpVals(edge, tmp_gmat1, g1_edge);
871 QuadExp::GetEdgeInterpVals(edge, tmp_gmat3, g3_edge);
872
873 for (i = 0; i < nquad0; ++i)
874 {
875 outarray[i] = sqrt(g1_edge[i] * g1_edge[i] +
876 g3_edge[i] * g3_edge[i]);
877 }
878 break;
879
880 case 1:
881 Vmath::Vmul(nqtot, &(df[0][0]), 1, &jac[0], 1,
882 &(tmp_gmat0[0]), 1);
883 Vmath::Vmul(nqtot, &(df[2][0]), 1, &jac[0], 1,
884 &(tmp_gmat2[0]), 1);
885 QuadExp::GetEdgeInterpVals(edge, tmp_gmat0, g0_edge);
886 QuadExp::GetEdgeInterpVals(edge, tmp_gmat2, g2_edge);
887
888 for (i = 0; i < nquad1; ++i)
889 {
890 outarray[i] = sqrt(g0_edge[i] * g0_edge[i] +
891 g2_edge[i] * g2_edge[i]);
892 }
893
894 break;
895 case 2:
896
897 Vmath::Vmul(nqtot, &(df[1][0]), 1, &jac[0], 1,
898 &(tmp_gmat1[0]), 1);
899 Vmath::Vmul(nqtot, &(df[3][0]), 1, &jac[0], 1,
900 &(tmp_gmat3[0]), 1);
901 QuadExp::GetEdgeInterpVals(edge, tmp_gmat1, g1_edge);
902 QuadExp::GetEdgeInterpVals(edge, tmp_gmat3, g3_edge);
903
904 for (i = 0; i < nquad0; ++i)
905 {
906 outarray[i] = sqrt(g1_edge[i] * g1_edge[i] +
907 g3_edge[i] * g3_edge[i]);
908 }
909
910 Vmath::Reverse(nquad0, &outarray[0], 1, &outarray[0], 1);
911
912 break;
913 case 3:
914 Vmath::Vmul(nqtot, &(df[0][0]), 1, &jac[0], 1,
915 &(tmp_gmat0[0]), 1);
916 Vmath::Vmul(nqtot, &(df[2][0]), 1, &jac[0], 1,
917 &(tmp_gmat2[0]), 1);
918 QuadExp::GetEdgeInterpVals(edge, tmp_gmat0, g0_edge);
919 QuadExp::GetEdgeInterpVals(edge, tmp_gmat2, g2_edge);
920
921 for (i = 0; i < nquad1; ++i)
922 {
923 outarray[i] = sqrt(g0_edge[i] * g0_edge[i] +
924 g2_edge[i] * g2_edge[i]);
925 }
926
927 Vmath::Reverse(nquad1, &outarray[0], 1, &outarray[0], 1);
928
929 break;
930 default:
931 ASSERTL0(false, "edge value (< 3) is out of range");
932 break;
933 }
934 }
935 }
936 else
937 {
938
939 switch (edge)
940 {
941 case 0:
942
943 for (i = 0; i < nquad0; ++i)
944 {
945 outarray[i] = jac[0] * sqrt(df[1][0] * df[1][0] +
946 df[3][0] * df[3][0]);
947 }
948 break;
949 case 1:
950 for (i = 0; i < nquad1; ++i)
951 {
952 outarray[i] = jac[0] * sqrt(df[0][0] * df[0][0] +
953 df[2][0] * df[2][0]);
954 }
955 break;
956 case 2:
957 for (i = 0; i < nquad0; ++i)
958 {
959 outarray[i] = jac[0] * sqrt(df[1][0] * df[1][0] +
960 df[3][0] * df[3][0]);
961 }
962 break;
963 case 3:
964 for (i = 0; i < nquad1; ++i)
965 {
966 outarray[i] = jac[0] * sqrt(df[0][0] * df[0][0] +
967 df[2][0] * df[2][0]);
968 }
969 break;
970 default:
971 ASSERTL0(false, "edge value (< 3) is out of range");
972 break;
973 }
974 }
975}
976
978{
979 int i;
980 const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
981 GetGeom()->GetMetricInfo();
982 SpatialDomains::GeomType type = geomFactors->GetGtype();
983
985 for (i = 0; i < ptsKeys.size(); ++i)
986 {
987 // Need at least 2 points for computing normals
988 if (ptsKeys[i].GetNumPoints() == 1)
989 {
990 LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
991 ptsKeys[i] = pKey;
992 }
993 }
994
996 geomFactors->GetDerivFactors(ptsKeys);
997 const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
998
999 // The points of normals should follow trace basis, not local basis.
1001
1002 int nqe = tobasis.GetNumPoints();
1003 int vCoordDim = GetCoordim();
1004
1007 for (i = 0; i < vCoordDim; ++i)
1008 {
1009 normal[i] = Array<OneD, NekDouble>(nqe);
1010 }
1011
1012 size_t nqb = nqe;
1013 size_t nbnd = edge;
1016
1017 // Regular geometry case
1018 if ((type == SpatialDomains::eRegular) ||
1020 {
1021 NekDouble fac;
1022 // Set up normals
1023 switch (edge)
1024 {
1025 case 0:
1026 for (i = 0; i < vCoordDim; ++i)
1027 {
1028 Vmath::Fill(nqe, -df[2 * i + 1][0], normal[i], 1);
1029 }
1030 break;
1031 case 1:
1032 for (i = 0; i < vCoordDim; ++i)
1033 {
1034 Vmath::Fill(nqe, df[2 * i][0], normal[i], 1);
1035 }
1036 break;
1037 case 2:
1038 for (i = 0; i < vCoordDim; ++i)
1039 {
1040 Vmath::Fill(nqe, df[2 * i + 1][0], normal[i], 1);
1041 }
1042 break;
1043 case 3:
1044 for (i = 0; i < vCoordDim; ++i)
1045 {
1046 Vmath::Fill(nqe, -df[2 * i][0], normal[i], 1);
1047 }
1048 break;
1049 default:
1050 ASSERTL0(false, "edge is out of range (edge < 4)");
1051 }
1052
1053 // normalise
1054 fac = 0.0;
1055 for (i = 0; i < vCoordDim; ++i)
1056 {
1057 fac += normal[i][0] * normal[i][0];
1058 }
1059 fac = 1.0 / sqrt(fac);
1060
1061 Vmath::Fill(nqb, fac, length, 1);
1062
1063 for (i = 0; i < vCoordDim; ++i)
1064 {
1065 Vmath::Smul(nqe, fac, normal[i], 1, normal[i], 1);
1066 }
1067 }
1068 else // Set up deformed normals
1069 {
1070 int j;
1071
1072 int nquad0 = ptsKeys[0].GetNumPoints();
1073 int nquad1 = ptsKeys[1].GetNumPoints();
1074
1075 LibUtilities::PointsKey from_key;
1076
1077 Array<OneD, NekDouble> normals(vCoordDim * max(nquad0, nquad1), 0.0);
1078 Array<OneD, NekDouble> edgejac(vCoordDim * max(nquad0, nquad1), 0.0);
1079
1080 // Extract Jacobian along edges and recover local
1081 // derivates (dx/dr) for polynomial interpolation by
1082 // multiplying m_gmat by jacobian
1083
1084 // Implementation for all the basis except Gauss points
1087 {
1088 switch (edge)
1089 {
1090 case 0:
1091 for (j = 0; j < nquad0; ++j)
1092 {
1093 edgejac[j] = jac[j];
1094 for (i = 0; i < vCoordDim; ++i)
1095 {
1096 normals[i * nquad0 + j] =
1097 -df[2 * i + 1][j] * edgejac[j];
1098 }
1099 }
1100 from_key = ptsKeys[0];
1101 break;
1102 case 1:
1103 for (j = 0; j < nquad1; ++j)
1104 {
1105 edgejac[j] = jac[nquad0 * j + nquad0 - 1];
1106 for (i = 0; i < vCoordDim; ++i)
1107 {
1108 normals[i * nquad1 + j] =
1109 df[2 * i][nquad0 * j + nquad0 - 1] * edgejac[j];
1110 }
1111 }
1112 from_key = ptsKeys[1];
1113 break;
1114 case 2:
1115 for (j = 0; j < nquad0; ++j)
1116 {
1117 edgejac[j] = jac[nquad0 * (nquad1 - 1) + j];
1118 for (i = 0; i < vCoordDim; ++i)
1119 {
1120 normals[i * nquad0 + j] =
1121 (df[2 * i + 1][nquad0 * (nquad1 - 1) + j]) *
1122 edgejac[j];
1123 }
1124 }
1125 from_key = ptsKeys[0];
1126 break;
1127 case 3:
1128 for (j = 0; j < nquad1; ++j)
1129 {
1130 edgejac[j] = jac[nquad0 * j];
1131 for (i = 0; i < vCoordDim; ++i)
1132 {
1133 normals[i * nquad1 + j] =
1134 -df[2 * i][nquad0 * j] * edgejac[j];
1135 }
1136 }
1137 from_key = ptsKeys[1];
1138 break;
1139 default:
1140 ASSERTL0(false, "edge is out of range (edge < 3)");
1141 }
1142 }
1143 else
1144 {
1145 int nqtot = nquad0 * nquad1;
1146 Array<OneD, NekDouble> tmp_gmat(nqtot, 0.0);
1147 Array<OneD, NekDouble> tmp_gmat_edge(nqe, 0.0);
1148
1149 switch (edge)
1150 {
1151 case 0:
1152 for (j = 0; j < nquad0; ++j)
1153 {
1154 for (i = 0; i < vCoordDim; ++i)
1155 {
1156 Vmath::Vmul(nqtot, &(df[2 * i + 1][0]), 1, &jac[0],
1157 1, &(tmp_gmat[0]), 1);
1158 QuadExp::GetEdgeInterpVals(edge, tmp_gmat,
1159 tmp_gmat_edge);
1160 normals[i * nquad0 + j] = -tmp_gmat_edge[j];
1161 }
1162 }
1163 from_key = ptsKeys[0];
1164 break;
1165 case 1:
1166 for (j = 0; j < nquad1; ++j)
1167 {
1168 for (i = 0; i < vCoordDim; ++i)
1169 {
1170 Vmath::Vmul(nqtot, &(df[2 * i][0]), 1, &jac[0], 1,
1171 &(tmp_gmat[0]), 1);
1172 QuadExp::GetEdgeInterpVals(edge, tmp_gmat,
1173 tmp_gmat_edge);
1174 normals[i * nquad1 + j] = tmp_gmat_edge[j];
1175 }
1176 }
1177 from_key = ptsKeys[1];
1178 break;
1179 case 2:
1180 for (j = 0; j < nquad0; ++j)
1181 {
1182 for (i = 0; i < vCoordDim; ++i)
1183 {
1184 Vmath::Vmul(nqtot, &(df[2 * i + 1][0]), 1, &jac[0],
1185 1, &(tmp_gmat[0]), 1);
1186 QuadExp::GetEdgeInterpVals(edge, tmp_gmat,
1187 tmp_gmat_edge);
1188 normals[i * nquad0 + j] = tmp_gmat_edge[j];
1189 }
1190 }
1191 from_key = ptsKeys[0];
1192 break;
1193 case 3:
1194 for (j = 0; j < nquad1; ++j)
1195 {
1196 for (i = 0; i < vCoordDim; ++i)
1197 {
1198 Vmath::Vmul(nqtot, &(df[2 * i][0]), 1, &jac[0], 1,
1199 &(tmp_gmat[0]), 1);
1200 QuadExp::GetEdgeInterpVals(edge, tmp_gmat,
1201 tmp_gmat_edge);
1202 normals[i * nquad1 + j] = -tmp_gmat_edge[j];
1203 }
1204 }
1205 from_key = ptsKeys[1];
1206 break;
1207 default:
1208 ASSERTL0(false, "edge is out of range (edge < 3)");
1209 }
1210 }
1211
1212 int nq = from_key.GetNumPoints();
1213 Array<OneD, NekDouble> work(nqe, 0.0);
1214
1215 // interpolate Jacobian and invert
1216 LibUtilities::Interp1D(from_key, jac, tobasis.GetPointsKey(), work);
1217 Vmath::Sdiv(nqe, 1.0, &work[0], 1, &work[0], 1);
1218
1219 // interpolate
1220 for (i = 0; i < GetCoordim(); ++i)
1221 {
1222 LibUtilities::Interp1D(from_key, &normals[i * nq],
1223 tobasis.GetPointsKey(), &normal[i][0]);
1224 Vmath::Vmul(nqe, work, 1, normal[i], 1, normal[i], 1);
1225 }
1226
1227 // normalise normal vectors
1228 Vmath::Zero(nqe, work, 1);
1229 for (i = 0; i < GetCoordim(); ++i)
1230 {
1231 Vmath::Vvtvp(nqe, normal[i], 1, normal[i], 1, work, 1, work, 1);
1232 }
1233
1234 Vmath::Vsqrt(nqe, work, 1, work, 1);
1235 Vmath::Sdiv(nqe, 1.0, work, 1, work, 1);
1236
1237 Vmath::Vcopy(nqb, work, 1, length, 1);
1238
1239 for (i = 0; i < GetCoordim(); ++i)
1240 {
1241 Vmath::Vmul(nqe, normal[i], 1, work, 1, normal[i], 1);
1242 }
1243 }
1244 if (GetGeom()->GetEorient(edge) == StdRegions::eBackwards)
1245 {
1246 for (i = 0; i < vCoordDim; ++i)
1247 {
1248 if (geomFactors->GetGtype() == SpatialDomains::eDeformed)
1249 {
1250 Vmath::Reverse(nqe, normal[i], 1, normal[i], 1);
1251 }
1252 }
1253 }
1254}
1255
1257 const NekDouble *data, const std::vector<unsigned int> &nummodes,
1258 int mode_offset, NekDouble *coeffs,
1259 std::vector<LibUtilities::BasisType> &fromType)
1260{
1261 int data_order0 = nummodes[mode_offset];
1262 int fillorder0 = std::min(m_base[0]->GetNumModes(), data_order0);
1263
1264 int data_order1 = nummodes[mode_offset + 1];
1265 int order1 = m_base[1]->GetNumModes();
1266 int fillorder1 = min(order1, data_order1);
1267
1268 // Check if same basis
1269 if (fromType[0] != m_base[0]->GetBasisType() ||
1270 fromType[1] != m_base[1]->GetBasisType())
1271 {
1272 // Construct a quad with the appropriate basis type at our
1273 // quadrature points, and one more to do a forwards
1274 // transform. We can then copy the output to coeffs.
1275 StdRegions::StdQuadExp tmpQuad(
1276 LibUtilities::BasisKey(fromType[0], data_order0,
1277 m_base[0]->GetPointsKey()),
1278 LibUtilities::BasisKey(fromType[1], data_order1,
1279 m_base[1]->GetPointsKey()));
1280 StdRegions::StdQuadExp tmpQuad2(m_base[0]->GetBasisKey(),
1281 m_base[1]->GetBasisKey());
1282
1283 Array<OneD, const NekDouble> tmpData(tmpQuad.GetNcoeffs(), data);
1284 Array<OneD, NekDouble> tmpBwd(tmpQuad2.GetTotPoints());
1285 Array<OneD, NekDouble> tmpOut(tmpQuad2.GetNcoeffs());
1286
1287 tmpQuad.BwdTrans(tmpData, tmpBwd);
1288 tmpQuad2.FwdTrans(tmpBwd, tmpOut);
1289 Vmath::Vcopy(tmpOut.size(), &tmpOut[0], 1, coeffs, 1);
1290
1291 return;
1292 }
1293
1294 switch (m_base[0]->GetBasisType())
1295 {
1297 {
1298 int i;
1299 int cnt = 0;
1300 int cnt1 = 0;
1301
1303 "Extraction routine not set up for this basis");
1304
1305 Vmath::Zero(m_ncoeffs, coeffs, 1);
1306 for (i = 0; i < fillorder0; ++i)
1307 {
1308 Vmath::Vcopy(fillorder1, data + cnt, 1, coeffs + cnt1, 1);
1309 cnt += data_order1;
1310 cnt1 += order1;
1311 }
1312 }
1313 break;
1315 {
1316 LibUtilities::PointsKey p0(nummodes[0],
1318 LibUtilities::PointsKey p1(nummodes[1],
1320 LibUtilities::PointsKey t0(m_base[0]->GetNumModes(),
1322 LibUtilities::PointsKey t1(m_base[1]->GetNumModes(),
1324 LibUtilities::Interp2D(p0, p1, data, t0, t1, coeffs);
1325 }
1326 break;
1328 {
1329 // Assume that input is also Gll_Lagrange but no way to check;
1330 LibUtilities::PointsKey p0(nummodes[0],
1332 LibUtilities::PointsKey p1(nummodes[1],
1334 LibUtilities::PointsKey t0(m_base[0]->GetNumModes(),
1336 LibUtilities::PointsKey t1(m_base[1]->GetNumModes(),
1338 LibUtilities::Interp2D(p0, p1, data, t0, t1, coeffs);
1339 }
1340 break;
1341 default:
1342 ASSERTL0(false, "basis is either not set up or not hierarchicial");
1343 }
1344}
1345
1347{
1348 return m_geom->GetEorient(edge);
1349}
1350
1352{
1353 DNekMatSharedPtr returnval;
1354 switch (mkey.GetMatrixType())
1355 {
1363 returnval = Expansion2D::v_GenMatrix(mkey);
1364 break;
1365 default:
1366 returnval = StdQuadExp::v_GenMatrix(mkey);
1367 }
1368 return returnval;
1369}
1370
1372 const StdRegions::StdMatrixKey &mkey)
1373{
1374 LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1375 LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1378 return tmp->GetStdMatrix(mkey);
1379}
1380
1382{
1383 return m_matrixManager[mkey];
1384}
1385
1387{
1388 m_matrixManager.DeleteObject(mkey);
1389}
1390
1392{
1393 return m_staticCondMatrixManager[mkey];
1394}
1395
1397{
1398 m_staticCondMatrixManager.DeleteObject(mkey);
1399}
1400
1402 Array<OneD, NekDouble> &outarray,
1403 const StdRegions::StdMatrixKey &mkey)
1404{
1405 StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
1406}
1407
1409 Array<OneD, NekDouble> &outarray,
1410 const StdRegions::StdMatrixKey &mkey)
1411{
1412 QuadExp::LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
1413}
1414
1415void QuadExp::v_LaplacianMatrixOp(const int k1, const int k2,
1416 const Array<OneD, const NekDouble> &inarray,
1417 Array<OneD, NekDouble> &outarray,
1418 const StdRegions::StdMatrixKey &mkey)
1419{
1420 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1421}
1422
1424 const Array<OneD, const NekDouble> &inarray,
1425 Array<OneD, NekDouble> &outarray,
1426 const StdRegions::StdMatrixKey &mkey)
1427{
1428 StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
1429}
1430
1432 const Array<OneD, const NekDouble> &inarray,
1433 Array<OneD, NekDouble> &outarray, const StdRegions::StdMatrixKey &mkey)
1434{
1435 StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray, outarray, mkey);
1436}
1437
1439 const Array<OneD, const NekDouble> &inarray,
1440 Array<OneD, NekDouble> &outarray, const StdRegions::StdMatrixKey &mkey)
1441{
1442 StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray, outarray, mkey);
1443}
1444
1446 Array<OneD, NekDouble> &outarray,
1447 const StdRegions::StdMatrixKey &mkey)
1448{
1449 QuadExp::HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
1450}
1451
1453 const Array<OneD, const NekDouble> &inarray,
1454 Array<OneD, NekDouble> &outarray)
1455{
1456 int n_coeffs = inarray.size();
1457
1458 Array<OneD, NekDouble> coeff(n_coeffs);
1459 Array<OneD, NekDouble> coeff_tmp(n_coeffs, 0.0);
1460 Array<OneD, NekDouble> tmp, tmp2;
1461
1462 int nmodes0 = m_base[0]->GetNumModes();
1463 int nmodes1 = m_base[1]->GetNumModes();
1464 int numMax = nmodes0;
1465
1466 Vmath::Vcopy(n_coeffs, inarray, 1, coeff_tmp, 1);
1467
1468 const LibUtilities::PointsKey Pkey0(nmodes0,
1470 const LibUtilities::PointsKey Pkey1(nmodes1,
1472 LibUtilities::BasisKey b0(m_base[0]->GetBasisType(), nmodes0, Pkey0);
1473 LibUtilities::BasisKey b1(m_base[1]->GetBasisType(), nmodes1, Pkey1);
1474 LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A, nmodes0, Pkey0);
1475 LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_A, nmodes1, Pkey1);
1476
1477 LibUtilities::InterpCoeff2D(b0, b1, coeff_tmp, bortho0, bortho1, coeff);
1478
1479 Vmath::Zero(n_coeffs, coeff_tmp, 1);
1480
1481 int cnt = 0;
1482 for (int i = 0; i < numMin + 1; ++i)
1483 {
1484 Vmath::Vcopy(numMin, tmp = coeff + cnt, 1, tmp2 = coeff_tmp + cnt, 1);
1485
1486 cnt = i * numMax;
1487 }
1488
1489 LibUtilities::InterpCoeff2D(bortho0, bortho1, coeff_tmp, b0, b1, outarray);
1490}
1491
1493 const Array<OneD, const NekDouble> &inarray,
1495{
1496 if (m_metrics.count(eMetricLaplacian00) == 0)
1497 {
1499 }
1500
1501 int nquad0 = m_base[0]->GetNumPoints();
1502 int nquad1 = m_base[1]->GetNumPoints();
1503 int nqtot = nquad0 * nquad1;
1504 int nmodes0 = m_base[0]->GetNumModes();
1505 int nmodes1 = m_base[1]->GetNumModes();
1506 int wspsize =
1507 max(max(max(nqtot, m_ncoeffs), nquad1 * nmodes0), nquad0 * nmodes1);
1508
1509 ASSERTL1(wsp.size() >= 3 * wspsize, "Workspace is of insufficient size.");
1510
1511 const Array<OneD, const NekDouble> &base0 = m_base[0]->GetBdata();
1512 const Array<OneD, const NekDouble> &base1 = m_base[1]->GetBdata();
1513 const Array<OneD, const NekDouble> &dbase0 = m_base[0]->GetDbdata();
1514 const Array<OneD, const NekDouble> &dbase1 = m_base[1]->GetDbdata();
1515 const Array<OneD, const NekDouble> &metric00 =
1517 const Array<OneD, const NekDouble> &metric01 =
1519 const Array<OneD, const NekDouble> &metric11 =
1521
1522 // Allocate temporary storage
1523 Array<OneD, NekDouble> wsp0(wsp);
1524 Array<OneD, NekDouble> wsp1(wsp + wspsize);
1525 Array<OneD, NekDouble> wsp2(wsp + 2 * wspsize);
1526
1527 StdExpansion2D::PhysTensorDeriv(inarray, wsp1, wsp2);
1528
1529 // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1530 // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1531 // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1532 // especially for this purpose
1533 Vmath::Vvtvvtp(nqtot, &metric00[0], 1, &wsp1[0], 1, &metric01[0], 1,
1534 &wsp2[0], 1, &wsp0[0], 1);
1535 Vmath::Vvtvvtp(nqtot, &metric01[0], 1, &wsp1[0], 1, &metric11[0], 1,
1536 &wsp2[0], 1, &wsp2[0], 1);
1537
1538 // outarray = m = (D_xi1 * B)^T * k
1539 // wsp1 = n = (D_xi2 * B)^T * l
1540 IProductWRTBase_SumFacKernel(dbase0, base1, wsp0, outarray, wsp1, false,
1541 true);
1542 IProductWRTBase_SumFacKernel(base0, dbase1, wsp2, wsp1, wsp0, true, false);
1543
1544 // outarray = outarray + wsp1
1545 // = L * u_hat
1546 Vmath::Vadd(m_ncoeffs, wsp1.get(), 1, outarray.get(), 1, outarray.get(), 1);
1547}
1548
1550{
1551 if (m_metrics.count(eMetricQuadrature) == 0)
1552 {
1554 }
1555
1556 const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
1557 const unsigned int nqtot = GetTotPoints();
1558 const unsigned int dim = 2;
1559 const MetricType m[3][3] = {
1563
1564 const Array<TwoD, const NekDouble> gmat =
1565 m_metricinfo->GetGmat(GetPointsKeys());
1566 for (unsigned int i = 0; i < dim; ++i)
1567 {
1568 for (unsigned int j = i; j < dim; ++j)
1569 {
1570 m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1571 if (type == SpatialDomains::eDeformed)
1572 {
1573 Vmath::Vcopy(nqtot, &gmat[i * dim + j][0], 1,
1574 &m_metrics[m[i][j]][0], 1);
1575 }
1576 else
1577 {
1578 Vmath::Fill(nqtot, gmat[i * dim + j][0], &m_metrics[m[i][j]][0],
1579 1);
1580 }
1582 }
1583 }
1584}
1585
1587 const StdRegions::StdMatrixKey &mkey)
1588{
1589 int nq = GetTotPoints();
1590
1591 // Calculate sqrt of the Jacobian
1593 Array<OneD, NekDouble> sqrt_jac(nq);
1594 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1595 {
1596 Vmath::Vsqrt(nq, jac, 1, sqrt_jac, 1);
1597 }
1598 else
1599 {
1600 Vmath::Fill(nq, sqrt(jac[0]), sqrt_jac, 1);
1601 }
1602
1603 // Multiply array by sqrt(Jac)
1604 Vmath::Vmul(nq, sqrt_jac, 1, array, 1, array, 1);
1605
1606 // Apply std region filter
1607 StdQuadExp::v_SVVLaplacianFilter(array, mkey);
1608
1609 // Divide by sqrt(Jac)
1610 Vmath::Vdiv(nq, array, 1, sqrt_jac, 1, array, 1);
1611}
1612
1613/** @brief: This method gets all of the factors which are
1614 required as part of the Gradient Jump Penalty (GJP)
1615 stabilisation and involves the product of the normal and
1616 geometric factors along the element trace.
1617*/
1619 Array<OneD, Array<OneD, NekDouble>> &d0factors,
1620 Array<OneD, Array<OneD, NekDouble>> &d1factors,
1621 [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &d2factors)
1622{
1623 int nquad0 = GetNumPoints(0);
1624 int nquad1 = GetNumPoints(1);
1625
1627 m_metricinfo->GetDerivFactors(GetPointsKeys());
1628
1629 if (d0factors.size() != 4)
1630 {
1631 d0factors = Array<OneD, Array<OneD, NekDouble>>(4);
1632 d1factors = Array<OneD, Array<OneD, NekDouble>>(4);
1633 }
1634
1635 if (d0factors[0].size() != nquad0)
1636 {
1637 d0factors[0] = Array<OneD, NekDouble>(nquad0);
1638 d0factors[2] = Array<OneD, NekDouble>(nquad0);
1639 d1factors[0] = Array<OneD, NekDouble>(nquad0);
1640 d1factors[2] = Array<OneD, NekDouble>(nquad0);
1641 }
1642
1643 if (d0factors[1].size() != nquad1)
1644 {
1645 d0factors[1] = Array<OneD, NekDouble>(nquad1);
1646 d0factors[3] = Array<OneD, NekDouble>(nquad1);
1647 d1factors[1] = Array<OneD, NekDouble>(nquad1);
1648 d1factors[3] = Array<OneD, NekDouble>(nquad1);
1649 }
1650
1651 // Outwards normals
1653 GetTraceNormal(0);
1655 GetTraceNormal(1);
1657 GetTraceNormal(2);
1659 GetTraceNormal(3);
1660
1661 int ncoords = normal_0.size();
1662
1663 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1664 {
1665 // needs checking for 3D coords
1666
1667 // factors 0 and 2
1668 for (int i = 0; i < nquad0; ++i)
1669 {
1670 d0factors[0][i] = df[0][i] * normal_0[0][i];
1671 d0factors[2][i] = df[0][nquad0 * (nquad1 - 1) + i] * normal_2[0][i];
1672
1673 d1factors[0][i] = df[1][i] * normal_0[0][i];
1674 d1factors[2][i] = df[1][nquad0 * (nquad1 - 1) + i] * normal_2[0][i];
1675 }
1676
1677 for (int n = 1; n < ncoords; ++n)
1678 {
1679 // d xi_1/dy n_y
1680 // needs checking for 3D coords
1681 for (int i = 0; i < nquad0; ++i)
1682 {
1683 d0factors[0][i] += df[2 * n][i] * normal_0[n][i];
1684 d0factors[2][i] +=
1685 df[2 * n][nquad0 * (nquad1 - 1) + i] * normal_2[n][i];
1686
1687 d1factors[0][i] += df[2 * n + 1][i] * normal_0[n][i];
1688 d1factors[2][i] +=
1689 df[2 * n + 1][nquad0 * (nquad1 - 1) + i] * normal_2[n][i];
1690 }
1691 }
1692
1693 // faces 1 and 3
1694 for (int i = 0; i < nquad1; ++i)
1695 {
1696 d0factors[1][i] = df[0][(i + 1) * nquad0 - 1] * normal_1[0][i];
1697 d0factors[3][i] = df[0][i * nquad0] * normal_3[0][i];
1698
1699 d1factors[1][i] = df[1][(i + 1) * nquad0 - 1] * normal_1[0][i];
1700 d1factors[3][i] = df[1][i * nquad0] * normal_3[0][i];
1701 }
1702
1703 for (int n = 1; n < ncoords; ++n)
1704 {
1705 for (int i = 0; i < nquad1; ++i)
1706 {
1707 d0factors[1][i] +=
1708 df[2 * n][(i + 1) * nquad0 - 1] * normal_1[n][i];
1709 d0factors[3][i] += df[2 * n][i * nquad0] * normal_3[n][i];
1710
1711 d1factors[1][i] +=
1712 df[2 * n + 1][(i + 1) * nquad0 - 1] * normal_1[n][i];
1713 d1factors[3][i] += df[2 * n + 1][i * nquad0] * normal_3[n][i];
1714 }
1715 }
1716 }
1717 else
1718 {
1719 // d xi_2/dx n_x
1720 for (int i = 0; i < nquad0; ++i)
1721 {
1722 d1factors[0][i] = df[1][0] * normal_0[0][i];
1723 d1factors[2][i] = df[1][0] * normal_2[0][i];
1724 }
1725
1726 // d xi_1/dx n_x
1727 for (int i = 0; i < nquad1; ++i)
1728 {
1729 d0factors[1][i] = df[0][0] * normal_1[0][i];
1730 d0factors[3][i] = df[0][0] * normal_3[0][i];
1731 }
1732
1733 for (int n = 1; n < ncoords; ++n)
1734 {
1735 // d xi_2/dy n_y
1736 // needs checking for 3D coords
1737 for (int i = 0; i < nquad0; ++i)
1738 {
1739 d1factors[0][i] += df[2 * n + 1][0] * normal_0[n][i];
1740 d1factors[2][i] += df[2 * n + 1][0] * normal_2[n][i];
1741 }
1742
1743 // d xi_1/dy n_y
1744 // needs checking for 3D coords
1745 for (int i = 0; i < nquad1; ++i)
1746 {
1747 d0factors[1][i] += df[2 * n][0] * normal_1[n][i];
1748 d0factors[3][i] += df[2 * n][0] * normal_3[n][i];
1749 }
1750 }
1751
1752 // d1factors
1753 // d xi_1/dx n_x
1754 for (int i = 0; i < nquad0; ++i)
1755 {
1756 d0factors[0][i] = df[0][0] * normal_0[0][i];
1757 d0factors[2][i] = df[0][0] * normal_2[0][i];
1758 }
1759
1760 // d xi_2/dx n_x
1761 for (int i = 0; i < nquad1; ++i)
1762 {
1763 d1factors[1][i] = df[1][0] * normal_1[0][i];
1764 d1factors[3][i] = df[1][0] * normal_3[0][i];
1765 }
1766
1767 for (int n = 1; n < ncoords; ++n)
1768 {
1769 // d xi_1/dy n_y
1770 // needs checking for 3D coords
1771 for (int i = 0; i < nquad0; ++i)
1772 {
1773 d0factors[0][i] += df[2 * n][0] * normal_0[n][i];
1774 d0factors[2][i] += df[2 * n][0] * normal_2[n][i];
1775 }
1776
1777 // d xi_2/dy n_y
1778 // needs checking for 3D coords
1779 for (int i = 0; i < nquad1; ++i)
1780 {
1781 d1factors[1][i] += df[2 * n + 1][0] * normal_1[n][i];
1782 d1factors[3][i] += df[2 * n + 1][0] * normal_3[n][i];
1783 }
1784 }
1785 }
1786}
1787} // 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 v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
Definition: Expansion.cpp:530
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_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1401
StdRegions::Orientation v_GetTraceOrient(int edge) override
Definition: QuadExp.cpp:1346
DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey) override
Definition: QuadExp.cpp:1391
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: QuadExp.cpp:227
DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey) override
Definition: QuadExp.cpp:1381
StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const override
Definition: QuadExp.cpp:527
void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1445
void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords) override
Definition: QuadExp.cpp:545
void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculate the inner product of inarray with respect to the basis B=base0*base1 and put into outarray.
Definition: QuadExp.cpp:370
void v_ComputeLaplacianMetric() override
Definition: QuadExp.cpp:1549
void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp) override
Definition: QuadExp.cpp:1492
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
Definition: QuadExp.cpp:390
void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray) override
Definition: QuadExp.cpp:452
void GetEdgeInterpVals(const int edge, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: QuadExp.cpp:657
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: QuadExp.cpp:418
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: QuadExp.cpp:486
void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1408
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: QuadExp.cpp:574
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: QuadExp.cpp:595
void v_WeakDirectionalDerivMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1431
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: QuadExp.cpp:92
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray) override
Integrates the specified function over the domain.
Definition: QuadExp.cpp:69
void v_DropLocMatrix(const MatrixKey &mkey) override
Definition: QuadExp.cpp:1386
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: QuadExp.h:244
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: QuadExp.cpp:1256
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
Definition: QuadExp.cpp:538
void v_ComputeTraceNormal(const int edge) override
Definition: QuadExp.cpp:977
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: QuadExp.cpp:181
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: QuadExp.h:246
StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const override
Definition: QuadExp.cpp:521
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1351
void v_DropLocStaticCondMatrix(const MatrixKey &mkey) override
Definition: QuadExp.cpp:1396
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: QuadExp.cpp:383
void v_MassLevelCurvatureMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1438
DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1371
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1586
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 (GJP) s...
Definition: QuadExp.cpp:1618
NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals) override
Definition: QuadExp.cpp:566
void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: QuadExp.cpp:250
void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: QuadExp.cpp:1452
void v_GetTracePhysMap(const int edge, Array< OneD, int > &outarray) override
Definition: QuadExp.cpp:721
void v_GetTraceQFactors(const int edge, Array< OneD, NekDouble > &outarray) override
Definition: QuadExp.cpp:763
void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1423
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 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
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:752
const LibUtilities::BasisKey GetTraceBasisKey(const int i, int k=-1) const
This function returns the basis key belonging to the i-th trace.
Definition: StdExpansion.h:299
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:613
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:528
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:684
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:367
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:674
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:424
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:205
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space.
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:218
void 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:723
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 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 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
void InterpCoeff2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
Definition: InterpCoeff.cpp:67
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
Definition: Interp.cpp:101
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:231
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:46
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
Definition: BasisType.h:57
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:42
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:56
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
std::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:248
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< StdQuadExp > StdQuadExpSharedPtr
Definition: StdQuadExp.h:224
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:402
StdRegions::ConstFactorMap factors
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
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294