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.data(), outarray.data() + 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]).data(),
299 (physEdge[i]).data() + 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.data(), 1, 0.0, result.data(), 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.data(), 1, tmp1.data(),
476 1);
477 Vmath::Vmul(nqtot, &df[2 * dir + 1][0], 1, inarray.data(), 1,
478 tmp2.data(), 1);
479 }
480 else
481 {
482 Vmath::Smul(nqtot, df[2 * dir][0], inarray.data(), 1, tmp1.data(), 1);
483 Vmath::Smul(nqtot, df[2 * dir + 1][0], inarray.data(), 1, tmp2.data(),
484 1);
485 }
486}
487
492{
493 int nq = m_base[0]->GetNumPoints() * m_base[1]->GetNumPoints();
495
497 GetLeftAdjacentElementExp()->GetTraceNormal(
499
500 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
501 {
502 Vmath::Vvtvvtp(nq, &normals[0][0], 1, &Fx[0], 1, &normals[1][0], 1,
503 &Fy[0], 1, &Fn[0], 1);
504 Vmath::Vvtvp(nq, &normals[2][0], 1, &Fz[0], 1, &Fn[0], 1, &Fn[0], 1);
505 }
506 else
507 {
508 Vmath::Svtsvtp(nq, normals[0][0], &Fx[0], 1, normals[1][0], &Fy[0], 1,
509 &Fn[0], 1);
510 Vmath::Svtvp(nq, normals[2][0], &Fz[0], 1, &Fn[0], 1, &Fn[0], 1);
511 }
512
513 IProductWRTBase(Fn, outarray);
514}
515
517 const Array<OneD, const Array<OneD, NekDouble>> &Fvec,
518 Array<OneD, NekDouble> &outarray)
519{
520 NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
521}
522
524{
526 m_base[0]->GetBasisKey(), m_base[1]->GetBasisKey());
527}
528
530{
532 m_base[0]->GetPointsKey());
534 m_base[1]->GetPointsKey());
535
537 bkey1);
538}
539
541 Array<OneD, NekDouble> &coords_1,
542 Array<OneD, NekDouble> &coords_2)
543{
544 Expansion::v_GetCoords(coords_0, coords_1, coords_2);
545}
546
549{
550 int i;
551
552 ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 && Lcoords[1] >= -1.0 &&
553 Lcoords[1] <= 1.0,
554 "Local coordinates are not in region [-1,1]");
555
556 m_geom->FillGeom();
557 for (i = 0; i < m_geom->GetCoordim(); ++i)
558 {
559 coords[i] = m_geom->GetCoord(i, Lcoords);
560 }
561}
562
563/**
564 * Given the local cartesian coordinate \a Lcoord evaluate the
565 * value of physvals at this point by calling through to the
566 * StdExpansion method
567 */
569 const Array<OneD, const NekDouble> &Lcoord,
570 const Array<OneD, const NekDouble> &physvals)
571{
572 // Evaluate point in local (eta) coordinates.
573 return StdExpansion2D::v_PhysEvaluate(Lcoord, physvals);
574}
575
577 const Array<OneD, const NekDouble> &physvals)
578{
580
581 ASSERTL0(m_geom, "m_geom not defined");
582 m_geom->GetLocCoords(coord, Lcoord);
583
584 return StdExpansion2D::v_PhysEvaluate(Lcoord, physvals);
585}
586
588 const Array<OneD, NekDouble> &coord,
589 const Array<OneD, const NekDouble> &inarray,
590 std::array<NekDouble, 3> &firstOrderDerivs)
591{
592 Array<OneD, NekDouble> Lcoord(2);
593 ASSERTL0(m_geom, "m_geom not defined");
594 m_geom->GetLocCoords(coord, Lcoord);
595 return StdQuadExp::v_PhysEvalFirstDeriv(Lcoord, inarray, firstOrderDerivs);
596}
597
599 const int edge, const StdRegions::StdExpansionSharedPtr &EdgeExp,
600 const Array<OneD, const NekDouble> &inarray,
602{
603 int nquad0 = m_base[0]->GetNumPoints();
604 int nquad1 = m_base[1]->GetNumPoints();
605
606 // Implementation for all the basis except Gauss points
609 {
610 switch (edge)
611 {
612 case 0:
613 Vmath::Vcopy(nquad0, &(inarray[0]), 1, &(outarray[0]), 1);
614 break;
615 case 1:
616 Vmath::Vcopy(nquad1, &(inarray[0]) + (nquad0 - 1), nquad0,
617 &(outarray[0]), 1);
618 break;
619 case 2:
620 Vmath::Vcopy(nquad0, &(inarray[0]) + nquad0 * (nquad1 - 1), 1,
621 &(outarray[0]), 1);
622 break;
623 case 3:
624 Vmath::Vcopy(nquad1, &(inarray[0]), nquad0, &(outarray[0]), 1);
625 break;
626 default:
627 ASSERTL0(false, "edge value (< 3) is out of range");
628 break;
629 }
630 }
631 else
632 {
633 QuadExp::GetEdgeInterpVals(edge, inarray, outarray);
634 }
635
636 // Interpolate if required
637 if (m_base[edge % 2]->GetPointsKey() !=
638 EdgeExp->GetBasis(0)->GetPointsKey())
639 {
640 Array<OneD, NekDouble> outtmp(max(nquad0, nquad1));
641
642 outtmp = outarray;
643
644 LibUtilities::Interp1D(m_base[edge % 2]->GetPointsKey(), outtmp,
645 EdgeExp->GetBasis(0)->GetPointsKey(), outarray);
646 }
647
648 if (orient == StdRegions::eNoOrientation)
649 {
650 orient = GetTraceOrient(edge);
651 }
652 // Reverse data if necessary
653 if (orient == StdRegions::eBackwards)
654 {
655 Vmath::Reverse(EdgeExp->GetNumPoints(0), &outarray[0], 1, &outarray[0],
656 1);
657 }
658}
659
660void QuadExp::GetEdgeInterpVals(const int edge,
661 const Array<OneD, const NekDouble> &inarray,
662 Array<OneD, NekDouble> &outarray)
663{
664 int i;
665 int nq0 = m_base[0]->GetNumPoints();
666 int nq1 = m_base[1]->GetNumPoints();
667
670
672 *this, factors);
673
674 DNekScalMatSharedPtr mat_gauss = m_matrixManager[key];
675
676 switch (edge)
677 {
678 case 0:
679 {
680 for (i = 0; i < nq0; i++)
681 {
682 outarray[i] = Blas::Ddot(
683 nq1, mat_gauss->GetOwnedMatrix()->GetPtr().data(), 1,
684 &inarray[i], nq0);
685 }
686 break;
687 }
688 case 1:
689 {
690 for (i = 0; i < nq1; i++)
691 {
692 outarray[i] = Blas::Ddot(
693 nq0, mat_gauss->GetOwnedMatrix()->GetPtr().data(), 1,
694 &inarray[i * nq0], 1);
695 }
696 break;
697 }
698 case 2:
699 {
700 for (i = 0; i < nq0; i++)
701 {
702 outarray[i] = Blas::Ddot(
703 nq1, mat_gauss->GetOwnedMatrix()->GetPtr().data(), 1,
704 &inarray[i], nq0);
705 }
706 break;
707 }
708 case 3:
709 {
710 for (i = 0; i < nq1; i++)
711 {
712 outarray[i] = Blas::Ddot(
713 nq0, mat_gauss->GetOwnedMatrix()->GetPtr().data(), 1,
714 &inarray[i * nq0], 1);
715 }
716 break;
717 }
718 default:
719 ASSERTL0(false, "edge value (< 3) is out of range");
720 break;
721 }
722}
723
724void QuadExp::v_GetTracePhysMap(const int edge, Array<OneD, int> &outarray)
725{
726 int nquad0 = m_base[0]->GetNumPoints();
727 int nquad1 = m_base[1]->GetNumPoints();
728
729 // Get points in Cartesian orientation
730 switch (edge)
731 {
732 case 0:
733 outarray = Array<OneD, int>(nquad0);
734 for (int i = 0; i < nquad0; ++i)
735 {
736 outarray[i] = i;
737 }
738 break;
739 case 1:
740 outarray = Array<OneD, int>(nquad1);
741 for (int i = 0; i < nquad1; ++i)
742 {
743 outarray[i] = (nquad0 - 1) + i * nquad0;
744 }
745 break;
746 case 2:
747 outarray = Array<OneD, int>(nquad0);
748 for (int i = 0; i < nquad0; ++i)
749 {
750 outarray[i] = i + nquad0 * (nquad1 - 1);
751 }
752 break;
753 case 3:
754 outarray = Array<OneD, int>(nquad1);
755 for (int i = 0; i < nquad1; ++i)
756 {
757 outarray[i] = i * nquad0;
758 }
759 break;
760 default:
761 ASSERTL0(false, "edge value (< 3) is out of range");
762 break;
763 }
764}
765
766void QuadExp::v_GetTraceQFactors(const int edge,
767 Array<OneD, NekDouble> &outarray)
768{
769 int i;
770 int nquad0 = m_base[0]->GetNumPoints();
771 int nquad1 = m_base[1]->GetNumPoints();
772
774 const Array<OneD, const NekDouble> &jac = m_metricinfo->GetJac(ptsKeys);
776 m_metricinfo->GetDerivFactors(ptsKeys);
777
778 Array<OneD, NekDouble> j(max(nquad0, nquad1), 0.0);
779 Array<OneD, NekDouble> g0(max(nquad0, nquad1), 0.0);
780 Array<OneD, NekDouble> g1(max(nquad0, nquad1), 0.0);
781 Array<OneD, NekDouble> g2(max(nquad0, nquad1), 0.0);
782 Array<OneD, NekDouble> g3(max(nquad0, nquad1), 0.0);
783
784 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
785 {
786 // Implementation for all the basis except Gauss points
789 {
790 switch (edge)
791 {
792 case 0:
793 Vmath::Vcopy(nquad0, &(df[1][0]), 1, &(g1[0]), 1);
794 Vmath::Vcopy(nquad0, &(df[3][0]), 1, &(g3[0]), 1);
795 Vmath::Vcopy(nquad0, &(jac[0]), 1, &(j[0]), 1);
796
797 for (i = 0; i < nquad0; ++i)
798 {
799 outarray[i] =
800 j[i] * sqrt(g1[i] * g1[i] + g3[i] * g3[i]);
801 }
802 break;
803 case 1:
804 Vmath::Vcopy(nquad1, &(df[0][0]) + (nquad0 - 1), nquad0,
805 &(g0[0]), 1);
806
807 Vmath::Vcopy(nquad1, &(df[2][0]) + (nquad0 - 1), nquad0,
808 &(g2[0]), 1);
809
810 Vmath::Vcopy(nquad1, &(jac[0]) + (nquad0 - 1), nquad0,
811 &(j[0]), 1);
812
813 for (i = 0; i < nquad1; ++i)
814 {
815 outarray[i] =
816 j[i] * sqrt(g0[i] * g0[i] + g2[i] * g2[i]);
817 }
818 break;
819 case 2:
820
821 Vmath::Vcopy(nquad0, &(df[1][0]) + (nquad0 * (nquad1 - 1)),
822 1, &(g1[0]), 1);
823
824 Vmath::Vcopy(nquad0, &(df[3][0]) + (nquad0 * (nquad1 - 1)),
825 1, &(g3[0]), 1);
826
827 Vmath::Vcopy(nquad0, &(jac[0]) + (nquad0 * (nquad1 - 1)), 1,
828 &(j[0]), 1);
829
830 for (i = 0; i < nquad0; ++i)
831 {
832 outarray[i] =
833 j[i] * sqrt(g1[i] * g1[i] + g3[i] * g3[i]);
834 }
835 break;
836 case 3:
837
838 Vmath::Vcopy(nquad1, &(df[0][0]), nquad0, &(g0[0]), 1);
839 Vmath::Vcopy(nquad1, &(df[2][0]), nquad0, &(g2[0]), 1);
840 Vmath::Vcopy(nquad1, &(jac[0]), nquad0, &(j[0]), 1);
841
842 for (i = 0; i < nquad1; ++i)
843 {
844 outarray[i] =
845 j[i] * sqrt(g0[i] * g0[i] + g2[i] * g2[i]);
846 }
847 break;
848 default:
849 ASSERTL0(false, "edge value (< 3) is out of range");
850 break;
851 }
852 }
853 else
854 {
855 int nqtot = nquad0 * nquad1;
856 Array<OneD, NekDouble> tmp_gmat0(nqtot, 0.0);
857 Array<OneD, NekDouble> tmp_gmat1(nqtot, 0.0);
858 Array<OneD, NekDouble> tmp_gmat2(nqtot, 0.0);
859 Array<OneD, NekDouble> tmp_gmat3(nqtot, 0.0);
860 Array<OneD, NekDouble> g0_edge(max(nquad0, nquad1), 0.0);
861 Array<OneD, NekDouble> g1_edge(max(nquad0, nquad1), 0.0);
862 Array<OneD, NekDouble> g2_edge(max(nquad0, nquad1), 0.0);
863 Array<OneD, NekDouble> g3_edge(max(nquad0, nquad1), 0.0);
864 Array<OneD, NekDouble> jac_edge(max(nquad0, nquad1), 0.0);
865
866 switch (edge)
867 {
868 case 0:
869 Vmath::Vmul(nqtot, &(df[1][0]), 1, &jac[0], 1,
870 &(tmp_gmat1[0]), 1);
871 Vmath::Vmul(nqtot, &(df[3][0]), 1, &jac[0], 1,
872 &(tmp_gmat3[0]), 1);
873 QuadExp::GetEdgeInterpVals(edge, tmp_gmat1, g1_edge);
874 QuadExp::GetEdgeInterpVals(edge, tmp_gmat3, g3_edge);
875
876 for (i = 0; i < nquad0; ++i)
877 {
878 outarray[i] = sqrt(g1_edge[i] * g1_edge[i] +
879 g3_edge[i] * g3_edge[i]);
880 }
881 break;
882
883 case 1:
884 Vmath::Vmul(nqtot, &(df[0][0]), 1, &jac[0], 1,
885 &(tmp_gmat0[0]), 1);
886 Vmath::Vmul(nqtot, &(df[2][0]), 1, &jac[0], 1,
887 &(tmp_gmat2[0]), 1);
888 QuadExp::GetEdgeInterpVals(edge, tmp_gmat0, g0_edge);
889 QuadExp::GetEdgeInterpVals(edge, tmp_gmat2, g2_edge);
890
891 for (i = 0; i < nquad1; ++i)
892 {
893 outarray[i] = sqrt(g0_edge[i] * g0_edge[i] +
894 g2_edge[i] * g2_edge[i]);
895 }
896
897 break;
898 case 2:
899
900 Vmath::Vmul(nqtot, &(df[1][0]), 1, &jac[0], 1,
901 &(tmp_gmat1[0]), 1);
902 Vmath::Vmul(nqtot, &(df[3][0]), 1, &jac[0], 1,
903 &(tmp_gmat3[0]), 1);
904 QuadExp::GetEdgeInterpVals(edge, tmp_gmat1, g1_edge);
905 QuadExp::GetEdgeInterpVals(edge, tmp_gmat3, g3_edge);
906
907 for (i = 0; i < nquad0; ++i)
908 {
909 outarray[i] = sqrt(g1_edge[i] * g1_edge[i] +
910 g3_edge[i] * g3_edge[i]);
911 }
912
913 Vmath::Reverse(nquad0, &outarray[0], 1, &outarray[0], 1);
914
915 break;
916 case 3:
917 Vmath::Vmul(nqtot, &(df[0][0]), 1, &jac[0], 1,
918 &(tmp_gmat0[0]), 1);
919 Vmath::Vmul(nqtot, &(df[2][0]), 1, &jac[0], 1,
920 &(tmp_gmat2[0]), 1);
921 QuadExp::GetEdgeInterpVals(edge, tmp_gmat0, g0_edge);
922 QuadExp::GetEdgeInterpVals(edge, tmp_gmat2, g2_edge);
923
924 for (i = 0; i < nquad1; ++i)
925 {
926 outarray[i] = sqrt(g0_edge[i] * g0_edge[i] +
927 g2_edge[i] * g2_edge[i]);
928 }
929
930 Vmath::Reverse(nquad1, &outarray[0], 1, &outarray[0], 1);
931
932 break;
933 default:
934 ASSERTL0(false, "edge value (< 3) is out of range");
935 break;
936 }
937 }
938 }
939 else
940 {
941
942 switch (edge)
943 {
944 case 0:
945
946 for (i = 0; i < nquad0; ++i)
947 {
948 outarray[i] = jac[0] * sqrt(df[1][0] * df[1][0] +
949 df[3][0] * df[3][0]);
950 }
951 break;
952 case 1:
953 for (i = 0; i < nquad1; ++i)
954 {
955 outarray[i] = jac[0] * sqrt(df[0][0] * df[0][0] +
956 df[2][0] * df[2][0]);
957 }
958 break;
959 case 2:
960 for (i = 0; i < nquad0; ++i)
961 {
962 outarray[i] = jac[0] * sqrt(df[1][0] * df[1][0] +
963 df[3][0] * df[3][0]);
964 }
965 break;
966 case 3:
967 for (i = 0; i < nquad1; ++i)
968 {
969 outarray[i] = jac[0] * sqrt(df[0][0] * df[0][0] +
970 df[2][0] * df[2][0]);
971 }
972 break;
973 default:
974 ASSERTL0(false, "edge value (< 3) is out of range");
975 break;
976 }
977 }
978}
979
981{
982 int i;
983 const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
984 GetGeom()->GetMetricInfo();
985 SpatialDomains::GeomType type = geomFactors->GetGtype();
986
988 for (i = 0; i < ptsKeys.size(); ++i)
989 {
990 // Need at least 2 points for computing normals
991 if (ptsKeys[i].GetNumPoints() == 1)
992 {
993 LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
994 ptsKeys[i] = pKey;
995 }
996 }
997
999 geomFactors->GetDerivFactors(ptsKeys);
1000 const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
1001
1002 // The points of normals should follow trace basis, not local basis.
1004
1005 int nqe = tobasis.GetNumPoints();
1006 int vCoordDim = GetCoordim();
1007
1010 for (i = 0; i < vCoordDim; ++i)
1011 {
1012 normal[i] = Array<OneD, NekDouble>(nqe);
1013 }
1014
1015 size_t nqb = nqe;
1016 size_t nbnd = edge;
1019
1020 // Regular geometry case
1021 if ((type == SpatialDomains::eRegular) ||
1023 {
1024 NekDouble fac;
1025 // Set up normals
1026 switch (edge)
1027 {
1028 case 0:
1029 for (i = 0; i < vCoordDim; ++i)
1030 {
1031 Vmath::Fill(nqe, -df[2 * i + 1][0], normal[i], 1);
1032 }
1033 break;
1034 case 1:
1035 for (i = 0; i < vCoordDim; ++i)
1036 {
1037 Vmath::Fill(nqe, df[2 * i][0], normal[i], 1);
1038 }
1039 break;
1040 case 2:
1041 for (i = 0; i < vCoordDim; ++i)
1042 {
1043 Vmath::Fill(nqe, df[2 * i + 1][0], normal[i], 1);
1044 }
1045 break;
1046 case 3:
1047 for (i = 0; i < vCoordDim; ++i)
1048 {
1049 Vmath::Fill(nqe, -df[2 * i][0], normal[i], 1);
1050 }
1051 break;
1052 default:
1053 ASSERTL0(false, "edge is out of range (edge < 4)");
1054 }
1055
1056 // normalise
1057 fac = 0.0;
1058 for (i = 0; i < vCoordDim; ++i)
1059 {
1060 fac += normal[i][0] * normal[i][0];
1061 }
1062 fac = 1.0 / sqrt(fac);
1063
1064 Vmath::Fill(nqb, fac, length, 1);
1065
1066 for (i = 0; i < vCoordDim; ++i)
1067 {
1068 Vmath::Smul(nqe, fac, normal[i], 1, normal[i], 1);
1069 }
1070 }
1071 else // Set up deformed normals
1072 {
1073 int j;
1074
1075 int nquad0 = ptsKeys[0].GetNumPoints();
1076 int nquad1 = ptsKeys[1].GetNumPoints();
1077
1078 LibUtilities::PointsKey from_key;
1079
1080 Array<OneD, NekDouble> normals(vCoordDim * max(nquad0, nquad1), 0.0);
1081 Array<OneD, NekDouble> edgejac(vCoordDim * max(nquad0, nquad1), 0.0);
1082
1083 // Extract Jacobian along edges and recover local
1084 // derivates (dx/dr) for polynomial interpolation by
1085 // multiplying m_gmat by jacobian
1086
1087 // Implementation for all the basis except Gauss points
1090 {
1091 switch (edge)
1092 {
1093 case 0:
1094 for (j = 0; j < nquad0; ++j)
1095 {
1096 edgejac[j] = jac[j];
1097 for (i = 0; i < vCoordDim; ++i)
1098 {
1099 normals[i * nquad0 + j] =
1100 -df[2 * i + 1][j] * edgejac[j];
1101 }
1102 }
1103 from_key = ptsKeys[0];
1104 break;
1105 case 1:
1106 for (j = 0; j < nquad1; ++j)
1107 {
1108 edgejac[j] = jac[nquad0 * j + nquad0 - 1];
1109 for (i = 0; i < vCoordDim; ++i)
1110 {
1111 normals[i * nquad1 + j] =
1112 df[2 * i][nquad0 * j + nquad0 - 1] * edgejac[j];
1113 }
1114 }
1115 from_key = ptsKeys[1];
1116 break;
1117 case 2:
1118 for (j = 0; j < nquad0; ++j)
1119 {
1120 edgejac[j] = jac[nquad0 * (nquad1 - 1) + j];
1121 for (i = 0; i < vCoordDim; ++i)
1122 {
1123 normals[i * nquad0 + j] =
1124 (df[2 * i + 1][nquad0 * (nquad1 - 1) + j]) *
1125 edgejac[j];
1126 }
1127 }
1128 from_key = ptsKeys[0];
1129 break;
1130 case 3:
1131 for (j = 0; j < nquad1; ++j)
1132 {
1133 edgejac[j] = jac[nquad0 * j];
1134 for (i = 0; i < vCoordDim; ++i)
1135 {
1136 normals[i * nquad1 + j] =
1137 -df[2 * i][nquad0 * j] * edgejac[j];
1138 }
1139 }
1140 from_key = ptsKeys[1];
1141 break;
1142 default:
1143 ASSERTL0(false, "edge is out of range (edge < 3)");
1144 }
1145 }
1146 else
1147 {
1148 int nqtot = nquad0 * nquad1;
1149 Array<OneD, NekDouble> tmp_gmat(nqtot, 0.0);
1150 Array<OneD, NekDouble> tmp_gmat_edge(nqe, 0.0);
1151
1152 switch (edge)
1153 {
1154 case 0:
1155 for (j = 0; j < nquad0; ++j)
1156 {
1157 for (i = 0; i < vCoordDim; ++i)
1158 {
1159 Vmath::Vmul(nqtot, &(df[2 * i + 1][0]), 1, &jac[0],
1160 1, &(tmp_gmat[0]), 1);
1161 QuadExp::GetEdgeInterpVals(edge, tmp_gmat,
1162 tmp_gmat_edge);
1163 normals[i * nquad0 + j] = -tmp_gmat_edge[j];
1164 }
1165 }
1166 from_key = ptsKeys[0];
1167 break;
1168 case 1:
1169 for (j = 0; j < nquad1; ++j)
1170 {
1171 for (i = 0; i < vCoordDim; ++i)
1172 {
1173 Vmath::Vmul(nqtot, &(df[2 * i][0]), 1, &jac[0], 1,
1174 &(tmp_gmat[0]), 1);
1175 QuadExp::GetEdgeInterpVals(edge, tmp_gmat,
1176 tmp_gmat_edge);
1177 normals[i * nquad1 + j] = tmp_gmat_edge[j];
1178 }
1179 }
1180 from_key = ptsKeys[1];
1181 break;
1182 case 2:
1183 for (j = 0; j < nquad0; ++j)
1184 {
1185 for (i = 0; i < vCoordDim; ++i)
1186 {
1187 Vmath::Vmul(nqtot, &(df[2 * i + 1][0]), 1, &jac[0],
1188 1, &(tmp_gmat[0]), 1);
1189 QuadExp::GetEdgeInterpVals(edge, tmp_gmat,
1190 tmp_gmat_edge);
1191 normals[i * nquad0 + j] = tmp_gmat_edge[j];
1192 }
1193 }
1194 from_key = ptsKeys[0];
1195 break;
1196 case 3:
1197 for (j = 0; j < nquad1; ++j)
1198 {
1199 for (i = 0; i < vCoordDim; ++i)
1200 {
1201 Vmath::Vmul(nqtot, &(df[2 * i][0]), 1, &jac[0], 1,
1202 &(tmp_gmat[0]), 1);
1203 QuadExp::GetEdgeInterpVals(edge, tmp_gmat,
1204 tmp_gmat_edge);
1205 normals[i * nquad1 + j] = -tmp_gmat_edge[j];
1206 }
1207 }
1208 from_key = ptsKeys[1];
1209 break;
1210 default:
1211 ASSERTL0(false, "edge is out of range (edge < 3)");
1212 }
1213 }
1214
1215 int nq = from_key.GetNumPoints();
1216 Array<OneD, NekDouble> work(nqe, 0.0);
1217
1218 // interpolate Jacobian and invert
1219 LibUtilities::Interp1D(from_key, jac, tobasis.GetPointsKey(), work);
1220 Vmath::Sdiv(nqe, 1.0, &work[0], 1, &work[0], 1);
1221
1222 // interpolate
1223 for (i = 0; i < GetCoordim(); ++i)
1224 {
1225 LibUtilities::Interp1D(from_key, &normals[i * nq],
1226 tobasis.GetPointsKey(), &normal[i][0]);
1227 Vmath::Vmul(nqe, work, 1, normal[i], 1, normal[i], 1);
1228 }
1229
1230 // normalise normal vectors
1231 Vmath::Zero(nqe, work, 1);
1232 for (i = 0; i < GetCoordim(); ++i)
1233 {
1234 Vmath::Vvtvp(nqe, normal[i], 1, normal[i], 1, work, 1, work, 1);
1235 }
1236
1237 Vmath::Vsqrt(nqe, work, 1, work, 1);
1238 Vmath::Sdiv(nqe, 1.0, work, 1, work, 1);
1239
1240 Vmath::Vcopy(nqb, work, 1, length, 1);
1241
1242 for (i = 0; i < GetCoordim(); ++i)
1243 {
1244 Vmath::Vmul(nqe, normal[i], 1, work, 1, normal[i], 1);
1245 }
1246 }
1247 if (GetGeom()->GetEorient(edge) == StdRegions::eBackwards)
1248 {
1249 for (i = 0; i < vCoordDim; ++i)
1250 {
1251 if (geomFactors->GetGtype() == SpatialDomains::eDeformed)
1252 {
1253 Vmath::Reverse(nqe, normal[i], 1, normal[i], 1);
1254 }
1255 }
1256 }
1257}
1258
1260 const NekDouble *data, const std::vector<unsigned int> &nummodes,
1261 int mode_offset, NekDouble *coeffs,
1262 std::vector<LibUtilities::BasisType> &fromType)
1263{
1264 int data_order0 = nummodes[mode_offset];
1265 int fillorder0 = std::min(m_base[0]->GetNumModes(), data_order0);
1266
1267 int data_order1 = nummodes[mode_offset + 1];
1268 int order1 = m_base[1]->GetNumModes();
1269 int fillorder1 = min(order1, data_order1);
1270
1271 // Check if same basis
1272 if (fromType[0] != m_base[0]->GetBasisType() ||
1273 fromType[1] != m_base[1]->GetBasisType())
1274 {
1275 // Construct a quad with the appropriate basis type at our
1276 // quadrature points, and one more to do a forwards
1277 // transform. We can then copy the output to coeffs.
1278 StdRegions::StdQuadExp tmpQuad(
1279 LibUtilities::BasisKey(fromType[0], data_order0,
1280 m_base[0]->GetPointsKey()),
1281 LibUtilities::BasisKey(fromType[1], data_order1,
1282 m_base[1]->GetPointsKey()));
1283 StdRegions::StdQuadExp tmpQuad2(m_base[0]->GetBasisKey(),
1284 m_base[1]->GetBasisKey());
1285
1286 Array<OneD, const NekDouble> tmpData(tmpQuad.GetNcoeffs(), data);
1287 Array<OneD, NekDouble> tmpBwd(tmpQuad2.GetTotPoints());
1288 Array<OneD, NekDouble> tmpOut(tmpQuad2.GetNcoeffs());
1289
1290 tmpQuad.BwdTrans(tmpData, tmpBwd);
1291 tmpQuad2.FwdTrans(tmpBwd, tmpOut);
1292 Vmath::Vcopy(tmpOut.size(), &tmpOut[0], 1, coeffs, 1);
1293
1294 return;
1295 }
1296
1297 switch (m_base[0]->GetBasisType())
1298 {
1300 {
1301 int i;
1302 int cnt = 0;
1303 int cnt1 = 0;
1304
1306 "Extraction routine not set up for this basis");
1307
1308 Vmath::Zero(m_ncoeffs, coeffs, 1);
1309 for (i = 0; i < fillorder0; ++i)
1310 {
1311 Vmath::Vcopy(fillorder1, data + cnt, 1, coeffs + cnt1, 1);
1312 cnt += data_order1;
1313 cnt1 += order1;
1314 }
1315 }
1316 break;
1318 {
1319 LibUtilities::PointsKey p0(nummodes[0],
1321 LibUtilities::PointsKey p1(nummodes[1],
1323 LibUtilities::PointsKey t0(m_base[0]->GetNumModes(),
1325 LibUtilities::PointsKey t1(m_base[1]->GetNumModes(),
1327 LibUtilities::Interp2D(p0, p1, data, t0, t1, coeffs);
1328 }
1329 break;
1331 {
1332 // Assume that input is also Gll_Lagrange but no way to check;
1333 LibUtilities::PointsKey p0(nummodes[0],
1335 LibUtilities::PointsKey p1(nummodes[1],
1337 LibUtilities::PointsKey t0(m_base[0]->GetNumModes(),
1339 LibUtilities::PointsKey t1(m_base[1]->GetNumModes(),
1341 LibUtilities::Interp2D(p0, p1, data, t0, t1, coeffs);
1342 }
1343 break;
1344 default:
1345 ASSERTL0(false, "basis is either not set up or not hierarchicial");
1346 }
1347}
1348
1350{
1351 return m_geom->GetEorient(edge);
1352}
1353
1355{
1356 DNekMatSharedPtr returnval;
1357 switch (mkey.GetMatrixType())
1358 {
1366 returnval = Expansion2D::v_GenMatrix(mkey);
1367 break;
1368 default:
1369 returnval = StdQuadExp::v_GenMatrix(mkey);
1370 }
1371 return returnval;
1372}
1373
1375 const StdRegions::StdMatrixKey &mkey)
1376{
1377 LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1378 LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1381 return tmp->GetStdMatrix(mkey);
1382}
1383
1385{
1386 return m_matrixManager[mkey];
1387}
1388
1390{
1391 m_matrixManager.DeleteObject(mkey);
1392}
1393
1395{
1396 return m_staticCondMatrixManager[mkey];
1397}
1398
1400{
1401 m_staticCondMatrixManager.DeleteObject(mkey);
1402}
1403
1405 Array<OneD, NekDouble> &outarray,
1406 const StdRegions::StdMatrixKey &mkey)
1407{
1408 StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
1409}
1410
1412 Array<OneD, NekDouble> &outarray,
1413 const StdRegions::StdMatrixKey &mkey)
1414{
1415 QuadExp::LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
1416}
1417
1418void QuadExp::v_LaplacianMatrixOp(const int k1, const int k2,
1419 const Array<OneD, const NekDouble> &inarray,
1420 Array<OneD, NekDouble> &outarray,
1421 const StdRegions::StdMatrixKey &mkey)
1422{
1423 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1424}
1425
1427 const Array<OneD, const NekDouble> &inarray,
1428 Array<OneD, NekDouble> &outarray,
1429 const StdRegions::StdMatrixKey &mkey)
1430{
1431 StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
1432}
1433
1435 const Array<OneD, const NekDouble> &inarray,
1436 Array<OneD, NekDouble> &outarray, const StdRegions::StdMatrixKey &mkey)
1437{
1438 StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray, outarray, mkey);
1439}
1440
1442 const Array<OneD, const NekDouble> &inarray,
1443 Array<OneD, NekDouble> &outarray, const StdRegions::StdMatrixKey &mkey)
1444{
1445 StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray, outarray, mkey);
1446}
1447
1449 Array<OneD, NekDouble> &outarray,
1450 const StdRegions::StdMatrixKey &mkey)
1451{
1452 QuadExp::HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
1453}
1454
1456 const Array<OneD, const NekDouble> &inarray,
1457 Array<OneD, NekDouble> &outarray)
1458{
1459 int n_coeffs = inarray.size();
1460
1461 Array<OneD, NekDouble> coeff(n_coeffs);
1462 Array<OneD, NekDouble> coeff_tmp(n_coeffs, 0.0);
1463 Array<OneD, NekDouble> tmp, tmp2;
1464
1465 int nmodes0 = m_base[0]->GetNumModes();
1466 int nmodes1 = m_base[1]->GetNumModes();
1467 int numMax = nmodes0;
1468
1469 Vmath::Vcopy(n_coeffs, inarray, 1, coeff_tmp, 1);
1470
1471 const LibUtilities::PointsKey Pkey0(nmodes0,
1473 const LibUtilities::PointsKey Pkey1(nmodes1,
1475 LibUtilities::BasisKey b0(m_base[0]->GetBasisType(), nmodes0, Pkey0);
1476 LibUtilities::BasisKey b1(m_base[1]->GetBasisType(), nmodes1, Pkey1);
1477 LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A, nmodes0, Pkey0);
1478 LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_A, nmodes1, Pkey1);
1479
1480 LibUtilities::InterpCoeff2D(b0, b1, coeff_tmp, bortho0, bortho1, coeff);
1481
1482 Vmath::Zero(n_coeffs, coeff_tmp, 1);
1483
1484 int cnt = 0;
1485 for (int i = 0; i < numMin + 1; ++i)
1486 {
1487 Vmath::Vcopy(numMin, tmp = coeff + cnt, 1, tmp2 = coeff_tmp + cnt, 1);
1488
1489 cnt = i * numMax;
1490 }
1491
1492 LibUtilities::InterpCoeff2D(bortho0, bortho1, coeff_tmp, b0, b1, outarray);
1493}
1494
1496 const Array<OneD, const NekDouble> &inarray,
1498{
1499 if (m_metrics.count(eMetricLaplacian00) == 0)
1500 {
1502 }
1503
1504 int nquad0 = m_base[0]->GetNumPoints();
1505 int nquad1 = m_base[1]->GetNumPoints();
1506 int nqtot = nquad0 * nquad1;
1507 int nmodes0 = m_base[0]->GetNumModes();
1508 int nmodes1 = m_base[1]->GetNumModes();
1509 int wspsize =
1510 max(max(max(nqtot, m_ncoeffs), nquad1 * nmodes0), nquad0 * nmodes1);
1511
1512 ASSERTL1(wsp.size() >= 3 * wspsize, "Workspace is of insufficient size.");
1513
1514 const Array<OneD, const NekDouble> &base0 = m_base[0]->GetBdata();
1515 const Array<OneD, const NekDouble> &base1 = m_base[1]->GetBdata();
1516 const Array<OneD, const NekDouble> &dbase0 = m_base[0]->GetDbdata();
1517 const Array<OneD, const NekDouble> &dbase1 = m_base[1]->GetDbdata();
1518 const Array<OneD, const NekDouble> &metric00 =
1520 const Array<OneD, const NekDouble> &metric01 =
1522 const Array<OneD, const NekDouble> &metric11 =
1524
1525 // Allocate temporary storage
1526 Array<OneD, NekDouble> wsp0(wsp);
1527 Array<OneD, NekDouble> wsp1(wsp + wspsize);
1528 Array<OneD, NekDouble> wsp2(wsp + 2 * wspsize);
1529
1530 StdExpansion2D::PhysTensorDeriv(inarray, wsp1, wsp2);
1531
1532 // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1533 // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1534 // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1535 // especially for this purpose
1536 Vmath::Vvtvvtp(nqtot, &metric00[0], 1, &wsp1[0], 1, &metric01[0], 1,
1537 &wsp2[0], 1, &wsp0[0], 1);
1538 Vmath::Vvtvvtp(nqtot, &metric01[0], 1, &wsp1[0], 1, &metric11[0], 1,
1539 &wsp2[0], 1, &wsp2[0], 1);
1540
1541 // outarray = m = (D_xi1 * B)^T * k
1542 // wsp1 = n = (D_xi2 * B)^T * l
1543 IProductWRTBase_SumFacKernel(dbase0, base1, wsp0, outarray, wsp1, false,
1544 true);
1545 IProductWRTBase_SumFacKernel(base0, dbase1, wsp2, wsp1, wsp0, true, false);
1546
1547 // outarray = outarray + wsp1
1548 // = L * u_hat
1549 Vmath::Vadd(m_ncoeffs, wsp1.data(), 1, outarray.data(), 1, outarray.data(),
1550 1);
1551}
1552
1554{
1555 if (m_metrics.count(eMetricQuadrature) == 0)
1556 {
1558 }
1559
1560 const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
1561 const unsigned int nqtot = GetTotPoints();
1562 const unsigned int dim = 2;
1563 const MetricType m[3][3] = {
1567
1568 const Array<TwoD, const NekDouble> gmat =
1569 m_metricinfo->GetGmat(GetPointsKeys());
1570 for (unsigned int i = 0; i < dim; ++i)
1571 {
1572 for (unsigned int j = i; j < dim; ++j)
1573 {
1574 m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1575 if (type == SpatialDomains::eDeformed)
1576 {
1577 Vmath::Vcopy(nqtot, &gmat[i * dim + j][0], 1,
1578 &m_metrics[m[i][j]][0], 1);
1579 }
1580 else
1581 {
1582 Vmath::Fill(nqtot, gmat[i * dim + j][0], &m_metrics[m[i][j]][0],
1583 1);
1584 }
1586 }
1587 }
1588}
1589
1591 const StdRegions::StdMatrixKey &mkey)
1592{
1593 int nq = GetTotPoints();
1594
1595 // Calculate sqrt of the Jacobian
1597 Array<OneD, NekDouble> sqrt_jac(nq);
1598 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1599 {
1600 Vmath::Vsqrt(nq, jac, 1, sqrt_jac, 1);
1601 }
1602 else
1603 {
1604 Vmath::Fill(nq, sqrt(jac[0]), sqrt_jac, 1);
1605 }
1606
1607 // Multiply array by sqrt(Jac)
1608 Vmath::Vmul(nq, sqrt_jac, 1, array, 1, array, 1);
1609
1610 // Apply std region filter
1611 StdQuadExp::v_SVVLaplacianFilter(array, mkey);
1612
1613 // Divide by sqrt(Jac)
1614 Vmath::Vdiv(nq, array, 1, sqrt_jac, 1, array, 1);
1615}
1616
1617/** @brief: This method gets all of the factors which are
1618 required as part of the Gradient Jump Penalty (GJP)
1619 stabilisation and involves the product of the normal and
1620 geometric factors along the element trace.
1621*/
1623 Array<OneD, Array<OneD, NekDouble>> &d0factors,
1624 Array<OneD, Array<OneD, NekDouble>> &d1factors,
1625 [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &d2factors)
1626{
1627 int nquad0 = GetNumPoints(0);
1628 int nquad1 = GetNumPoints(1);
1629
1631 m_metricinfo->GetDerivFactors(GetPointsKeys());
1632
1633 if (d0factors.size() != 4)
1634 {
1635 d0factors = Array<OneD, Array<OneD, NekDouble>>(4);
1636 d1factors = Array<OneD, Array<OneD, NekDouble>>(4);
1637 }
1638
1639 if (d0factors[0].size() != nquad0)
1640 {
1641 d0factors[0] = Array<OneD, NekDouble>(nquad0);
1642 d0factors[2] = Array<OneD, NekDouble>(nquad0);
1643 d1factors[0] = Array<OneD, NekDouble>(nquad0);
1644 d1factors[2] = Array<OneD, NekDouble>(nquad0);
1645 }
1646
1647 if (d0factors[1].size() != nquad1)
1648 {
1649 d0factors[1] = Array<OneD, NekDouble>(nquad1);
1650 d0factors[3] = Array<OneD, NekDouble>(nquad1);
1651 d1factors[1] = Array<OneD, NekDouble>(nquad1);
1652 d1factors[3] = Array<OneD, NekDouble>(nquad1);
1653 }
1654
1655 // Outwards normals
1657 GetTraceNormal(0);
1659 GetTraceNormal(1);
1661 GetTraceNormal(2);
1663 GetTraceNormal(3);
1664
1665 int ncoords = normal_0.size();
1666
1667 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1668 {
1669 // needs checking for 3D coords
1670
1671 // factors 0 and 2
1672 for (int i = 0; i < nquad0; ++i)
1673 {
1674 d0factors[0][i] = df[0][i] * normal_0[0][i];
1675 d0factors[2][i] = df[0][nquad0 * (nquad1 - 1) + i] * normal_2[0][i];
1676
1677 d1factors[0][i] = df[1][i] * normal_0[0][i];
1678 d1factors[2][i] = df[1][nquad0 * (nquad1 - 1) + i] * normal_2[0][i];
1679 }
1680
1681 for (int n = 1; n < ncoords; ++n)
1682 {
1683 // d xi_1/dy n_y
1684 // needs checking for 3D coords
1685 for (int i = 0; i < nquad0; ++i)
1686 {
1687 d0factors[0][i] += df[2 * n][i] * normal_0[n][i];
1688 d0factors[2][i] +=
1689 df[2 * n][nquad0 * (nquad1 - 1) + i] * normal_2[n][i];
1690
1691 d1factors[0][i] += df[2 * n + 1][i] * normal_0[n][i];
1692 d1factors[2][i] +=
1693 df[2 * n + 1][nquad0 * (nquad1 - 1) + i] * normal_2[n][i];
1694 }
1695 }
1696
1697 // faces 1 and 3
1698 for (int i = 0; i < nquad1; ++i)
1699 {
1700 d0factors[1][i] = df[0][(i + 1) * nquad0 - 1] * normal_1[0][i];
1701 d0factors[3][i] = df[0][i * nquad0] * normal_3[0][i];
1702
1703 d1factors[1][i] = df[1][(i + 1) * nquad0 - 1] * normal_1[0][i];
1704 d1factors[3][i] = df[1][i * nquad0] * normal_3[0][i];
1705 }
1706
1707 for (int n = 1; n < ncoords; ++n)
1708 {
1709 for (int i = 0; i < nquad1; ++i)
1710 {
1711 d0factors[1][i] +=
1712 df[2 * n][(i + 1) * nquad0 - 1] * normal_1[n][i];
1713 d0factors[3][i] += df[2 * n][i * nquad0] * normal_3[n][i];
1714
1715 d1factors[1][i] +=
1716 df[2 * n + 1][(i + 1) * nquad0 - 1] * normal_1[n][i];
1717 d1factors[3][i] += df[2 * n + 1][i * nquad0] * normal_3[n][i];
1718 }
1719 }
1720 }
1721 else
1722 {
1723 // d xi_2/dx n_x
1724 for (int i = 0; i < nquad0; ++i)
1725 {
1726 d1factors[0][i] = df[1][0] * normal_0[0][i];
1727 d1factors[2][i] = df[1][0] * normal_2[0][i];
1728 }
1729
1730 // d xi_1/dx n_x
1731 for (int i = 0; i < nquad1; ++i)
1732 {
1733 d0factors[1][i] = df[0][0] * normal_1[0][i];
1734 d0factors[3][i] = df[0][0] * normal_3[0][i];
1735 }
1736
1737 for (int n = 1; n < ncoords; ++n)
1738 {
1739 // d xi_2/dy n_y
1740 // needs checking for 3D coords
1741 for (int i = 0; i < nquad0; ++i)
1742 {
1743 d1factors[0][i] += df[2 * n + 1][0] * normal_0[n][i];
1744 d1factors[2][i] += df[2 * n + 1][0] * normal_2[n][i];
1745 }
1746
1747 // d xi_1/dy n_y
1748 // needs checking for 3D coords
1749 for (int i = 0; i < nquad1; ++i)
1750 {
1751 d0factors[1][i] += df[2 * n][0] * normal_1[n][i];
1752 d0factors[3][i] += df[2 * n][0] * normal_3[n][i];
1753 }
1754 }
1755
1756 // d1factors
1757 // d xi_1/dx n_x
1758 for (int i = 0; i < nquad0; ++i)
1759 {
1760 d0factors[0][i] = df[0][0] * normal_0[0][i];
1761 d0factors[2][i] = df[0][0] * normal_2[0][i];
1762 }
1763
1764 // d xi_2/dx n_x
1765 for (int i = 0; i < nquad1; ++i)
1766 {
1767 d1factors[1][i] = df[1][0] * normal_1[0][i];
1768 d1factors[3][i] = df[1][0] * normal_3[0][i];
1769 }
1770
1771 for (int n = 1; n < ncoords; ++n)
1772 {
1773 // d xi_1/dy n_y
1774 // needs checking for 3D coords
1775 for (int i = 0; i < nquad0; ++i)
1776 {
1777 d0factors[0][i] += df[2 * n][0] * normal_0[n][i];
1778 d0factors[2][i] += df[2 * n][0] * normal_2[n][i];
1779 }
1780
1781 // d xi_2/dy n_y
1782 // needs checking for 3D coords
1783 for (int i = 0; i < nquad1; ++i)
1784 {
1785 d1factors[1][i] += df[2 * n + 1][0] * normal_1[n][i];
1786 d1factors[3][i] += df[2 * n + 1][0] * normal_3[n][i];
1787 }
1788 }
1789 }
1790}
1791} // 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: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
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
Definition: QuadExp.cpp:587
void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1404
StdRegions::Orientation v_GetTraceOrient(int edge) override
Definition: QuadExp.cpp:1349
DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey) override
Definition: QuadExp.cpp:1394
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:1384
StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const override
Definition: QuadExp.cpp:529
void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1448
void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords) override
Definition: QuadExp.cpp:547
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:1553
void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp) override
Definition: QuadExp.cpp:1495
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:660
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:488
void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1411
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:576
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:598
void v_WeakDirectionalDerivMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1434
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:1389
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:1259
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
Definition: QuadExp.cpp:540
void v_ComputeTraceNormal(const int edge) override
Definition: QuadExp.cpp:980
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:523
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1354
void v_DropLocStaticCondMatrix(const MatrixKey &mkey) override
Definition: QuadExp.cpp:1399
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:1441
DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1374
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1590
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:1622
NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals) override
Definition: QuadExp.cpp:568
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:1455
void v_GetTracePhysMap(const int edge, Array< OneD, int > &outarray) override
Definition: QuadExp.cpp:724
void v_GetTraceQFactors(const int edge, Array< OneD, NekDouble > &outarray) override
Definition: QuadExp.cpp:766
void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: QuadExp.cpp:1426
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: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
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space.
Definition: StdExpansion.h:433
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: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
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: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< StdQuadExp > StdQuadExpSharedPtr
Definition: StdQuadExp.h:221
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:430
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
STL namespace.
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285