Nektar++
Loading...
Searching...
No Matches
Expansion2D.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: Expansion2D.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: File for Expansion2D routines
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <LibUtilities/Foundations/Interp.h>
36#include <LibUtilities/Foundations/InterpCoeff.h>
37#include <LibUtilities/Foundations/ManagerAccess.h>
38#include <LocalRegions/Expansion1D.h>
39#include <LocalRegions/Expansion2D.h>
40#include <LocalRegions/Expansion3D.h>
41#include <LocalRegions/MatrixKey.h>
42#include <LocalRegions/SegExp.h>
43#include <SpatialDomains/Geometry.h>
44#include <SpatialDomains/Geometry2D.h>
45
46using namespace std;
47
48namespace Nektar::LocalRegions
49{
50Expansion2D::Expansion2D(SpatialDomains::Geometry2D *pGeom)
51 : StdExpansion(), Expansion(pGeom), StdExpansion2D()
52{
53}
54
55DNekScalMatSharedPtr Expansion2D::CreateMatrix(const MatrixKey &mkey)
56{
57 DNekScalMatSharedPtr returnval;
58 LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
59
60 ASSERTL2(m_geomFactors->GetGtype() != SpatialDomains::eNoGeomType,
61 "Geometric information is not set up");
62
63 switch (mkey.GetMatrixType())
64 {
65 case StdRegions::eMass:
66 {
67 if ((m_geomFactors->GetGtype() == SpatialDomains::eDeformed) ||
68 mkey.HasVarCoeffForMatrixType(StdRegions::eMass))
69 {
70 NekDouble one = 1.0;
71 DNekMatSharedPtr mat = GenMatrix(mkey);
72
73 returnval =
74 MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
75 }
76 else
77 {
78 NekDouble jac = (m_geomFactors->GetJac())[0];
79 DNekMatSharedPtr mat = GetStdMatrix(mkey);
80
81 returnval =
82 MemoryManager<DNekScalMat>::AllocateSharedPtr(jac, mat);
83 }
84 }
85 break;
86 case StdRegions::eMassGJP:
87 {
88 MatrixKey masskey(mkey, StdRegions::eMass);
89 DNekScalMat &MassMat = *GetLocMatrix(masskey);
90
91 // Generate a local copy of traceMat
92 MatrixKey key(mkey, StdRegions::eNormDerivOnTrace);
93 DNekMatSharedPtr NDTraceMat = Expansion2D::v_GenMatrix(key);
94
95 ASSERTL1(mkey.ConstFactorExists(StdRegions::eFactorGJP),
96 "Need to specify eFactorGJP to construct "
97 "a MassGJP matrix");
98
99 NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorGJP);
100
101 factor /= MassMat.Scale();
102
103 int ntot = MassMat.GetRows() * MassMat.GetColumns();
104
105 Vmath::Svtvp(ntot, factor, &NDTraceMat->GetPtr()[0], 1,
106 MassMat.GetRawPtr(), 1, &NDTraceMat->GetPtr()[0], 1);
107
108 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(
109 MassMat.Scale(), NDTraceMat);
110 }
111 break;
112 case StdRegions::eInvMass:
113 {
114 if (m_geomFactors->GetGtype() == SpatialDomains::eDeformed)
115 {
116 NekDouble one = 1.0;
117 StdRegions::StdMatrixKey masskey(StdRegions::eMass,
118 DetShapeType(), *this);
119 DNekMatSharedPtr mat = GenMatrix(masskey);
120 mat->Invert();
121
122 returnval =
123 MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
124 }
125 else
126 {
127 NekDouble fac = 1.0 / (m_geomFactors->GetJac())[0];
128 DNekMatSharedPtr mat = GetStdMatrix(mkey);
129
130 returnval =
131 MemoryManager<DNekScalMat>::AllocateSharedPtr(fac, mat);
132 }
133 }
134 break;
135 case StdRegions::eWeakDeriv0:
136 case StdRegions::eWeakDeriv1:
137 case StdRegions::eWeakDeriv2:
138 {
139 if (m_geomFactors->GetGtype() == SpatialDomains::eDeformed ||
140 mkey.HasVarCoeffForMatrixType(StdRegions::eLaplacian))
141 {
142 NekDouble one = 1.0;
143 DNekMatSharedPtr mat = GenMatrix(mkey);
144
145 returnval =
146 MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
147 }
148 else
149 {
150 NekDouble jac = (m_geomFactors->GetJac())[0];
151 Array<TwoD, const NekDouble> df =
152 m_geomFactors->GetDerivFactors();
153 int dir = 0;
154 if (mkey.GetMatrixType() == StdRegions::eWeakDeriv0)
155 {
156 dir = 0;
157 }
158 if (mkey.GetMatrixType() == StdRegions::eWeakDeriv1)
159 {
160 dir = 1;
161 }
162 if (mkey.GetMatrixType() == StdRegions::eWeakDeriv2)
163 {
164 dir = 2;
165 }
166
167 MatrixKey deriv0key(StdRegions::eWeakDeriv0,
168 mkey.GetShapeType(), *this);
169 MatrixKey deriv1key(StdRegions::eWeakDeriv1,
170 mkey.GetShapeType(), *this);
171
172 DNekMat &deriv0 = *GetStdMatrix(deriv0key);
173 DNekMat &deriv1 = *GetStdMatrix(deriv1key);
174
175 int rows = deriv0.GetRows();
176 int cols = deriv1.GetColumns();
177
178 DNekMatSharedPtr WeakDeriv =
179 MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
180 (*WeakDeriv) =
181 df[2 * dir][0] * deriv0 + df[2 * dir + 1][0] * deriv1;
182
183 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(
184 jac, WeakDeriv);
185 }
186 }
187 break;
188 case StdRegions::eWeakDirectionalDeriv:
189 {
190 if (m_geomFactors->GetGtype() == SpatialDomains::eDeformed ||
191 mkey.GetNVarCoeff())
192 {
193 NekDouble one = 1.0;
194 DNekMatSharedPtr mat = GenMatrix(mkey);
195
196 returnval =
197 MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
198 }
199 else
200 {
201 int shapedim = 2;
202
203 // dfdirxi = tan_{xi_x} * d \xi/dx
204 // + tan_{xi_y} * d \xi/dy
205 // + tan_{xi_z} * d \xi/dz
206 // dfdireta = tan_{eta_x} * d \eta/dx
207 // + tan_{xi_y} * d \xi/dy
208 // + tan_{xi_z} * d \xi/dz
209 NekDouble jac = (m_geomFactors->GetJac())[0];
210 Array<TwoD, const NekDouble> df =
211 m_geomFactors->GetDerivFactors();
212
213 Array<OneD, NekDouble> direction =
214 mkey.GetVarCoeff(StdRegions::eVarCoeffMF);
215
216 // d / dx = df[0]*deriv0 + df[1]*deriv1
217 // d / dy = df[2]*deriv0 + df[3]*deriv1
218 // d / dz = df[4]*deriv0 + df[5]*deriv1
219
220 // dfdir[dir] = e \cdot (d/dx, d/dy, d/dz)
221 // = (e^0 * df[0] + e^1 * df[2]
222 // + e^2 * df[4]) * deriv0
223 // + (e^0 * df[1] + e^1 * df[3]
224 // + e^2 * df[5]) * deriv1
225 // dfdir[dir] = e^0 * df[2 * 0 + dir]
226 // + e^1 * df[2 * 1 + dir]
227 // + e^2 * df [ 2 * 2 + dir]
228 Array<OneD, Array<OneD, NekDouble>> dfdir(shapedim);
229 Expansion::ComputeGmatcdotMF(df, direction, dfdir);
230
231 StdRegions::VarCoeffMap dfdirxi;
232 StdRegions::VarCoeffMap dfdireta;
233
234 dfdirxi[StdRegions::eVarCoeffWeakDeriv] = dfdir[0];
235 dfdireta[StdRegions::eVarCoeffWeakDeriv] = dfdir[1];
236
237 MatrixKey derivxikey(StdRegions::eWeakDeriv0,
238 mkey.GetShapeType(), *this,
239 StdRegions::NullConstFactorMap, dfdirxi);
240 MatrixKey derivetakey(StdRegions::eWeakDeriv1,
241 mkey.GetShapeType(), *this,
242 StdRegions::NullConstFactorMap, dfdireta);
243
244 DNekMat &derivxi = *GetStdMatrix(derivxikey);
245 DNekMat &deriveta = *GetStdMatrix(derivetakey);
246
247 int rows = derivxi.GetRows();
248 int cols = deriveta.GetColumns();
249
250 DNekMatSharedPtr WeakDirDeriv =
251 MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
252
253 (*WeakDirDeriv) = derivxi + deriveta;
254
255 // Add (\nabla \cdot e^k ) Mass
256 StdRegions::VarCoeffMap DiveMass;
257 DiveMass[StdRegions::eVarCoeffMass] =
258 mkey.GetVarCoeff(StdRegions::eVarCoeffMFDiv);
259 StdRegions::StdMatrixKey stdmasskey(
260 StdRegions::eMass, mkey.GetShapeType(), *this,
261 StdRegions::NullConstFactorMap, DiveMass);
262
263 DNekMatSharedPtr DiveMassmat = GetStdMatrix(stdmasskey);
264
265 (*WeakDirDeriv) = (*WeakDirDeriv) + (*DiveMassmat);
266
267 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(
268 jac, WeakDirDeriv);
269 }
270 break;
271 }
272 case StdRegions::eLaplacian:
273 {
274 if (m_geomFactors->GetGtype() == SpatialDomains::eDeformed ||
275 mkey.HasVarCoeffForMatrixType(StdRegions::eLaplacian) ||
276 mkey.ConstFactorExists(StdRegions::eFactorSVVDiffCoeff))
277 {
278 NekDouble one = 1.0;
279 DNekMatSharedPtr mat = GenMatrix(mkey);
280
281 returnval =
282 MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
283 }
284 else
285 {
286 MatrixKey lap00key(StdRegions::eLaplacian00,
287 mkey.GetShapeType(), *this);
288 MatrixKey lap01key(StdRegions::eLaplacian01,
289 mkey.GetShapeType(), *this);
290 MatrixKey lap11key(StdRegions::eLaplacian11,
291 mkey.GetShapeType(), *this);
292
293 DNekMat &lap00 = *GetStdMatrix(lap00key);
294 DNekMat &lap01 = *GetStdMatrix(lap01key);
295 DNekMat &lap11 = *GetStdMatrix(lap11key);
296
297 NekDouble jac = (m_geomFactors->GetJac())[0];
298 Array<TwoD, const NekDouble> gmat =
299 m_geomFactors->GetGmat(ptsKeys);
300
301 int rows = lap00.GetRows();
302 int cols = lap00.GetColumns();
303
304 DNekMatSharedPtr lap =
305 MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
306
307 (*lap) = gmat[0][0] * lap00 +
308 gmat[1][0] * (lap01 + Transpose(lap01)) +
309 gmat[3][0] * lap11;
310
311 returnval =
312 MemoryManager<DNekScalMat>::AllocateSharedPtr(jac, lap);
313 }
314 }
315 break;
316 case StdRegions::eInvLaplacianWithUnityMean:
317 {
318 DNekMatSharedPtr mat = GenMatrix(mkey);
319
320 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, mat);
321 }
322 break;
323 case StdRegions::eHelmholtz:
324 {
325 NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorLambda);
326
327 // Construct mass matrix
328 // Check for mass-specific varcoeffs to avoid unncessary
329 // re-computation of the elemental matrix every time step
330 StdRegions::VarCoeffMap massVarcoeffs = StdRegions::NullVarCoeffMap;
331 if (mkey.HasVarCoeffForMatrixType(StdRegions::eMass))
332 {
333 massVarcoeffs[StdRegions::eVarCoeffMass] =
334 mkey.GetVarCoeff(StdRegions::eVarCoeffMass);
335 }
336 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this,
337 mkey.GetConstFactors(), massVarcoeffs);
338 DNekScalMat &MassMat = *GetLocMatrix(masskey);
339
340 // Construct laplacian matrix
341 // Take all varcoeffs if one or more are detected
342 // use mapping from MatrixType to Vector of Varcoeffs
343 StdRegions::VarCoeffMap lapVarcoeffs = StdRegions::NullVarCoeffMap;
344 if (mkey.HasVarCoeffForMatrixType(StdRegions::eLaplacian))
345 {
346 lapVarcoeffs = mkey.GetVarCoeffs();
347 }
348 MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
349 mkey.GetConstFactors(), lapVarcoeffs);
350 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
351
352 int rows = LapMat.GetRows();
353 int cols = LapMat.GetColumns();
354
355 DNekMatSharedPtr helm =
356 MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
357
358 NekDouble one = 1.0;
359 (*helm) = LapMat + factor * MassMat;
360
361 returnval =
362 MemoryManager<DNekScalMat>::AllocateSharedPtr(one, helm);
363
364 // Only drop matrix if time-dependence possible
365 if (!massVarcoeffs.empty())
366 {
367 DropLocMatrix(masskey);
368 }
369 if (!lapVarcoeffs.empty())
370 {
371 DropLocMatrix(lapkey);
372 }
373 }
374 break;
375 case StdRegions::eHelmholtzGJP:
376 {
377 MatrixKey helmkey(mkey, StdRegions::eHelmholtz);
378 DNekScalMat &HelmMat = *GetLocMatrix(helmkey);
379
380 // Generate a local copy of traceMat
381 MatrixKey key(mkey, StdRegions::eNormDerivOnTrace);
382 DNekMatSharedPtr NDTraceMat = Expansion2D::v_GenMatrix(key);
383
384 ASSERTL1(mkey.ConstFactorExists(StdRegions::eFactorGJP),
385 "Need to specify eFactorGJP to construct "
386 "a HelmholtzGJP matrix");
387
388 NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorGJP);
389
390 factor /= HelmMat.Scale();
391
392 int ntot = HelmMat.GetRows() * HelmMat.GetColumns();
393
394 Vmath::Svtvp(ntot, factor, &NDTraceMat->GetPtr()[0], 1,
395 HelmMat.GetRawPtr(), 1, &NDTraceMat->GetPtr()[0], 1);
396
397 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(
398 HelmMat.Scale(), NDTraceMat);
399 }
400 break;
401 case StdRegions::eLinearAdvectionReaction:
402 {
403 NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
404
405 // Construct mass matrix
406 // Check for mass-specific varcoeffs to avoid unncessary
407 // re-computation of the elemental matrix every time step
408 StdRegions::VarCoeffMap massVarcoeffs = StdRegions::NullVarCoeffMap;
409 if (mkey.HasVarCoeffForMatrixType(StdRegions::eMass))
410 {
411 massVarcoeffs[StdRegions::eVarCoeffMass] =
412 mkey.GetVarCoeff(StdRegions::eVarCoeffMass);
413 }
414 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this,
415 mkey.GetConstFactors(), massVarcoeffs);
416 DNekScalMat &MassMat = *GetLocMatrix(masskey);
417
418 // Construct advection matrix
419 // Check for varcoeffs not required;
420 // assume advection velocity is always time-dependent
421 MatrixKey advkey(mkey, StdRegions::eLinearAdvection);
422 DNekScalMat &AdvMat = *GetLocMatrix(advkey);
423
424 int rows = MassMat.GetRows();
425 int cols = MassMat.GetColumns();
426
427 DNekMatSharedPtr adr =
428 MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
429
430 NekDouble one = 1.0;
431 (*adr) = -lambda * MassMat + AdvMat;
432
433 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, adr);
434
435 // Clear memory for time-dependent matrices
436 DropLocMatrix(advkey);
437 if (!massVarcoeffs.empty())
438 {
439 DropLocMatrix(masskey);
440 }
441 }
442 break;
443 case StdRegions::eLinearAdvectionDiffusionReaction:
444 {
445 NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
446
447 // Construct mass matrix
448 // Check for mass-specific varcoeffs to avoid unncessary
449 // re-computation of the elemental matrix every time step
450 StdRegions::VarCoeffMap massVarcoeffs = StdRegions::NullVarCoeffMap;
451 if (mkey.HasVarCoeff(StdRegions::eVarCoeffMass))
452 {
453 massVarcoeffs[StdRegions::eVarCoeffMass] =
454 mkey.GetVarCoeff(StdRegions::eVarCoeffMass);
455 }
456 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this,
457 mkey.GetConstFactors(), massVarcoeffs);
458 DNekScalMat &MassMat = *GetLocMatrix(masskey);
459
460 // Construct laplacian matrix (Check for varcoeffs)
461 // Take all varcoeffs if one or more are detected
462 // TODO We might want to have a map
463 // from MatrixType to Vector of Varcoeffs and vice-versa
464 StdRegions::VarCoeffMap lapVarcoeffs = StdRegions::NullVarCoeffMap;
465 if (mkey.HasVarCoeffForMatrixType(StdRegions::eLaplacian))
466 {
467 lapVarcoeffs = mkey.GetVarCoeffs();
468 }
469 MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
470 mkey.GetConstFactors(), lapVarcoeffs);
471 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
472
473 // Construct advection matrix
474 // Check for varcoeffs not required;
475 // assume advection velocity is always time-dependent
476 MatrixKey advkey(mkey, StdRegions::eLinearAdvection);
477 DNekScalMat &AdvMat = *GetLocMatrix(advkey);
478
479 int rows = LapMat.GetRows();
480 int cols = LapMat.GetColumns();
481
482 DNekMatSharedPtr adr =
483 MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
484
485 NekDouble one = 1.0;
486 (*adr) = LapMat - lambda * MassMat + AdvMat;
487
488 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, adr);
489
490 // Clear memory for time-dependent matrices
491 DropLocMatrix(advkey);
492 if (!massVarcoeffs.empty())
493 {
494 DropLocMatrix(masskey);
495 }
496 if (!lapVarcoeffs.empty())
497 {
498 DropLocMatrix(lapkey);
499 }
500 }
501 break;
502 case StdRegions::eLinearAdvectionDiffusionReactionGJP:
503 {
504 // Copied mostly from ADR solve to have fine-grain control
505 // over updating only advection matrix, relevant for performance!
506 NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
507
508 // Construct mass matrix (Check for varcoeffs)
509 StdRegions::VarCoeffMap massVarcoeffs = StdRegions::NullVarCoeffMap;
510 if (mkey.HasVarCoeffForMatrixType(StdRegions::eMass))
511 {
512 massVarcoeffs[StdRegions::eVarCoeffMass] =
513 mkey.GetVarCoeff(StdRegions::eVarCoeffMass);
514 }
515 MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this,
516 mkey.GetConstFactors(), massVarcoeffs);
517 DNekScalMat &MassMat = *GetLocMatrix(masskey);
518
519 // Construct laplacian matrix (Check for varcoeffs)
520 StdRegions::VarCoeffMap lapVarcoeffs = StdRegions::NullVarCoeffMap;
521 if (mkey.HasVarCoeffForMatrixType(StdRegions::eLaplacian))
522 {
523 lapVarcoeffs = mkey.GetVarCoeffs();
524 }
525 MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
526 mkey.GetConstFactors(), lapVarcoeffs);
527 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
528
529 // Construct advection matrix
530 // (assume advection velocity defined and non-zero)
531 // Could check L2(AdvectionVelocity) or HasVarCoeff
532 MatrixKey advkey(mkey, StdRegions::eLinearAdvection);
533 DNekScalMat &AdvMat = *GetLocMatrix(advkey);
534
535 // Generate a local copy of traceMat
536 MatrixKey gjpkey(StdRegions::eNormDerivOnTrace, mkey.GetShapeType(),
537 *this, mkey.GetConstFactors());
538 DNekScalMat &NDTraceMat = *GetLocMatrix(gjpkey);
539
540 NekDouble gjpfactor = mkey.GetConstFactor(StdRegions::eFactorGJP);
541 ASSERTL1(mkey.ConstFactorExists(StdRegions::eFactorGJP),
542 "Need to specify eFactorGJP to construct "
543 "a LinearAdvectionDiffusionReactionGJP matrix");
544
545 int rows = LapMat.GetRows();
546 int cols = LapMat.GetColumns();
547
548 DNekMatSharedPtr adr =
549 MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
550
551 NekDouble one = 1.0;
552 (*adr) =
553 LapMat - lambda * MassMat + AdvMat + gjpfactor * NDTraceMat;
554
555 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, adr);
556
557 // Clear memory
558 DropLocMatrix(advkey);
559 DropLocMatrix(masskey);
560 DropLocMatrix(lapkey);
561 }
562 break;
563 case StdRegions::eNormDerivOnTrace:
564 {
565 NekDouble one = 1.0;
566 DNekMatSharedPtr mat = Expansion2D::v_GenMatrix(mkey);
567
568 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
569 }
570 break;
571 case StdRegions::eIProductWRTBase:
572 {
573 if (m_geomFactors->GetGtype() == SpatialDomains::eDeformed)
574 {
575 NekDouble one = 1.0;
576 DNekMatSharedPtr mat = GenMatrix(mkey);
577
578 returnval =
579 MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
580 }
581 else
582 {
583 NekDouble jac = (m_geomFactors->GetJac())[0];
584 DNekMatSharedPtr mat = GetStdMatrix(mkey);
585
586 returnval =
587 MemoryManager<DNekScalMat>::AllocateSharedPtr(jac, mat);
588 }
589 }
590 break;
591 case StdRegions::eIProductWRTDerivBase0:
592 case StdRegions::eIProductWRTDerivBase1:
593 case StdRegions::eIProductWRTDerivBase2:
594 {
595 if (m_geomFactors->GetGtype() == SpatialDomains::eDeformed)
596 {
597 NekDouble one = 1.0;
598 DNekMatSharedPtr mat = GenMatrix(mkey);
599
600 returnval =
601 MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
602 }
603 else
604 {
605 NekDouble jac = (m_geomFactors->GetJac())[0];
606
607 const Array<TwoD, const NekDouble> &df =
608 m_geomFactors->GetDerivFactors();
609 int dir = 0;
610 if (mkey.GetMatrixType() == StdRegions::eIProductWRTDerivBase0)
611 {
612 dir = 0;
613 }
614 if (mkey.GetMatrixType() == StdRegions::eIProductWRTDerivBase1)
615 {
616 dir = 1;
617 }
618 if (mkey.GetMatrixType() == StdRegions::eIProductWRTDerivBase2)
619 {
620 dir = 2;
621 }
622
623 MatrixKey iProdDeriv0Key(StdRegions::eIProductWRTDerivBase0,
624 mkey.GetShapeType(), *this);
625 MatrixKey iProdDeriv1Key(StdRegions::eIProductWRTDerivBase1,
626 mkey.GetShapeType(), *this);
627
628 DNekMat &stdiprod0 = *GetStdMatrix(iProdDeriv0Key);
629 DNekMat &stdiprod1 = *GetStdMatrix(iProdDeriv0Key);
630
631 int rows = stdiprod0.GetRows();
632 int cols = stdiprod1.GetColumns();
633
634 DNekMatSharedPtr mat =
635 MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols);
636 (*mat) =
637 df[2 * dir][0] * stdiprod0 + df[2 * dir + 1][0] * stdiprod1;
638
639 returnval =
640 MemoryManager<DNekScalMat>::AllocateSharedPtr(jac, mat);
641 }
642 }
643 break;
644
645 case StdRegions::eInvHybridDGHelmholtz:
646 {
647 NekDouble one = 1.0;
648
649 MatrixKey hkey(StdRegions::eHybridDGHelmholtz, DetShapeType(),
650 *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
651
652 DNekMatSharedPtr mat = GenMatrix(hkey);
653
654 mat->Invert();
655
656 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
657 }
658 break;
659 case StdRegions::eInterpGauss:
660 {
661 ASSERTL1(DetShapeType() == LibUtilities::eQuadrilateral,
662 "Matrix only setup for quad elements currently");
663 DNekMatSharedPtr m_Ix;
664 Array<OneD, NekDouble> coords(1, 0.0);
665 StdRegions::ConstFactorMap factors = mkey.GetConstFactors();
666 int edge = static_cast<int>(factors[StdRegions::eFactorGaussEdge]);
667
668 coords[0] = (edge == 0 || edge == 3) ? -1.0 : 1.0;
669
670 m_Ix = m_base[(edge + 1) % 2]->GetI(coords);
671
672 returnval =
673 MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, m_Ix);
674 }
675 break;
676 case StdRegions::ePreconLinearSpace:
677 {
678 NekDouble one = 1.0;
679 MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(),
680 *this, mkey.GetConstFactors(),
681 mkey.GetVarCoeffs());
682 DNekScalBlkMatSharedPtr helmStatCond =
683 GetLocStaticCondMatrix(helmkey);
684 DNekScalMatSharedPtr A = helmStatCond->GetBlock(0, 0);
685 DNekMatSharedPtr R = BuildVertexMatrix(A);
686
687 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, R);
688 }
689 break;
690 default:
691 {
692 NekDouble one = 1.0;
693 DNekMatSharedPtr mat = GenMatrix(mkey);
694
695 returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
696 }
697 break;
698 }
699
700 return returnval;
701}
702
703void Expansion2D::v_PhysDeriv(const int dir,
704 const Array<OneD, const NekDouble> &inarray,
705 Array<OneD, NekDouble> &outarray)
706{
707 switch (dir)
708 {
709 case 0:
710 {
711 v_PhysDeriv(inarray, outarray, NullNekDouble1DArray,
712 NullNekDouble1DArray);
713 break;
714 }
715
716 case 1:
717 {
718 v_PhysDeriv(inarray, NullNekDouble1DArray, outarray,
719 NullNekDouble1DArray);
720 break;
721 }
722
723 case 2:
724 {
725 v_PhysDeriv(inarray, NullNekDouble1DArray, NullNekDouble1DArray,
726 outarray);
727 break;
728 }
729
730 default:
731 {
732 ASSERTL1(false, "input dir is out of range");
733 }
734 break;
735 }
736}
737
738void Expansion2D::v_PhysDeriv(const Array<OneD, const NekDouble> &inarray,
739 Array<OneD, NekDouble> &out_d0,
740 Array<OneD, NekDouble> &out_d1,
741 Array<OneD, NekDouble> &out_d2)
742{
743 int nquad0 = m_base[0]->GetNumPoints();
744 int nquad1 = m_base[1]->GetNumPoints();
745 int nqtot = nquad0 * nquad1;
746 const Array<TwoD, const NekDouble> &df = m_geomFactors->GetDerivFactors();
747 Array<OneD, NekDouble> diff0(2 * nqtot);
748 Array<OneD, NekDouble> diff1(diff0 + nqtot);
749
750 v_StdPhysDeriv(inarray, diff0, diff1, NullNekDouble1DArray);
751
752 if (m_geomFactors->GetGtype() == SpatialDomains::eDeformed)
753 {
754 if (out_d0.size())
755 {
756 Vmath::Vmul(nqtot, df[0], 1, diff0, 1, out_d0, 1);
757 Vmath::Vvtvp(nqtot, df[1], 1, diff1, 1, out_d0, 1, out_d0, 1);
758 }
759
760 if (out_d1.size())
761 {
762 Vmath::Vmul(nqtot, df[2], 1, diff0, 1, out_d1, 1);
763 Vmath::Vvtvp(nqtot, df[3], 1, diff1, 1, out_d1, 1, out_d1, 1);
764 }
765
766 if (out_d2.size())
767 {
768 Vmath::Vmul(nqtot, df[4], 1, diff0, 1, out_d2, 1);
769 Vmath::Vvtvp(nqtot, df[5], 1, diff1, 1, out_d2, 1, out_d2, 1);
770 }
771 }
772 else // regular geometry
773 {
774 if (out_d0.size())
775 {
776 Vmath::Smul(nqtot, df[0][0], diff0, 1, out_d0, 1);
777 Blas::Daxpy(nqtot, df[1][0], diff1, 1, out_d0, 1);
778 }
779
780 if (out_d1.size())
781 {
782 Vmath::Smul(nqtot, df[2][0], diff0, 1, out_d1, 1);
783 Blas::Daxpy(nqtot, df[3][0], diff1, 1, out_d1, 1);
784 }
785
786 if (out_d2.size())
787 {
788 Vmath::Smul(nqtot, df[4][0], diff0, 1, out_d2, 1);
789 Blas::Daxpy(nqtot, df[5][0], diff1, 1, out_d2, 1);
790 }
791 }
792}
793
794void Expansion2D::v_PhysDirectionalDeriv(
795 const Array<OneD, const NekDouble> &inarray,
796 const Array<OneD, const NekDouble> &direction,
797 Array<OneD, NekDouble> &outarray)
798{
799 int nquad0 = m_base[0]->GetNumPoints();
800 int nquad1 = m_base[1]->GetNumPoints();
801 int nqtot = nquad0 * nquad1;
802
803 const Array<TwoD, const NekDouble> &df = m_geomFactors->GetDerivFactors();
804
805 Array<OneD, NekDouble> diff0(2 * nqtot);
806 Array<OneD, NekDouble> diff1(diff0 + nqtot);
807
808 // diff0 = du/d_xi, diff1 = du/d_eta
809 v_StdPhysDeriv(inarray, diff0, diff1, NullNekDouble1DArray);
810
811 if (m_geomFactors->GetGtype() == SpatialDomains::eDeformed)
812 {
813 Array<OneD, Array<OneD, NekDouble>> tangmat(2);
814
815 // D^v_xi = v_x*d_xi/dx + v_y*d_xi/dy + v_z*d_xi/dz
816 // D^v_eta = v_x*d_eta/dx + v_y*d_eta/dy + v_z*d_eta/dz
817 for (int i = 0; i < 2; ++i)
818 {
819 tangmat[i] = Array<OneD, NekDouble>(nqtot, 0.0);
820 for (int k = 0; k < (m_geom->GetCoordim()); ++k)
821 {
822 Vmath::Vvtvp(nqtot, &df[2 * k + i][0], 1, &direction[k * nqtot],
823 1, &tangmat[i][0], 1, &tangmat[i][0], 1);
824 }
825 }
826
827 /// D_v = D^v_xi * du/d_xi + D^v_eta * du/d_eta
828 Vmath::Vmul(nqtot, &tangmat[0][0], 1, &diff0[0], 1, &outarray[0], 1);
829 Vmath::Vvtvp(nqtot, &tangmat[1][0], 1, &diff1[0], 1, &outarray[0], 1,
830 &outarray[0], 1);
831 }
832 else
833 {
834 Array<OneD, Array<OneD, NekDouble>> tangmat(2);
835
836 for (int i = 0; i < 2; ++i)
837 {
838 tangmat[i] = Array<OneD, NekDouble>(nqtot, 0.0);
839 for (int k = 0; k < (m_geom->GetCoordim()); ++k)
840 {
841 Vmath::Svtvp(nqtot, df[2 * k + i][0], &direction[k * nqtot], 1,
842 &tangmat[i][0], 1, &tangmat[i][0], 1);
843 }
844 }
845
846 /// D_v = D^v_xi * du/d_xi + D^v_eta * du/d_eta
847 Vmath::Vmul(nqtot, &tangmat[0][0], 1, &diff0[0], 1, &outarray[0], 1);
848
849 Vmath::Vvtvp(nqtot, &tangmat[1][0], 1, &diff1[0], 1, &outarray[0], 1,
850 &outarray[0], 1);
851 }
852}
853
854void Expansion2D::v_IProductWRTBase(const Array<OneD, const NekDouble> &inarray,
855 Array<OneD, NekDouble> &outarray)
856{
857 const bool CollDir0 = m_base[0]->Collocation();
858 const bool CollDir1 = m_base[1]->Collocation();
859
860 const Array<OneD, const NekDouble> &jac = m_geomFactors->GetJac();
861 bool Deformed = (m_geomFactors->GetGtype() == SpatialDomains::eDeformed);
862
863 if (v_IsCollocatedBasis())
864 {
865 int nqtot = GetTotPoints();
866 if (Deformed)
867 {
868 Vmath::Vmul(nqtot, jac, 1, inarray, 1, outarray, 1);
869 }
870 else
871 {
872 Vmath::Smul(nqtot, jac[0], inarray, 1, outarray, 1);
873 }
874 v_MultiplyByStdQuadratureMetric(outarray, outarray);
875 }
876 else
877 {
878 v_IProductWRTBaseKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(),
879 inarray, outarray, jac, Deformed, CollDir0,
880 CollDir1);
881 }
882}
883
884void Expansion2D::v_AddEdgeNormBoundaryInt(
885 const int edge, const ExpansionSharedPtr &EdgeExp,
886 const Array<OneD, const NekDouble> &Fx,
887 const Array<OneD, const NekDouble> &Fy, Array<OneD, NekDouble> &outarray)
888{
889 ASSERTL1(GetCoordim() == 2, "Routine only set up for two-dimensions");
890
891 const Array<OneD, const Array<OneD, NekDouble>> normals =
892 GetTraceNormal(edge);
893
894 if (m_requireNeg.size() == 0)
895 {
896 int nedges = GetNtraces();
897 m_requireNeg.resize(nedges);
898
899 for (int i = 0; i < nedges; ++i)
900 {
901 m_requireNeg[i] = false;
902
903 ExpansionSharedPtr edgeExp = m_traceExp[i].lock();
904
905 if (edgeExp->GetRightAdjacentElementExp())
906 {
907 if (edgeExp->GetRightAdjacentElementExp()
908 ->GetGeom()
909 ->GetGlobalID() == GetGeom()->GetGlobalID())
910 {
911 m_requireNeg[i] = true;
912 }
913 }
914 }
915 }
916
917 // We allow the case of mixed polynomial order by supporting only
918 // those modes on the edge common to both adjoining elements. This
919 // is enforced here by taking the minimum size and padding with
920 // zeros.
921 int nquad_e = min(EdgeExp->GetNumPoints(0), int(normals[0].size()));
922
923 int nEdgePts = EdgeExp->GetTotPoints();
924 Array<OneD, NekDouble> edgePhys(nEdgePts);
925 Vmath::Vmul(nquad_e, normals[0], 1, Fx, 1, edgePhys, 1);
926 Vmath::Vvtvp(nquad_e, normals[1], 1, Fy, 1, edgePhys, 1, edgePhys, 1);
927
928 Expansion1DSharedPtr locExp = EdgeExp->as<Expansion1D>();
929
930 if (m_requireNeg[edge])
931 {
932 if (locExp->GetRightAdjacentElementExp()->GetGeom()->GetGlobalID() ==
933 m_geom->GetGlobalID())
934 {
935 Vmath::Neg(nquad_e, edgePhys, 1);
936 }
937 }
938
939 AddEdgeNormBoundaryInt(edge, EdgeExp, edgePhys, outarray);
940}
941
942void Expansion2D::v_AddEdgeNormBoundaryInt(
943 const int edge, const ExpansionSharedPtr &EdgeExp,
944 const Array<OneD, const NekDouble> &Fn, Array<OneD, NekDouble> &outarray)
945{
946 int i;
947
948 if (m_requireNeg.size() == 0)
949 {
950 int nedges = GetNtraces();
951 m_requireNeg.resize(nedges);
952
953 for (i = 0; i < nedges; ++i)
954 {
955 m_requireNeg[i] = false;
956
957 ExpansionSharedPtr edgeExp = m_traceExp[i].lock();
958
959 if (edgeExp->GetRightAdjacentElementExp())
960 {
961 if (edgeExp->GetRightAdjacentElementExp()
962 ->GetGeom()
963 ->GetGlobalID() == GetGeom()->GetGlobalID())
964 {
965 m_requireNeg[i] = true;
966 }
967 }
968 }
969 }
970
971 IndexMapKey ikey(eEdgeToElement, DetShapeType(), GetBasisNumModes(0),
972 GetBasisNumModes(1), 0, edge, GetTraceOrient(edge));
973
974 IndexMapValuesSharedPtr map = GetIndexMap(ikey);
975
976 // Order of the element
977 int order_e = map->size();
978 // Order of the trace
979 int n_coeffs = EdgeExp->GetNcoeffs();
980
981 Array<OneD, NekDouble> edgeCoeffs(n_coeffs);
982 if (n_coeffs != order_e) // Going to orthogonal space
983 {
984 EdgeExp->FwdTrans(Fn, edgeCoeffs);
985 Expansion1DSharedPtr locExp = EdgeExp->as<Expansion1D>();
986
987 if (m_requireNeg[edge])
988 {
989 Vmath::Neg(n_coeffs, edgeCoeffs, 1);
990 }
991
992 Array<OneD, NekDouble> coeff(n_coeffs, 0.0);
993 LibUtilities::BasisType btype =
994 ((LibUtilities::BasisType)1); // 1-->Ortho_A
995 LibUtilities::BasisKey bkey_ortho(btype,
996 EdgeExp->GetBasis(0)->GetNumModes(),
997 EdgeExp->GetBasis(0)->GetPointsKey());
998 LibUtilities::BasisKey bkey(EdgeExp->GetBasis(0)->GetBasisType(),
999 EdgeExp->GetBasis(0)->GetNumModes(),
1000 EdgeExp->GetBasis(0)->GetPointsKey());
1001 LibUtilities::InterpCoeff1D(bkey, edgeCoeffs, bkey_ortho, coeff);
1002
1003 // Cutting high frequencies
1004 for (i = order_e; i < n_coeffs; i++)
1005 {
1006 coeff[i] = 0.0;
1007 }
1008
1009 LibUtilities::InterpCoeff1D(bkey_ortho, coeff, bkey, edgeCoeffs);
1010
1011 StdRegions::StdMatrixKey masskey(StdRegions::eMass,
1012 LibUtilities::eSegment, *EdgeExp);
1013 EdgeExp->MassMatrixOp(edgeCoeffs, edgeCoeffs, masskey);
1014 }
1015 else
1016 {
1017 EdgeExp->IProductWRTBase(Fn, edgeCoeffs);
1018
1019 Expansion1DSharedPtr locExp = EdgeExp->as<Expansion1D>();
1020
1021 if (m_requireNeg[edge])
1022 {
1023 Vmath::Neg(n_coeffs, edgeCoeffs, 1);
1024 }
1025 }
1026
1027 // Implementation for all the basis except Gauss points
1028 if (EdgeExp->GetBasis(0)->GetBasisType() != LibUtilities::eGauss_Lagrange)
1029 {
1030 // add data to outarray if forward edge normal is outwards
1031 for (i = 0; i < order_e; ++i)
1032 {
1033 outarray[(*map)[i].index] += (*map)[i].sign * edgeCoeffs[i];
1034 }
1035 }
1036 else
1037 {
1038 int nCoeffs0, nCoeffs1;
1039 int j;
1040
1041 StdRegions::ConstFactorMap factors;
1042 factors[StdRegions::eFactorGaussEdge] = edge;
1043 StdRegions::StdMatrixKey key(StdRegions::eGaussDG, DetShapeType(),
1044 *this, factors);
1045
1046 DNekMatSharedPtr mat_gauss = m_stdMatrixManager[key];
1047
1048 switch (edge)
1049 {
1050 case 0:
1051 {
1052 nCoeffs1 = m_base[1]->GetNumModes();
1053
1054 for (i = 0; i < order_e; ++i)
1055 {
1056 for (j = 0; j < nCoeffs1; j++)
1057 {
1058 outarray[(*map)[i].index + j * order_e] +=
1059 mat_gauss->GetPtr()[j] * (*map)[i].sign *
1060 edgeCoeffs[i];
1061 }
1062 }
1063 break;
1064 }
1065 case 1:
1066 {
1067 nCoeffs0 = m_base[0]->GetNumModes();
1068
1069 for (i = 0; i < order_e; ++i)
1070 {
1071 for (j = 0; j < nCoeffs0; j++)
1072 {
1073 outarray[(*map)[i].index - j] +=
1074 mat_gauss->GetPtr()[order_e - 1 - j] *
1075 (*map)[i].sign * edgeCoeffs[i];
1076 }
1077 }
1078 break;
1079 }
1080 case 2:
1081 {
1082 nCoeffs1 = m_base[1]->GetNumModes();
1083
1084 for (i = 0; i < order_e; ++i)
1085 {
1086 for (j = 0; j < nCoeffs1; j++)
1087 {
1088 outarray[(*map)[i].index - j * order_e] +=
1089 mat_gauss->GetPtr()[order_e - 1 - j] *
1090 (*map)[i].sign * edgeCoeffs[i];
1091 }
1092 }
1093 break;
1094 }
1095 case 3:
1096 {
1097 nCoeffs0 = m_base[0]->GetNumModes();
1098
1099 for (i = 0; i < order_e; ++i)
1100 {
1101 for (j = 0; j < nCoeffs0; j++)
1102 {
1103 outarray[(*map)[i].index + j] +=
1104 mat_gauss->GetPtr()[j] * (*map)[i].sign *
1105 edgeCoeffs[i];
1106 }
1107 }
1108 break;
1109 }
1110 default:
1111 ASSERTL0(false, "edge value (< 3) is out of range");
1112 break;
1113 }
1114 }
1115}
1116
1117void Expansion2D::SetTraceToGeomOrientation(
1118 Array<OneD, ExpansionSharedPtr> &EdgeExp, Array<OneD, NekDouble> &inout)
1119{
1120 int i, cnt = 0;
1121 int nedges = GetNtraces();
1122 Array<OneD, NekDouble> e_tmp;
1123
1124 for (i = 0; i < nedges; ++i)
1125 {
1126 EdgeExp[i]->SetCoeffsToOrientation(
1127 GetTraceOrient(i), e_tmp = inout + cnt, e_tmp = inout + cnt);
1128 cnt += GetTraceNcoeffs(i);
1129 }
1130}
1131
1132/**
1133 * Computes the C matrix entries due to the presence of the identity
1134 * matrix in Eqn. 32.
1135 */
1136void Expansion2D::AddNormTraceInt(const int dir,
1137 Array<OneD, const NekDouble> &inarray,
1138 Array<OneD, ExpansionSharedPtr> &EdgeExp,
1139 Array<OneD, NekDouble> &outarray,
1140 const StdRegions::VarCoeffMap &varcoeffs)
1141{
1142 int i, e, cnt;
1143 int order_e, nquad_e;
1144 int nedges = GetNtraces();
1145
1146 cnt = 0;
1147 for (e = 0; e < nedges; ++e)
1148 {
1149 order_e = EdgeExp[e]->GetNcoeffs();
1150 nquad_e = EdgeExp[e]->GetNumPoints(0);
1151
1152 const Array<OneD, const Array<OneD, NekDouble>> &normals =
1153 GetTraceNormal(e);
1154 Array<OneD, NekDouble> edgeCoeffs(order_e);
1155 Array<OneD, NekDouble> edgePhys(nquad_e);
1156
1157 for (i = 0; i < order_e; ++i)
1158 {
1159 edgeCoeffs[i] = inarray[i + cnt];
1160 }
1161 cnt += order_e;
1162
1163 EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1164
1165 // Multiply by variable coefficient
1166 /// @TODO: Document this
1167 // StdRegions::VarCoeffType VarCoeff[3] = {StdRegions::eVarCoeffD00,
1168 // StdRegions::eVarCoeffD11,
1169 // StdRegions::eVarCoeffD22};
1170 // StdRegions::VarCoeffMap::const_iterator x;
1171 // Array<OneD, NekDouble> varcoeff_work(nquad_e);
1172
1173 // if ((x = varcoeffs.find(VarCoeff[dir])) != varcoeffs.end())
1174 // {
1175 // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1176 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1177 // }
1178
1179 if (varcoeffs.find(StdRegions::eVarCoeffMF1x) != varcoeffs.end())
1180 {
1181 // MMF case
1182 Array<OneD, NekDouble> ncdotMF_e =
1183 GetnEdgecdotMF(dir, e, EdgeExp[e], normals, varcoeffs);
1184
1185 Vmath::Vmul(nquad_e, ncdotMF_e, 1, edgePhys, 1, edgePhys, 1);
1186 }
1187 else
1188 {
1189 Vmath::Vmul(nquad_e, normals[dir], 1, edgePhys, 1, edgePhys, 1);
1190 }
1191
1192 AddEdgeBoundaryInt(e, EdgeExp[e], edgePhys, outarray, varcoeffs);
1193 }
1194}
1195
1196void Expansion2D::AddNormTraceInt(
1197 const int dir, Array<OneD, ExpansionSharedPtr> &EdgeExp,
1198 Array<OneD, Array<OneD, NekDouble>> &edgeCoeffs,
1199 Array<OneD, NekDouble> &outarray)
1200{
1201 int e;
1202 int nquad_e;
1203 int nedges = GetNtraces();
1204
1205 for (e = 0; e < nedges; ++e)
1206 {
1207 nquad_e = EdgeExp[e]->GetNumPoints(0);
1208
1209 Array<OneD, NekDouble> edgePhys(nquad_e);
1210 const Array<OneD, const Array<OneD, NekDouble>> &normals =
1211 GetTraceNormal(e);
1212
1213 EdgeExp[e]->BwdTrans(edgeCoeffs[e], edgePhys);
1214
1215 Vmath::Vmul(nquad_e, normals[dir], 1, edgePhys, 1, edgePhys, 1);
1216
1217 AddEdgeBoundaryInt(e, EdgeExp[e], edgePhys, outarray);
1218 }
1219}
1220
1221/**
1222 * For a given edge add the \tilde{F}_1j contributions
1223 */
1224void Expansion2D::AddEdgeBoundaryInt(const int edge,
1225 ExpansionSharedPtr &EdgeExp,
1226 Array<OneD, NekDouble> &edgePhys,
1227 Array<OneD, NekDouble> &outarray,
1228 const StdRegions::VarCoeffMap &varcoeffs)
1229{
1230 int i;
1231 int order_e = EdgeExp->GetNcoeffs();
1232 int nquad_e = EdgeExp->GetNumPoints(0);
1233 Array<OneD, unsigned int> map;
1234 Array<OneD, int> sign;
1235 Array<OneD, NekDouble> coeff(order_e);
1236
1237 GetTraceToElementMap(edge, map, sign, v_GetTraceOrient(edge));
1238
1239 StdRegions::VarCoeffType VarCoeff[3] = {StdRegions::eVarCoeffD00,
1240 StdRegions::eVarCoeffD11,
1241 StdRegions::eVarCoeffD22};
1242 StdRegions::VarCoeffMap::const_iterator x;
1243
1244 /// @TODO Variable coeffs
1245 if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
1246 {
1247 Array<OneD, NekDouble> work(nquad_e);
1248 GetPhysEdgeVarCoeffsFromElement(edge, EdgeExp, x->second.GetValue(),
1249 work);
1250 Vmath::Vmul(nquad_e, work, 1, edgePhys, 1, edgePhys, 1);
1251 }
1252
1253 EdgeExp->IProductWRTBase(edgePhys, coeff);
1254
1255 // add data to out array
1256 for (i = 0; i < order_e; ++i)
1257 {
1258 outarray[map[i]] += sign[i] * coeff[i];
1259 }
1260}
1261
1262// This method assumes that data in EdgeExp is orientated according to
1263// elemental counter clockwise format AddHDGHelmholtzTraceTerms with
1264// directions
1265void Expansion2D::AddHDGHelmholtzTraceTerms(
1266 const NekDouble tau, const Array<OneD, const NekDouble> &inarray,
1267 Array<OneD, ExpansionSharedPtr> &EdgeExp,
1268 const StdRegions::VarCoeffMap &dirForcing, Array<OneD, NekDouble> &outarray)
1269{
1270 ASSERTL0(&inarray[0] != &outarray[0],
1271 "Input and output arrays use the same memory");
1272
1273 int e, cnt, order_e, nedges = GetNtraces();
1274 Array<OneD, const NekDouble> tmp;
1275
1276 cnt = 0;
1277
1278 for (e = 0; e < nedges; ++e)
1279 {
1280 order_e = EdgeExp[e]->GetNcoeffs();
1281 Array<OneD, NekDouble> edgeCoeffs(order_e);
1282 Array<OneD, NekDouble> edgePhys(EdgeExp[e]->GetTotPoints());
1283
1284 Vmath::Vcopy(order_e, tmp = inarray + cnt, 1, edgeCoeffs, 1);
1285 EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1286 AddHDGHelmholtzEdgeTerms(tau, e, EdgeExp, edgePhys, dirForcing,
1287 outarray);
1288
1289 cnt += order_e;
1290 }
1291}
1292
1293// evaluate additional terms in HDG edges. Not that this assumes that
1294// edges are unpacked into local cartesian order.
1295void Expansion2D::AddHDGHelmholtzEdgeTerms(
1296 const NekDouble tau, const int edge,
1297 Array<OneD, ExpansionSharedPtr> &EdgeExp, Array<OneD, NekDouble> &edgePhys,
1298 const StdRegions::VarCoeffMap &varcoeffs, Array<OneD, NekDouble> &outarray)
1299{
1300 bool mmf = (varcoeffs.find(StdRegions::eVarCoeffMF1x) != varcoeffs.end());
1301 int i, j, n;
1302 int nquad_e = EdgeExp[edge]->GetNumPoints(0);
1303 int order_e = EdgeExp[edge]->GetNcoeffs();
1304 int coordim = mmf ? 2 : GetCoordim();
1305 int ncoeffs = GetNcoeffs();
1306
1307 Array<OneD, NekDouble> inval(nquad_e);
1308 Array<OneD, NekDouble> outcoeff(order_e);
1309 Array<OneD, NekDouble> tmpcoeff(ncoeffs);
1310
1311 const Array<OneD, const Array<OneD, NekDouble>> &normals =
1312 GetTraceNormal(edge);
1313
1314 Array<OneD, unsigned int> emap;
1315 Array<OneD, int> sign;
1316
1317 StdRegions::Orientation edgedir = GetTraceOrient(edge);
1318
1319 DNekVec Coeffs(ncoeffs, outarray, eWrapper);
1320 DNekVec Tmpcoeff(ncoeffs, tmpcoeff, eWrapper);
1321
1322 GetTraceToElementMap(edge, emap, sign, edgedir);
1323
1324 StdRegions::MatrixType DerivType[3] = {StdRegions::eWeakDeriv0,
1325 StdRegions::eWeakDeriv1,
1326 StdRegions::eWeakDeriv2};
1327
1328 StdRegions::VarCoeffType VarCoeff[3] = {StdRegions::eVarCoeffD00,
1329 StdRegions::eVarCoeffD11,
1330 StdRegions::eVarCoeffD22};
1331
1332 StdRegions::VarCoeffMap::const_iterator x;
1333 /// @TODO: What direction to use here??
1334 if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
1335 {
1336 Array<OneD, NekDouble> work(nquad_e);
1337 GetPhysEdgeVarCoeffsFromElement(edge, EdgeExp[edge],
1338 x->second.GetValue(), work);
1339 Vmath::Vmul(nquad_e, work, 1, edgePhys, 1, edgePhys, 1);
1340 }
1341
1342 //================================================================
1343 // Add F = \tau <phi_i,in_phys>
1344 // Fill edge and take inner product
1345 EdgeExp[edge]->IProductWRTBase(edgePhys, outcoeff);
1346 // add data to out array
1347 for (i = 0; i < order_e; ++i)
1348 {
1349 outarray[emap[i]] += sign[i] * tau * outcoeff[i];
1350 }
1351 //================================================================
1352
1353 //===============================================================
1354 // Add -\sum_i D_i^T M^{-1} G_i + E_i M^{-1} G_i =
1355 // \sum_i D_i M^{-1} G_i term
1356
1357 // Two independent direction
1358 DNekScalMatSharedPtr invMass;
1359 for (n = 0; n < coordim; ++n)
1360 {
1361 if (mmf)
1362 {
1363 StdRegions::VarCoeffMap Weight;
1364 Weight[StdRegions::eVarCoeffMass] = GetMFMag(n, varcoeffs);
1365
1366 MatrixKey invMasskey(StdRegions::eInvMass, DetShapeType(), *this,
1367 StdRegions::NullConstFactorMap, Weight);
1368
1369 invMass = GetLocMatrix(invMasskey);
1370
1371 Array<OneD, NekDouble> ncdotMF_e =
1372 GetnEdgecdotMF(n, edge, EdgeExp[edge], normals, varcoeffs);
1373
1374 Vmath::Vmul(nquad_e, ncdotMF_e, 1, edgePhys, 1, inval, 1);
1375 }
1376 else
1377 {
1378 Vmath::Vmul(nquad_e, normals[n], 1, edgePhys, 1, inval, 1);
1379 invMass = GetLocMatrix(StdRegions::eInvMass);
1380 }
1381
1382 // Multiply by variable coefficient
1383 /// @TODO: Document this (probably not needed)
1384 // StdRegions::VarCoeffMap::const_iterator x;
1385 // if ((x = varcoeffs.find(VarCoeff[n])) !=
1386 // varcoeffs.end())
1387 // {
1388 // GetPhysEdgeVarCoeffsFromElement(edge,EdgeExp[edge],x->second,varcoeff_work);
1389 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[edge]->GetPhys(),1,EdgeExp[edge]->UpdatePhys(),1);
1390 // }
1391
1392 EdgeExp[edge]->IProductWRTBase(inval, outcoeff);
1393
1394 // M^{-1} G
1395 for (i = 0; i < ncoeffs; ++i)
1396 {
1397 tmpcoeff[i] = 0;
1398 for (j = 0; j < order_e; ++j)
1399 {
1400 tmpcoeff[i] += (*invMass)(i, emap[j]) * sign[j] * outcoeff[j];
1401 }
1402 }
1403
1404 if (mmf)
1405 {
1406 StdRegions::VarCoeffMap VarCoeffDirDeriv;
1407 VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
1408 GetMF(n, coordim, varcoeffs);
1409 VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
1410 GetMFDiv(n, varcoeffs);
1411
1412 MatrixKey Dmatkey(StdRegions::eWeakDirectionalDeriv, DetShapeType(),
1413 *this, StdRegions::NullConstFactorMap,
1414 VarCoeffDirDeriv);
1415
1416 DNekScalMat &Dmat = *GetLocMatrix(Dmatkey);
1417
1418 Coeffs = Coeffs + Dmat * Tmpcoeff;
1419 }
1420 else
1421 {
1422 if (varcoeffs.find(VarCoeff[n]) != varcoeffs.end())
1423 {
1424 MatrixKey mkey(DerivType[n], DetShapeType(), *this,
1425 StdRegions::NullConstFactorMap, varcoeffs);
1426
1427 DNekScalMat &Dmat = *GetLocMatrix(mkey);
1428 Coeffs = Coeffs + Dmat * Tmpcoeff;
1429 }
1430 else
1431 {
1432 DNekScalMat &Dmat = *GetLocMatrix(DerivType[n]);
1433 Coeffs = Coeffs + Dmat * Tmpcoeff;
1434 }
1435 }
1436 }
1437}
1438
1439/**
1440 * Extracts the variable coefficients along an edge
1441 */
1442void Expansion2D::GetPhysEdgeVarCoeffsFromElement(
1443 const int edge, ExpansionSharedPtr &EdgeExp,
1444 const Array<OneD, const NekDouble> &varcoeff,
1445 Array<OneD, NekDouble> &outarray)
1446{
1447 Array<OneD, NekDouble> tmp(GetNcoeffs());
1448 Array<OneD, NekDouble> edgetmp(EdgeExp->GetNcoeffs());
1449
1450 // FwdTrans varcoeffs
1451 FwdTrans(varcoeff, tmp);
1452
1453 // Map to edge
1454 Array<OneD, unsigned int> emap;
1455 Array<OneD, int> sign;
1456 StdRegions::Orientation edgedir = GetTraceOrient(edge);
1457 GetTraceToElementMap(edge, emap, sign, edgedir);
1458
1459 for (unsigned int i = 0; i < EdgeExp->GetNcoeffs(); ++i)
1460 {
1461 edgetmp[i] = tmp[emap[i]];
1462 }
1463
1464 // BwdTrans
1465 EdgeExp->BwdTrans(edgetmp, outarray);
1466}
1467
1468/**
1469 * Computes matrices needed for the HDG formulation. References to
1470 * equations relate to the following paper:
1471 * R. M. Kirby, S. J. Sherwin, B. Cockburn, To CG or to HDG: A
1472 * Comparative Study, J. Sci. Comp P1-30
1473 * DOI 10.1007/s10915-011-9501-7
1474 */
1475DNekMatSharedPtr Expansion2D::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
1476{
1477 DNekMatSharedPtr returnval;
1478
1479 switch (mkey.GetMatrixType())
1480 {
1481 // (Z^e)^{-1} (Eqn. 33, P22)
1482 case StdRegions::eHybridDGHelmholtz:
1483 {
1484 ASSERTL1(IsBoundaryInteriorExpansion(),
1485 "HybridDGHelmholtz matrix not set up "
1486 "for non boundary-interior expansions");
1487
1488 int i, j, k;
1489 NekDouble lambdaval =
1490 mkey.GetConstFactor(StdRegions::eFactorLambda);
1491 NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
1492 int ncoeffs = GetNcoeffs();
1493 int nedges = GetNtraces();
1494 int shapedim = 2;
1495 const StdRegions::VarCoeffMap &varcoeffs = mkey.GetVarCoeffs();
1496 bool mmf =
1497 (varcoeffs.find(StdRegions::eVarCoeffMF1x) != varcoeffs.end());
1498
1499 Array<OneD, unsigned int> emap;
1500 Array<OneD, int> sign;
1501 StdRegions::Orientation edgedir = StdRegions::eForwards;
1502 ExpansionSharedPtr EdgeExp;
1503
1504 int order_e, coordim = GetCoordim();
1505 DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
1506 StdRegions::MatrixType DerivType[3] = {StdRegions::eWeakDeriv0,
1507 StdRegions::eWeakDeriv1,
1508 StdRegions::eWeakDeriv2};
1509 DNekMat LocMat(ncoeffs, ncoeffs);
1510
1511 returnval =
1512 MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs, ncoeffs);
1513 DNekMat &Mat = *returnval;
1514 Vmath::Zero(ncoeffs * ncoeffs, Mat.GetPtr(), 1);
1515
1516 StdRegions::VarCoeffType Coeffs[3] = {StdRegions::eVarCoeffD00,
1517 StdRegions::eVarCoeffD11,
1518 StdRegions::eVarCoeffD22};
1519
1520 StdRegions::VarCoeffMap::const_iterator x;
1521
1522 for (i = 0; i < coordim; ++i)
1523 {
1524 if (mmf)
1525 {
1526 if (i < shapedim)
1527 {
1528 StdRegions::VarCoeffMap VarCoeffDirDeriv;
1529 VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
1530 GetMF(i, shapedim, varcoeffs);
1531 VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
1532 GetMFDiv(i, varcoeffs);
1533
1534 MatrixKey Dmatkey(StdRegions::eWeakDirectionalDeriv,
1535 DetShapeType(), *this,
1536 StdRegions::NullConstFactorMap,
1537 VarCoeffDirDeriv);
1538
1539 DNekScalMat &Dmat = *GetLocMatrix(Dmatkey);
1540
1541 StdRegions::VarCoeffMap Weight;
1542 Weight[StdRegions::eVarCoeffMass] =
1543 GetMFMag(i, mkey.GetVarCoeffs());
1544
1545 MatrixKey invMasskey(
1546 StdRegions::eInvMass, DetShapeType(), *this,
1547 StdRegions::NullConstFactorMap, Weight);
1548
1549 DNekScalMat &invMass = *GetLocMatrix(invMasskey);
1550
1551 Mat = Mat + Dmat * invMass * Transpose(Dmat);
1552 }
1553 }
1554 else if (mkey.HasVarCoeff(Coeffs[i]))
1555 {
1556 MatrixKey DmatkeyL(DerivType[i], DetShapeType(), *this,
1557 StdRegions::NullConstFactorMap,
1558 mkey.GetVarCoeffAsMap(Coeffs[i]));
1559
1560 MatrixKey DmatkeyR(DerivType[i], DetShapeType(), *this);
1561
1562 DNekScalMat &DmatL = *GetLocMatrix(DmatkeyL);
1563 DNekScalMat &DmatR = *GetLocMatrix(DmatkeyR);
1564 Mat = Mat + DmatL * invMass * Transpose(DmatR);
1565 }
1566 else
1567 {
1568 DNekScalMat &Dmat = *GetLocMatrix(DerivType[i]);
1569 Mat = Mat + Dmat * invMass * Transpose(Dmat);
1570 }
1571 }
1572
1573 // Add Mass Matrix Contribution for Helmholtz problem
1574 DNekScalMat &Mass = *GetLocMatrix(StdRegions::eMass);
1575 Mat = Mat + lambdaval * Mass;
1576
1577 // Add tau*E_l using elemental mass matrices on each edge
1578 for (i = 0; i < nedges; ++i)
1579 {
1580 EdgeExp = GetTraceExp(i);
1581 order_e = EdgeExp->GetNcoeffs();
1582
1583 int nq = EdgeExp->GetNumPoints(0);
1584 GetTraceToElementMap(i, emap, sign, edgedir);
1585
1586 // @TODO: Document
1587 StdRegions::VarCoeffMap edgeVarCoeffs;
1588 if (mkey.HasVarCoeff(StdRegions::eVarCoeffD00))
1589 {
1590 Array<OneD, NekDouble> mu(nq);
1591 GetPhysEdgeVarCoeffsFromElement(
1592 i, EdgeExp, mkey.GetVarCoeff(StdRegions::eVarCoeffD00),
1593 mu);
1594 edgeVarCoeffs[StdRegions::eVarCoeffMass] = mu;
1595 }
1596 DNekScalMat &eMass = *EdgeExp->GetLocMatrix(
1597 StdRegions::eMass, StdRegions::NullConstFactorMap,
1598 edgeVarCoeffs);
1599 // DNekScalMat &eMass =
1600 // *EdgeExp->GetLocMatrix(StdRegions::eMass);
1601
1602 for (j = 0; j < order_e; ++j)
1603 {
1604 for (k = 0; k < order_e; ++k)
1605 {
1606 Mat(emap[j], emap[k]) =
1607 Mat(emap[j], emap[k]) +
1608 tau * sign[j] * sign[k] * eMass(j, k);
1609 }
1610 }
1611 }
1612 }
1613 break;
1614 // U^e (P22)
1615 case StdRegions::eHybridDGLamToU:
1616 {
1617 int i, j, k;
1618 int nbndry = NumDGBndryCoeffs();
1619 int ncoeffs = GetNcoeffs();
1620 int nedges = GetNtraces();
1621 NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
1622
1623 Array<OneD, NekDouble> lambda(nbndry);
1624 DNekVec Lambda(nbndry, lambda, eWrapper);
1625 Array<OneD, NekDouble> ulam(ncoeffs);
1626 DNekVec Ulam(ncoeffs, ulam, eWrapper);
1627 Array<OneD, NekDouble> f(ncoeffs);
1628 DNekVec F(ncoeffs, f, eWrapper);
1629
1630 Array<OneD, ExpansionSharedPtr> EdgeExp(nedges);
1631 // declare matrix space
1632 returnval =
1633 MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs, nbndry);
1634 DNekMat &Umat = *returnval;
1635
1636 // Z^e matrix
1637 MatrixKey newkey(StdRegions::eInvHybridDGHelmholtz, DetShapeType(),
1638 *this, mkey.GetConstFactors(),
1639 mkey.GetVarCoeffs());
1640 DNekScalMat &invHmat = *GetLocMatrix(newkey);
1641
1642 Array<OneD, unsigned int> emap;
1643 Array<OneD, int> sign;
1644
1645 for (i = 0; i < nedges; ++i)
1646 {
1647 EdgeExp[i] = GetTraceExp(i);
1648 }
1649
1650 // for each degree of freedom of the lambda space
1651 // calculate Umat entry
1652 // Generate Lambda to U_lambda matrix
1653 for (j = 0; j < nbndry; ++j)
1654 {
1655 // standard basis vectors e_j
1656 Vmath::Zero(nbndry, &lambda[0], 1);
1657 Vmath::Zero(ncoeffs, &f[0], 1);
1658 lambda[j] = 1.0;
1659
1660 SetTraceToGeomOrientation(EdgeExp, lambda);
1661
1662 // Compute F = [I D_1 M^{-1} D_2 M^{-1}] C e_j
1663 AddHDGHelmholtzTraceTerms(tau, lambda, EdgeExp,
1664 mkey.GetVarCoeffs(), f);
1665
1666 // Compute U^e_j
1667 Ulam = invHmat * F; // generate Ulam from lambda
1668
1669 // fill column of matrix
1670 for (k = 0; k < ncoeffs; ++k)
1671 {
1672 Umat(k, j) = Ulam[k];
1673 }
1674 }
1675 }
1676 break;
1677 // Q_0, Q_1, Q_2 matrices (P23)
1678 // Each are a product of a row of Eqn 32 with the C matrix.
1679 // Rather than explicitly computing all of Eqn 32, we note each
1680 // row is almost a multiple of U^e, so use that as our starting
1681 // point.
1682 case StdRegions::eHybridDGLamToQ0:
1683 case StdRegions::eHybridDGLamToQ1:
1684 case StdRegions::eHybridDGLamToQ2:
1685 {
1686 int i = 0;
1687 int j = 0;
1688 int k = 0;
1689 int dir = 0;
1690 int nbndry = NumDGBndryCoeffs();
1691 int ncoeffs = GetNcoeffs();
1692 int nedges = GetNtraces();
1693 int shapedim = 2;
1694
1695 Array<OneD, NekDouble> lambda(nbndry);
1696 DNekVec Lambda(nbndry, lambda, eWrapper);
1697 Array<OneD, ExpansionSharedPtr> EdgeExp(nedges);
1698
1699 Array<OneD, NekDouble> ulam(ncoeffs);
1700 DNekVec Ulam(ncoeffs, ulam, eWrapper);
1701 Array<OneD, NekDouble> f(ncoeffs);
1702 DNekVec F(ncoeffs, f, eWrapper);
1703
1704 // declare matrix space
1705 returnval =
1706 MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs, nbndry);
1707 DNekMat &Qmat = *returnval;
1708
1709 // Lambda to U matrix
1710 MatrixKey lamToUkey(StdRegions::eHybridDGLamToU, DetShapeType(),
1711 *this, mkey.GetConstFactors(),
1712 mkey.GetVarCoeffs());
1713 DNekScalMat &lamToU = *GetLocMatrix(lamToUkey);
1714
1715 // Inverse mass matrix
1716 DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
1717
1718 for (i = 0; i < nedges; ++i)
1719 {
1720 EdgeExp[i] = GetTraceExp(i);
1721 }
1722
1723 // Weak Derivative matrix
1724 DNekScalMatSharedPtr Dmat;
1725 switch (mkey.GetMatrixType())
1726 {
1727 case StdRegions::eHybridDGLamToQ0:
1728 dir = 0;
1729 Dmat = GetLocMatrix(StdRegions::eWeakDeriv0);
1730 break;
1731 case StdRegions::eHybridDGLamToQ1:
1732 dir = 1;
1733 Dmat = GetLocMatrix(StdRegions::eWeakDeriv1);
1734 break;
1735 case StdRegions::eHybridDGLamToQ2:
1736 dir = 2;
1737 Dmat = GetLocMatrix(StdRegions::eWeakDeriv2);
1738 break;
1739 default:
1740 ASSERTL0(false, "Direction not known");
1741 break;
1742 }
1743
1744 const StdRegions::VarCoeffMap &varcoeffs = mkey.GetVarCoeffs();
1745 if (varcoeffs.find(StdRegions::eVarCoeffMF1x) != varcoeffs.end())
1746 {
1747 StdRegions::VarCoeffMap VarCoeffDirDeriv;
1748 VarCoeffDirDeriv[StdRegions::eVarCoeffMF] =
1749 GetMF(dir, shapedim, varcoeffs);
1750 VarCoeffDirDeriv[StdRegions::eVarCoeffMFDiv] =
1751 GetMFDiv(dir, varcoeffs);
1752
1753 MatrixKey Dmatkey(
1754 StdRegions::eWeakDirectionalDeriv, DetShapeType(), *this,
1755 StdRegions::NullConstFactorMap, VarCoeffDirDeriv);
1756
1757 Dmat = GetLocMatrix(Dmatkey);
1758
1759 StdRegions::VarCoeffMap Weight;
1760 Weight[StdRegions::eVarCoeffMass] =
1761 GetMFMag(dir, mkey.GetVarCoeffs());
1762
1763 MatrixKey invMasskey(StdRegions::eInvMass, DetShapeType(),
1764 *this, StdRegions::NullConstFactorMap,
1765 Weight);
1766
1767 invMass = *GetLocMatrix(invMasskey);
1768 }
1769 else
1770 {
1771 StdRegions::MatrixType DerivType[3] = {StdRegions::eWeakDeriv0,
1772 StdRegions::eWeakDeriv1,
1773 StdRegions::eWeakDeriv2};
1774
1775 Dmat = GetLocMatrix(DerivType[dir]);
1776
1777 MatrixKey invMasskey(StdRegions::eInvMass, DetShapeType(),
1778 *this);
1779 invMass = *GetLocMatrix(invMasskey);
1780 }
1781
1782 // for each degree of freedom of the lambda space
1783 // calculate Qmat entry
1784 // Generate Lambda to Q_lambda matrix
1785 for (j = 0; j < nbndry; ++j)
1786 {
1787 Vmath::Zero(nbndry, &lambda[0], 1);
1788 lambda[j] = 1.0;
1789
1790 // for lambda[j] = 1 this is the solution to ulam
1791 for (k = 0; k < ncoeffs; ++k)
1792 {
1793 Ulam[k] = lamToU(k, j);
1794 }
1795
1796 // -D^T ulam
1797 Vmath::Neg(ncoeffs, &ulam[0], 1);
1798 F = Transpose(*Dmat) * Ulam;
1799
1800 SetTraceToGeomOrientation(EdgeExp, lambda);
1801
1802 // Add the C terms resulting from the I's on the
1803 // diagonals of Eqn 32
1804 AddNormTraceInt(dir, lambda, EdgeExp, f, mkey.GetVarCoeffs());
1805
1806 // finally multiply by inverse mass matrix
1807 Ulam = invMass * F;
1808
1809 // fill column of matrix (Qmat is in column major format)
1810 Vmath::Vcopy(ncoeffs, &ulam[0], 1,
1811 &(Qmat.GetPtr())[0] + j * ncoeffs, 1);
1812 }
1813 }
1814 break;
1815 // Matrix K (P23)
1816 case StdRegions::eHybridDGHelmBndLam:
1817 {
1818 int i, j, e, cnt;
1819 int order_e, nquad_e;
1820 int nbndry = NumDGBndryCoeffs();
1821 int coordim = GetCoordim();
1822 int nedges = GetNtraces();
1823 NekDouble tau = mkey.GetConstFactor(StdRegions::eFactorTau);
1824 StdRegions::VarCoeffMap::const_iterator x;
1825 const StdRegions::VarCoeffMap &varcoeffs = mkey.GetVarCoeffs();
1826 bool mmf =
1827 (varcoeffs.find(StdRegions::eVarCoeffMF1x) != varcoeffs.end());
1828
1829 Array<OneD, NekDouble> work, varcoeff_work;
1830 Array<OneD, const Array<OneD, NekDouble>> normals;
1831 Array<OneD, ExpansionSharedPtr> EdgeExp(nedges);
1832 Array<OneD, NekDouble> lam(nbndry);
1833
1834 Array<OneD, unsigned int> emap;
1835 Array<OneD, int> sign;
1836 StdRegions::Orientation edgedir;
1837
1838 // declare matrix space
1839 returnval =
1840 MemoryManager<DNekMat>::AllocateSharedPtr(nbndry, nbndry);
1841 DNekMat &BndMat = *returnval;
1842
1843 DNekScalMatSharedPtr LamToQ[3];
1844
1845 // Matrix to map Lambda to U
1846 MatrixKey LamToUkey(StdRegions::eHybridDGLamToU, DetShapeType(),
1847 *this, mkey.GetConstFactors(),
1848 mkey.GetVarCoeffs());
1849 DNekScalMat &LamToU = *GetLocMatrix(LamToUkey);
1850
1851 // Matrix to map Lambda to Q0
1852 MatrixKey LamToQ0key(StdRegions::eHybridDGLamToQ0, DetShapeType(),
1853 *this, mkey.GetConstFactors(),
1854 mkey.GetVarCoeffs());
1855 LamToQ[0] = GetLocMatrix(LamToQ0key);
1856
1857 // Matrix to map Lambda to Q1
1858 MatrixKey LamToQ1key(StdRegions::eHybridDGLamToQ1, DetShapeType(),
1859 *this, mkey.GetConstFactors(),
1860 mkey.GetVarCoeffs());
1861 LamToQ[1] = GetLocMatrix(LamToQ1key);
1862
1863 // Matrix to map Lambda to Q2 for 3D coordinates
1864 if (coordim == 3)
1865 {
1866 MatrixKey LamToQ2key(
1867 StdRegions::eHybridDGLamToQ2, DetShapeType(), *this,
1868 mkey.GetConstFactors(), mkey.GetVarCoeffs());
1869 LamToQ[2] = GetLocMatrix(LamToQ2key);
1870 }
1871
1872 // Set up edge segment expansions from local geom info
1873 for (i = 0; i < nedges; ++i)
1874 {
1875 EdgeExp[i] = GetTraceExp(i);
1876 }
1877
1878 // Set up matrix derived from <mu, Q_lam.n - \tau (U_lam - Lam) >
1879 for (i = 0; i < nbndry; ++i)
1880 {
1881 cnt = 0;
1882
1883 Vmath::Zero(nbndry, lam, 1);
1884 lam[i] = 1.0;
1885 SetTraceToGeomOrientation(EdgeExp, lam);
1886
1887 for (e = 0; e < nedges; ++e)
1888 {
1889 order_e = EdgeExp[e]->GetNcoeffs();
1890 nquad_e = EdgeExp[e]->GetNumPoints(0);
1891
1892 normals = GetTraceNormal(e);
1893 edgedir = GetTraceOrient(e);
1894
1895 work = Array<OneD, NekDouble>(nquad_e);
1896 varcoeff_work = Array<OneD, NekDouble>(nquad_e);
1897
1898 GetTraceToElementMap(e, emap, sign, edgedir);
1899
1900 StdRegions::VarCoeffType VarCoeff[3] = {
1901 StdRegions::eVarCoeffD00, StdRegions::eVarCoeffD11,
1902 StdRegions::eVarCoeffD22};
1903
1904 // Q0 * n0 (BQ_0 terms)
1905 Array<OneD, NekDouble> edgeCoeffs(order_e);
1906 Array<OneD, NekDouble> edgePhys(nquad_e);
1907 for (j = 0; j < order_e; ++j)
1908 {
1909 edgeCoeffs[j] = sign[j] * (*LamToQ[0])(emap[j], i);
1910 }
1911
1912 EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1913 // @TODO Var coeffs
1914 // Multiply by variable coefficient
1915 // if ((x =
1916 // varcoeffs.find(VarCoeff[0]))
1917 // != varcoeffs.end())
1918 // {
1919 // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1920 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1921 // }
1922 if (mmf)
1923 {
1924 Array<OneD, NekDouble> ncdotMF = GetnEdgecdotMF(
1925 0, e, EdgeExp[e], normals, varcoeffs);
1926 Vmath::Vmul(nquad_e, ncdotMF, 1, edgePhys, 1, work, 1);
1927 }
1928 else
1929 {
1930 Vmath::Vmul(nquad_e, normals[0], 1, edgePhys, 1, work,
1931 1);
1932 }
1933
1934 // Q1 * n1 (BQ_1 terms)
1935 for (j = 0; j < order_e; ++j)
1936 {
1937 edgeCoeffs[j] = sign[j] * (*LamToQ[1])(emap[j], i);
1938 }
1939
1940 EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1941
1942 // @TODO var coeffs
1943 // Multiply by variable coefficients
1944 // if ((x =
1945 // varcoeffs.find(VarCoeff[1]))
1946 // != varcoeffs.end())
1947 // {
1948 // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1949 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1950 // }
1951
1952 if (mmf)
1953 {
1954 Array<OneD, NekDouble> ncdotMF = GetnEdgecdotMF(
1955 1, e, EdgeExp[e], normals, varcoeffs);
1956 Vmath::Vvtvp(nquad_e, ncdotMF, 1, edgePhys, 1, work, 1,
1957 work, 1);
1958 }
1959 else
1960 {
1961 Vmath::Vvtvp(nquad_e, normals[1], 1, edgePhys, 1, work,
1962 1, work, 1);
1963 }
1964
1965 // Q2 * n2 (BQ_2 terms)
1966 if (coordim == 3)
1967 {
1968 for (j = 0; j < order_e; ++j)
1969 {
1970 edgeCoeffs[j] = sign[j] * (*LamToQ[2])(emap[j], i);
1971 }
1972
1973 EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1974 // @TODO var coeffs
1975 // Multiply by variable coefficients
1976 // if ((x =
1977 // varcoeffs.find(VarCoeff[2]))
1978 // != varcoeffs.end())
1979 // {
1980 // GetPhysEdgeVarCoeffsFromElement(e,EdgeExp[e],x->second,varcoeff_work);
1981 // Vmath::Vmul(nquad_e,varcoeff_work,1,EdgeExp[e]->GetPhys(),1,EdgeExp[e]->UpdatePhys(),1);
1982 // }
1983
1984 Vmath::Vvtvp(nquad_e, normals[2], 1, edgePhys, 1, work,
1985 1, work, 1);
1986 }
1987
1988 // - tau (ulam - lam)
1989 // Corresponds to the G and BU terms.
1990 for (j = 0; j < order_e; ++j)
1991 {
1992 edgeCoeffs[j] =
1993 sign[j] * LamToU(emap[j], i) - lam[cnt + j];
1994 }
1995
1996 EdgeExp[e]->BwdTrans(edgeCoeffs, edgePhys);
1997
1998 // Multiply by variable coefficients
1999 if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
2000 {
2001 GetPhysEdgeVarCoeffsFromElement(
2002 e, EdgeExp[e], x->second.GetValue(), varcoeff_work);
2003 Vmath::Vmul(nquad_e, varcoeff_work, 1, edgePhys, 1,
2004 edgePhys, 1);
2005 }
2006
2007 Vmath::Svtvp(nquad_e, -tau, edgePhys, 1, work, 1, work, 1);
2008 /// TODO: Add variable coeffs
2009 EdgeExp[e]->IProductWRTBase(work, edgeCoeffs);
2010
2011 EdgeExp[e]->SetCoeffsToOrientation(edgedir, edgeCoeffs,
2012 edgeCoeffs);
2013
2014 for (j = 0; j < order_e; ++j)
2015 {
2016 BndMat(cnt + j, i) = edgeCoeffs[j];
2017 }
2018
2019 cnt += order_e;
2020 }
2021 }
2022 }
2023 break;
2024 // HDG postprocessing
2025 case StdRegions::eInvLaplacianWithUnityMean:
2026 {
2027 MatrixKey lapkey(StdRegions::eLaplacian, DetShapeType(), *this,
2028 mkey.GetConstFactors(), mkey.GetVarCoeffs());
2029 DNekScalMat &LapMat = *GetLocMatrix(lapkey);
2030
2031 returnval = MemoryManager<DNekMat>::AllocateSharedPtr(
2032 LapMat.GetRows(), LapMat.GetColumns());
2033 DNekMatSharedPtr lmat = returnval;
2034
2035 (*lmat) = LapMat;
2036
2037 // replace first column with inner product wrt 1
2038 int nq = GetTotPoints();
2039 Array<OneD, NekDouble> tmp(nq);
2040 Array<OneD, NekDouble> outarray(m_ncoeffs);
2041 Vmath::Fill(nq, 1.0, tmp, 1);
2042 IProductWRTBase(tmp, outarray);
2043
2044 Vmath::Vcopy(m_ncoeffs, &outarray[0], 1, &(lmat->GetPtr())[0], 1);
2045 lmat->Invert();
2046 }
2047 break;
2048 case StdRegions::eNormDerivOnTrace:
2049 {
2050 int ntraces = GetNtraces();
2051 int ncoords = GetCoordim();
2052 int nphys = GetTotPoints();
2053 Array<OneD, const Array<OneD, NekDouble>> normals;
2054 Array<OneD, NekDouble> phys(nphys);
2055 returnval =
2056 MemoryManager<DNekMat>::AllocateSharedPtr(m_ncoeffs, m_ncoeffs);
2057 DNekMat &Mat = *returnval;
2058 Vmath::Zero(m_ncoeffs * m_ncoeffs, Mat.GetPtr(), 1);
2059
2060 Array<OneD, Array<OneD, NekDouble>> Deriv(3, NullNekDouble1DArray);
2061
2062 for (int d = 0; d < ncoords; ++d)
2063 {
2064 Deriv[d] = Array<OneD, NekDouble>(nphys);
2065 }
2066
2067 Array<OneD, int> tracepts(ntraces);
2068 Array<OneD, ExpansionSharedPtr> traceExp(ntraces);
2069 int maxtpts = 0;
2070 for (int t = 0; t < ntraces; ++t)
2071 {
2072 traceExp[t] = GetTraceExp(t);
2073 tracepts[t] = traceExp[t]->GetTotPoints();
2074 maxtpts = (maxtpts > tracepts[t]) ? maxtpts : tracepts[t];
2075 }
2076
2077 Array<OneD, NekDouble> val(maxtpts), tmp, tmp1;
2078
2079 Array<OneD, Array<OneD, NekDouble>> dphidn(ntraces);
2080 for (int t = 0; t < ntraces; ++t)
2081 {
2082 dphidn[t] =
2083 Array<OneD, NekDouble>(m_ncoeffs * tracepts[t], 0.0);
2084 }
2085
2086 for (int i = 0; i < m_ncoeffs; ++i)
2087 {
2088 FillMode(i, phys);
2089 v_PhysDeriv(phys, Deriv[0], Deriv[1], Deriv[2]);
2090
2091 for (int t = 0; t < ntraces; ++t)
2092 {
2093 const NormalVector norm = GetTraceNormal(t);
2094
2095 LibUtilities::BasisKey fromkey = GetTraceBasisKey(t);
2096 LibUtilities::BasisKey tokey =
2097 traceExp[t]->GetBasis(0)->GetBasisKey();
2098 bool DoInterp = (fromkey != tokey);
2099
2100 Array<OneD, NekDouble> n(tracepts[t]);
2101 ;
2102 for (int d = 0; d < ncoords; ++d)
2103 {
2104 // if variable p may need to interpolate
2105 if (DoInterp)
2106 {
2107 LibUtilities::Interp1D(fromkey, norm[d], tokey, n);
2108 }
2109 else
2110 {
2111 n = norm[d];
2112 }
2113
2114 GetLocTracePhysVals(t, traceExp[t], Deriv[d], val);
2115
2116 Vmath::Vvtvp(tracepts[t], n, 1, val, 1,
2117 tmp = dphidn[t] + i * tracepts[t], 1,
2118 tmp1 = dphidn[t] + i * tracepts[t], 1);
2119 }
2120 }
2121 }
2122
2123 for (int t = 0; t < ntraces; ++t)
2124 {
2125 int nt = tracepts[t];
2126 NekDouble h, p;
2127 TraceNormLen(t, h, p);
2128
2129 // scaling of trace
2130 ASSERTL1(mkey.HasVarFactors(StdRegions::eFactorGJPTraceWeight),
2131 "Cannot find TraceWeights in key");
2132 NekDouble scale =
2133 mkey.GetVarFactors(StdRegions::eFactorGJPTraceWeight)[t];
2134
2135 for (int i = 0; i < m_ncoeffs; ++i)
2136 {
2137 for (int j = i; j < m_ncoeffs; ++j)
2138 {
2139 Vmath::Vmul(nt, dphidn[t] + i * nt, 1,
2140 dphidn[t] + j * nt, 1, val, 1);
2141 Mat(i, j) =
2142 Mat(i, j) + scale * traceExp[t]->Integral(val);
2143 }
2144 }
2145 }
2146
2147 // fill in symmetric components.
2148 for (int i = 0; i < m_ncoeffs; ++i)
2149 {
2150 for (int j = 0; j < i; ++j)
2151 {
2152 Mat(i, j) = Mat(j, i);
2153 }
2154 }
2155 }
2156 break;
2157 default:
2158 ASSERTL0(false,
2159 "This matrix type cannot be generated from this class");
2160 break;
2161 }
2162
2163 return returnval;
2164}
2165
2166// Evaluate Coefficients of weak deriviative in the direction dir
2167// given the input coefficicents incoeffs and the imposed
2168// boundary values in EdgeExp (which will have its phys space updated);
2169void Expansion2D::v_DGDeriv(int dir,
2170 const Array<OneD, const NekDouble> &incoeffs,
2171 Array<OneD, ExpansionSharedPtr> &EdgeExp,
2172 Array<OneD, Array<OneD, NekDouble>> &edgeCoeffs,
2173 Array<OneD, NekDouble> &out_d)
2174{
2175 StdRegions::MatrixType DerivType[3] = {StdRegions::eWeakDeriv0,
2176 StdRegions::eWeakDeriv1,
2177 StdRegions::eWeakDeriv2};
2178
2179 int ncoeffs = GetNcoeffs();
2180
2181 DNekScalMat &InvMass = *GetLocMatrix(StdRegions::eInvMass);
2182 DNekScalMat &Dmat = *GetLocMatrix(DerivType[dir]);
2183
2184 Array<OneD, NekDouble> coeffs = incoeffs;
2185 DNekVec Coeffs(ncoeffs, coeffs, eWrapper);
2186
2187 Coeffs = Transpose(Dmat) * Coeffs;
2188 Vmath::Neg(ncoeffs, coeffs, 1);
2189
2190 // Add the boundary integral including the relevant part of
2191 // the normal
2192 AddNormTraceInt(dir, EdgeExp, edgeCoeffs, coeffs);
2193
2194 DNekVec Out_d(ncoeffs, out_d, eWrapper);
2195
2196 Out_d = InvMass * Coeffs;
2197}
2198
2205
2206void Expansion2D::v_AddRobinMassMatrix(
2207 const int edge, const Array<OneD, const NekDouble> &primCoeffs,
2208 DNekMatSharedPtr &inoutmat)
2209{
2210 ASSERTL1(IsBoundaryInteriorExpansion(),
2211 "Not set up for non boundary-interior expansions");
2212 ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
2213 "Assuming that input matrix was square");
2214 int i, j;
2215 int id1, id2;
2216 ExpansionSharedPtr edgeExp = m_traceExp[edge].lock();
2217 int order_e = edgeExp->GetNcoeffs();
2218
2219 Array<OneD, unsigned int> map;
2220 Array<OneD, int> sign;
2221
2222 StdRegions::VarCoeffMap varcoeffs;
2223 varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
2224
2225 LocalRegions::MatrixKey mkey(StdRegions::eMass, LibUtilities::eSegment,
2226 *edgeExp, StdRegions::NullConstFactorMap,
2227 varcoeffs);
2228 DNekScalMat &edgemat = *edgeExp->GetLocMatrix(mkey);
2229
2230 // Now need to identify a map which takes the local edge
2231 // mass matrix to the matrix stored in inoutmat;
2232 // This can currently be deduced from the size of the matrix
2233
2234 // - if inoutmat.m_rows() == v_NCoeffs() it is a full
2235 // matrix system
2236
2237 // - if inoutmat.m_rows() == v_NumBndCoeffs() it is a
2238 // boundary CG system
2239
2240 // - if inoutmat.m_rows() == v_NumDGBndCoeffs() it is a
2241 // trace DG system
2242 int rows = inoutmat->GetRows();
2243
2244 if (rows == GetNcoeffs())
2245 {
2246 GetTraceToElementMap(edge, map, sign, v_GetTraceOrient(edge));
2247 }
2248 else if (rows == NumBndryCoeffs())
2249 {
2250 int nbndry = NumBndryCoeffs();
2251 Array<OneD, unsigned int> bmap(nbndry);
2252
2253 GetTraceToElementMap(edge, map, sign, v_GetTraceOrient(edge));
2254
2255 GetBoundaryMap(bmap);
2256
2257 for (i = 0; i < order_e; ++i)
2258 {
2259 for (j = 0; j < nbndry; ++j)
2260 {
2261 if (map[i] == bmap[j])
2262 {
2263 map[i] = j;
2264 break;
2265 }
2266 }
2267 ASSERTL1(j != nbndry, "Did not find number in map");
2268 }
2269 }
2270 else if (rows == NumDGBndryCoeffs())
2271 {
2272 // possibly this should be a separate method
2273 int cnt = 0;
2274 map = Array<OneD, unsigned int>(order_e);
2275 sign = Array<OneD, int>(order_e, 1);
2276
2277 for (i = 0; i < edge; ++i)
2278 {
2279 cnt += GetTraceNcoeffs(i);
2280 }
2281
2282 for (i = 0; i < order_e; ++i)
2283 {
2284 map[i] = cnt++;
2285 }
2286 // check for mapping reversal
2287 if (GetTraceOrient(edge) == StdRegions::eBackwards)
2288 {
2289 switch (edgeExp->GetBasis(0)->GetBasisType())
2290 {
2291 case LibUtilities::eGauss_Lagrange:
2292 reverse(map.data(), map.data() + order_e);
2293 break;
2294 case LibUtilities::eGLL_Lagrange:
2295 reverse(map.data(), map.data() + order_e);
2296 break;
2297 case LibUtilities::eModified_A:
2298 {
2299 swap(map[0], map[1]);
2300 for (i = 3; i < order_e; i += 2)
2301 {
2302 sign[i] = -1;
2303 }
2304 }
2305 break;
2306 default:
2307 ASSERTL0(false,
2308 "Edge boundary type not valid for this method");
2309 }
2310 }
2311 }
2312 else
2313 {
2314 ASSERTL0(false, "Could not identify matrix type from dimension");
2315 }
2316
2317 for (i = 0; i < order_e; ++i)
2318 {
2319 id1 = map[i];
2320 for (j = 0; j < order_e; ++j)
2321 {
2322 id2 = map[j];
2323 (*inoutmat)(id1, id2) += edgemat(i, j) * sign[i] * sign[j];
2324 }
2325 }
2326}
2327
2328/**
2329 * Given an edge and vector of element coefficients:
2330 * - maps those elemental coefficients corresponding to the edge into
2331 * an edge-vector.
2332 * - resets the element coefficients
2333 * - multiplies the edge vector by the edge mass matrix
2334 * - maps the edge coefficients back onto the elemental coefficients
2335 */
2336void Expansion2D::v_AddRobinTraceContribution(
2337 const int edgeid, const Array<OneD, const NekDouble> &primCoeffs,
2338 const Array<OneD, NekDouble> &incoeffs, Array<OneD, NekDouble> &coeffs)
2339{
2340 ASSERTL1(IsBoundaryInteriorExpansion(),
2341 "Not set up for non boundary-interior expansions");
2342 int i;
2343 ExpansionSharedPtr edgeExp = m_traceExp[edgeid].lock();
2344 int order_e = edgeExp->GetNcoeffs();
2345
2346 Array<OneD, unsigned int> map;
2347 Array<OneD, int> sign;
2348
2349 StdRegions::VarCoeffMap varcoeffs;
2350 varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
2351
2352 LocalRegions::MatrixKey mkey(StdRegions::eMass, LibUtilities::eSegment,
2353 *edgeExp, StdRegions::NullConstFactorMap,
2354 varcoeffs);
2355 DNekScalMat &edgemat = *edgeExp->GetLocMatrix(mkey);
2356
2357 NekVector<NekDouble> vEdgeCoeffs(order_e);
2358
2359 GetTraceToElementMap(edgeid, map, sign, v_GetTraceOrient(edgeid));
2360
2361 for (i = 0; i < order_e; ++i)
2362 {
2363 vEdgeCoeffs[i] = incoeffs[map[i]] * sign[i];
2364 }
2365
2366 vEdgeCoeffs = edgemat * vEdgeCoeffs;
2367
2368 for (i = 0; i < order_e; ++i)
2369 {
2370 coeffs[map[i]] += vEdgeCoeffs[i] * sign[i];
2371 }
2372}
2373
2374DNekMatSharedPtr Expansion2D::v_BuildVertexMatrix(
2375 const DNekScalMatSharedPtr &r_bnd)
2376{
2377 MatrixStorage storage = eFULL;
2378 DNekMatSharedPtr m_vertexmatrix;
2379
2380 int nVerts, vid1, vid2, vMap1, vMap2;
2381 NekDouble VertexValue;
2382
2383 nVerts = GetNverts();
2384
2385 m_vertexmatrix =
2386 MemoryManager<DNekMat>::AllocateSharedPtr(nVerts, nVerts, 0.0, storage);
2387 DNekMat &VertexMat = (*m_vertexmatrix);
2388
2389 for (vid1 = 0; vid1 < nVerts; ++vid1)
2390 {
2391 vMap1 = GetVertexMap(vid1);
2392
2393 for (vid2 = 0; vid2 < nVerts; ++vid2)
2394 {
2395 vMap2 = GetVertexMap(vid2);
2396 VertexValue = (*r_bnd)(vMap1, vMap2);
2397 VertexMat.SetValue(vid1, vid2, VertexValue);
2398 }
2399 }
2400
2401 return m_vertexmatrix;
2402}
2403
2409
2410Array<OneD, unsigned int> Expansion2D::GetTraceInverseBoundaryMap(int eid)
2411{
2412 int n, j;
2413 int nEdgeCoeffs;
2414 int nBndCoeffs = NumBndryCoeffs();
2415
2416 Array<OneD, unsigned int> bmap(nBndCoeffs);
2417 GetBoundaryMap(bmap);
2418
2419 // Map from full system to statically condensed system (i.e reverse
2420 // GetBoundaryMap)
2421 map<int, int> invmap;
2422 for (j = 0; j < nBndCoeffs; ++j)
2423 {
2424 invmap[bmap[j]] = j;
2425 }
2426
2427 // Number of interior edge coefficients
2428 nEdgeCoeffs = GetTraceNcoeffs(eid) - 2;
2429
2430 const SpatialDomains::Geometry2D *geom = GetGeom2D();
2431
2432 Array<OneD, unsigned int> edgemaparray(nEdgeCoeffs);
2433 Array<OneD, unsigned int> maparray(nEdgeCoeffs);
2434 Array<OneD, int> signarray(nEdgeCoeffs, 1);
2435 StdRegions::Orientation eOrient = geom->GetEorient(eid);
2436
2437 // maparray is the location of the edge within the matrix
2438 GetTraceInteriorToElementMap(eid, maparray, signarray, eOrient);
2439
2440 for (n = 0; n < nEdgeCoeffs; ++n)
2441 {
2442 edgemaparray[n] = invmap[maparray[n]];
2443 }
2444
2445 return edgemaparray;
2446}
2447
2448void Expansion2D::v_SetUpPhysNormals(const int edge)
2449{
2450 v_ComputeTraceNormal(edge);
2451}
2452
2453void Expansion2D::v_ReOrientTracePhysVals(
2454 const StdRegions::Orientation orient,
2455 const Array<OneD, const NekDouble> &in, Array<OneD, NekDouble> &out,
2456 const int nq0, [[maybe_unused]] const int nq1,
2457 [[maybe_unused]] bool Forwards)
2458{
2459 switch (orient)
2460 {
2461 case StdRegions::eForwards:
2462 // Fwd
2463 for (int i = 0; i < nq0; ++i)
2464 {
2465 out[i] = in[i];
2466 }
2467 break;
2468 case StdRegions::eBackwards:
2469 {
2470 // Bwd
2471 Vmath::Reverse(nq0, &in[0], 1, &out[0], 1);
2472 }
2473 break;
2474 default:
2475 ASSERTL0(false, "Unknown orientation");
2476 break;
2477 }
2478}
2479
2480// Compute edgenormal \cdot vector
2481Array<OneD, NekDouble> Expansion2D::GetnEdgecdotMF(
2482 const int dir, const int edge, ExpansionSharedPtr &EdgeExp_e,
2483 const Array<OneD, const Array<OneD, NekDouble>> &normals,
2484 const StdRegions::VarCoeffMap &varcoeffs)
2485{
2486 int nquad_e = EdgeExp_e->GetNumPoints(0);
2487 int coordim = GetCoordim();
2488 int nquad0 = m_base[0]->GetNumPoints();
2489 int nquad1 = m_base[1]->GetNumPoints();
2490 int nqtot = nquad0 * nquad1;
2491
2492 StdRegions::VarCoeffType MMFCoeffs[15] = {
2493 StdRegions::eVarCoeffMF1x, StdRegions::eVarCoeffMF1y,
2494 StdRegions::eVarCoeffMF1z, StdRegions::eVarCoeffMF1Div,
2495 StdRegions::eVarCoeffMF1Mag, StdRegions::eVarCoeffMF2x,
2496 StdRegions::eVarCoeffMF2y, StdRegions::eVarCoeffMF2z,
2497 StdRegions::eVarCoeffMF2Div, StdRegions::eVarCoeffMF2Mag,
2498 StdRegions::eVarCoeffMF3x, StdRegions::eVarCoeffMF3y,
2499 StdRegions::eVarCoeffMF3z, StdRegions::eVarCoeffMF3Div,
2500 StdRegions::eVarCoeffMF3Mag};
2501
2502 StdRegions::VarCoeffMap::const_iterator MFdir;
2503
2504 Array<OneD, NekDouble> ncdotMF(nqtot, 0.0);
2505 Array<OneD, NekDouble> tmp(nqtot);
2506 Array<OneD, NekDouble> tmp_e(nquad_e);
2507 for (int k = 0; k < coordim; k++)
2508 {
2509 MFdir = varcoeffs.find(MMFCoeffs[dir * 5 + k]);
2510 tmp = MFdir->second.GetValue();
2511
2512 GetPhysEdgeVarCoeffsFromElement(edge, EdgeExp_e, tmp, tmp_e);
2513
2514 Vmath::Vvtvp(nquad_e, &tmp_e[0], 1, &normals[k][0], 1, &ncdotMF[0], 1,
2515 &ncdotMF[0], 1);
2516 }
2517 return ncdotMF;
2518}
2519
2520NekDouble Expansion2D::v_VectorFlux(
2521 const Array<OneD, Array<OneD, NekDouble>> &vec)
2522{
2523 const Array<OneD, const Array<OneD, NekDouble>> &normals =
2524 GetLeftAdjacentElementExp()->GetTraceNormal(
2525 GetLeftAdjacentElementTrace());
2526
2527 int nq = GetTotPoints();
2528 Array<OneD, NekDouble> Fn(nq);
2529 Vmath::Vmul(nq, &vec[0][0], 1, &normals[0][0], 1, &Fn[0], 1);
2530 Vmath::Vvtvp(nq, &vec[1][0], 1, &normals[1][0], 1, &Fn[0], 1, &Fn[0], 1);
2531 Vmath::Vvtvp(nq, &vec[2][0], 1, &normals[2][0], 1, &Fn[0], 1, &Fn[0], 1);
2532
2533 return StdExpansion::Integral(Fn);
2534}
2535
2536void Expansion2D::v_TraceNormLen(const int traceid, NekDouble &h, NekDouble &p)
2537{
2538 SpatialDomains::Geometry *geom = GetGeom();
2539
2540 int nverts = geom->GetNumVerts();
2541
2542 // vertices on edges
2543 SpatialDomains::PointGeom ev0 = *geom->GetVertex(traceid);
2544 SpatialDomains::PointGeom ev1 = *geom->GetVertex((traceid + 1) % nverts);
2545
2546 // vertex on adjacent edge to ev0
2547 SpatialDomains::PointGeom vadj =
2548 *geom->GetVertex((traceid + (nverts - 1)) % nverts);
2549
2550 // calculate perpendicular distance of normal length
2551 // from first vertex
2552 NekDouble h1 = ev0.dist(vadj);
2553 SpatialDomains::PointGeom Dx, Dx1;
2554
2555 Dx.Sub(ev1, ev0);
2556 Dx1.Sub(vadj, ev0);
2557
2558 NekDouble d1 = Dx.dot(Dx1);
2559 NekDouble lenDx = Dx.dot(Dx);
2560 h = sqrt(h1 * h1 - d1 * d1 / lenDx);
2561
2562 // perpendicular distanace from second vertex
2563 SpatialDomains::PointGeom vadj1 = *geom->GetVertex((traceid + 2) % nverts);
2564
2565 h1 = ev1.dist(vadj1);
2566 Dx1.Sub(vadj1, ev1);
2567 d1 = Dx.dot(Dx1);
2568
2569 h = (h + sqrt(h1 * h1 - d1 * d1 / lenDx)) * 0.5;
2570
2571 int dirn = (geom->GetDir(traceid) == 0) ? 1 : 0;
2572
2573 p = (NekDouble)(GetBasisNumModes(dirn) - 1);
2574}
2575
2576/** @brief: This method gets all of the factors which are
2577 required as part of the Gradient Jump Penalty (GJP)
2578 stabilisation and involves the product of the normal and
2579 geometric factors along the element trace.
2580*/
2581void Expansion2D::v_NormalTraceDerivFactors(
2582 Array<OneD, Array<OneD, NekDouble>> &d0factors,
2583 Array<OneD, Array<OneD, NekDouble>> &d1factors,
2584 [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &d2factors)
2585{
2586 const Array<TwoD, const NekDouble> &df = m_geomFactors->GetDerivFactors();
2587 const Array<OneD, const NekDouble> &Jac = m_geomFactors->GetJac();
2588
2589 unsigned ntrace = GetNtraces();
2590
2591 if (d0factors.size() != ntrace)
2592 {
2593 d0factors = Array<OneD, Array<OneD, NekDouble>>(ntrace);
2594 d1factors = Array<OneD, Array<OneD, NekDouble>>(ntrace);
2595 }
2596
2597 Array<OneD, ExpansionSharedPtr> traceExp(ntrace);
2598 Array<OneD, unsigned> nq_edge(ntrace);
2599 unsigned nq_max = 0;
2600 for (int i = 0; i < ntrace; ++i)
2601 {
2602 // Note we are using GenTraceExp to ensure we have local trace expansion
2603 // not ont from shared trace which can happe if we use GetTraceExp since
2604 // it can be set in DisContField::SetupDG
2605 v_GenTraceExp(i, traceExp[i]);
2606 nq_edge[i] = traceExp[i]->GetTotPoints();
2607 if (d0factors[i].size() != nq_edge[i])
2608 {
2609 d0factors[i] = Array<OneD, NekDouble>(nq_edge[i]);
2610 d1factors[i] = Array<OneD, NekDouble>(nq_edge[i]);
2611 }
2612 nq_max = max(nq_max, nq_edge[i]);
2613 }
2614 Array<OneD, NekDouble> norm(nq_max);
2615
2616 const std::map<int, NormalVector> &normals = GetTraceNormals();
2617
2618 int ncoords = normals.find(0)->second.size();
2619
2620 if (m_geomFactors->GetGtype() == SpatialDomains::eDeformed)
2621 {
2622 Array<OneD, Array<OneD, NekDouble>> fac(2);
2623 for (int i = 0; i < 2; ++i)
2624 {
2625 fac[i] = Array<OneD, NekDouble>(nq_max);
2626 }
2627 Array<OneD, NekDouble> jac(nq_max);
2628 // construct local copy of df multipled by jacobian so that
2629 // interpolation is of a polynomial function to be accurate
2630 Array<OneD, Array<OneD, NekDouble>> dfdj(2 * ncoords);
2631 unsigned nqtot = GetTotPoints();
2632 for (unsigned i = 0; i < 2 * ncoords; ++i)
2633 {
2634 dfdj[i] = Array<OneD, NekDouble>(nqtot);
2635 Vmath::Vmul(nqtot, &(df[i][0]), 1, &(Jac[0]), 1, &(dfdj[i][0]), 1);
2636 }
2637
2638 // needs checking for 3D coords
2639 for (unsigned e = 0; e < ntrace; ++e)
2640 {
2641 // edge "e"
2642 v_GetLocTracePhysVals(e, traceExp[e], &(Jac[0]), jac);
2643 Vmath::Sdiv(nq_edge[e], 1.0, jac, 1, jac, 1);
2644 v_GetLocTracePhysVals(e, traceExp[e], &(dfdj[0][0]), fac[0]);
2645 v_GetLocTracePhysVals(e, traceExp[e], &(dfdj[1][0]), fac[1]);
2646
2647 norm = normals.find(e)->second[0];
2648 for (int i = 0; i < nq_edge[e]; ++i)
2649 {
2650 d0factors[e][i] = fac[0][i] * norm[i] * jac[i];
2651 d1factors[e][i] = fac[1][i] * norm[i] * jac[i];
2652 }
2653 // needs checking for 3D coords
2654 for (int n = 1; n < ncoords; ++n)
2655 {
2656 v_GetLocTracePhysVals(e, traceExp[e], &(dfdj[2 * n][0]),
2657 fac[0]);
2658 v_GetLocTracePhysVals(e, traceExp[e], &(dfdj[2 * n + 1][0]),
2659 fac[1]);
2660
2661 norm = normals.find(e)->second[n];
2662 for (int i = 0; i < nq_edge[e]; ++i)
2663 {
2664 d0factors[e][i] += fac[0][i] * norm[i] * jac[i];
2665 d1factors[e][i] += fac[1][i] * norm[i] * jac[i];
2666 }
2667 }
2668 }
2669 }
2670 else
2671 {
2672 for (unsigned e = 0; e < ntrace; ++e)
2673 {
2674 norm = normals.find(e)->second[0];
2675 for (int i = 0; i < nq_edge[e]; ++i)
2676 {
2677 d0factors[e][i] = df[0][0] * norm[i];
2678 d1factors[e][i] = df[1][0] * norm[i];
2679 }
2680
2681 for (int n = 1; n < ncoords; ++n)
2682 {
2683 norm = normals.find(e)->second[n];
2684 for (int i = 0; i < nq_edge[e]; ++i)
2685 {
2686 d0factors[e][i] += df[2 * n][0] * norm[i];
2687 d1factors[e][i] += df[2 * n + 1][0] * norm[i];
2688 }
2689 }
2690 }
2691 }
2692}
2693} // namespace Nektar::LocalRegions
ASSERTL0
#define ASSERTL0(condition, msg)
Definition ErrorUtil.hpp:216
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition ErrorUtil.hpp:250
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition ErrorUtil.hpp:273
sign
#define sign(a, b)
return the sign(b)*a
Definition Polylib.cpp:47
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition Basis.h:45
Nektar::LocalRegions::Expansion2D::v_AddEdgeNormBoundaryInt
void v_AddEdgeNormBoundaryInt(const int edge, const ExpansionSharedPtr &EdgeExp, const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray) override
Definition Expansion2D.cpp:884
Nektar::LocalRegions::Expansion2D::v_GenMatrix
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
Definition Expansion2D.cpp:1475
Nektar::LocalRegions::Expansion2D::GetGeom2D
SpatialDomains::Geometry2D * GetGeom2D() const
Definition Expansion2D.h:190
Nektar::LocalRegions::Expansion2D::m_requireNeg
std::vector< bool > m_requireNeg
Definition Expansion2D.h:114
Nektar::LocalRegions::Expansion2D::v_AddRobinTraceContribution
void v_AddRobinTraceContribution(const int traceid, const Array< OneD, const NekDouble > &primCoeffs, const Array< OneD, NekDouble > &incoeffs, Array< OneD, NekDouble > &coeffs) override
Definition Expansion2D.cpp:2336
Nektar::LocalRegions::Expansion2D::v_SetUpPhysNormals
void v_SetUpPhysNormals(const int edge) override
Definition Expansion2D.cpp:2448
Nektar::LocalRegions::Expansion2D::v_AddRobinMassMatrix
void v_AddRobinMassMatrix(const int edgeid, const Array< OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat) override
Definition Expansion2D.cpp:2206
Nektar::LocalRegions::Expansion2D::v_ReOrientTracePhysVals
void v_ReOrientTracePhysVals(const StdRegions::Orientation orient, const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out, const int nq0, const int nq1, bool Forwards) override
Definition Expansion2D.cpp:2453
Nektar::LocalRegions::Expansion2D::GetnEdgecdotMF
Array< OneD, NekDouble > GetnEdgecdotMF(const int dir, const int edge, ExpansionSharedPtr &EdgeExp_e, const Array< OneD, const Array< OneD, NekDouble > > &normals, const StdRegions::VarCoeffMap &varcoeffs)
Definition Expansion2D.cpp:2481
Nektar::LocalRegions::Expansion2D::v_PhysDirectionalDeriv
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 Expansion2D.cpp:794
Nektar::LocalRegions::Expansion2D::SetTraceToGeomOrientation
void SetTraceToGeomOrientation(Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, NekDouble > &inout)
Definition Expansion2D.cpp:1117
Nektar::LocalRegions::Expansion2D::v_BuildVertexMatrix
DNekMatSharedPtr v_BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd) override
Definition Expansion2D.cpp:2374
Nektar::LocalRegions::Expansion2D::v_VectorFlux
NekDouble v_VectorFlux(const Array< OneD, Array< OneD, NekDouble > > &vec) override
Definition Expansion2D.cpp:2520
Nektar::LocalRegions::Expansion2D::CreateMatrix
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition Expansion2D.cpp:55
Nektar::LocalRegions::Expansion2D::GetTraceInverseBoundaryMap
Array< OneD, unsigned int > GetTraceInverseBoundaryMap(int eid)
Definition Expansion2D.cpp:2410
Nektar::LocalRegions::Expansion2D::Expansion2D
Expansion2D(SpatialDomains::Geometry2D *pGeom)
Definition Expansion2D.cpp:50
Nektar::LocalRegions::Expansion2D::v_NormalTraceDerivFactors
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 Expansion2D.cpp:2581
Nektar::LocalRegions::Expansion2D::AddNormTraceInt
void AddNormTraceInt(const int dir, Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, Array< OneD, NekDouble > > &edgeCoeffs, Array< OneD, NekDouble > &outarray)
Definition Expansion2D.cpp:1196
Nektar::LocalRegions::Expansion2D::AddHDGHelmholtzEdgeTerms
void AddHDGHelmholtzEdgeTerms(const NekDouble tau, const int edge, Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, NekDouble > &edgePhys, const StdRegions::VarCoeffMap &dirForcing, Array< OneD, NekDouble > &outarray)
Definition Expansion2D.cpp:1295
Nektar::LocalRegions::Expansion2D::v_IProductWRTBase
void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculates the inner product of a given function f with the different modes of the expansion.
Definition Expansion2D.cpp:854
Nektar::LocalRegions::Expansion2D::v_DGDeriv
void v_DGDeriv(const int dir, const Array< OneD, const NekDouble > &incoeffs, Array< OneD, ExpansionSharedPtr > &EdgeExp, Array< OneD, Array< OneD, NekDouble > > &edgeCoeffs, Array< OneD, NekDouble > &out_d) override
Definition Expansion2D.cpp:2169
Nektar::LocalRegions::Expansion2D::v_PhysDeriv
void v_PhysDeriv(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculate the derivative of the physical points in a given direction.
Definition Expansion2D.cpp:703
Nektar::LocalRegions::Expansion2D::v_TraceNormLen
void v_TraceNormLen(const int traceid, NekDouble &h, NekDouble &p) override
Definition Expansion2D.cpp:2536
Nektar::LocalRegions::Expansion2D::GetPhysEdgeVarCoeffsFromElement
void GetPhysEdgeVarCoeffsFromElement(const int edge, ExpansionSharedPtr &EdgeExp, const Array< OneD, const NekDouble > &varcoeff, Array< OneD, NekDouble > &outarray)
Definition Expansion2D.cpp:1442
Nektar::LocalRegions::Expansion2D::AddEdgeBoundaryInt
void AddEdgeBoundaryInt(const int edge, ExpansionSharedPtr &EdgeExp, Array< OneD, NekDouble > &edgePhys, Array< OneD, NekDouble > &outarray, const StdRegions::VarCoeffMap &varcoeffs=StdRegions::NullVarCoeffMap)
Definition Expansion2D.cpp:1224
Nektar::LocalRegions::Expansion2D::AddHDGHelmholtzTraceTerms
void AddHDGHelmholtzTraceTerms(const NekDouble tau, const Array< OneD, const NekDouble > &inarray, Array< OneD, ExpansionSharedPtr > &EdgeExp, const StdRegions::VarCoeffMap &dirForcing, Array< OneD, NekDouble > &outarray)
Definition Expansion2D.cpp:1265
Nektar::LocalRegions::Expansion2D::v_GenTraceExp
void v_GenTraceExp(const int traceid, ExpansionSharedPtr &exp) override
Definition Expansion2D.cpp:2404
Nektar::LocalRegions::Expansion::DropLocMatrix
void DropLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition Expansion.cpp:94
Nektar::LocalRegions::Expansion::GetGeom
SpatialDomains::Geometry * GetGeom() const
Definition Expansion.cpp:171
Nektar::LocalRegions::Expansion::GetLeftAdjacentElementExp
ExpansionSharedPtr GetLeftAdjacentElementExp() const
Definition Expansion.h:531
Nektar::LocalRegions::Expansion::BuildVertexMatrix
DNekMatSharedPtr BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
Definition Expansion.cpp:105
Nektar::LocalRegions::Expansion::m_geom
SpatialDomains::Geometry * m_geom
Definition Expansion.h:306
Nektar::LocalRegions::Expansion::ComputeGmatcdotMF
void ComputeGmatcdotMF(const Array< TwoD, const NekDouble > &df, const Array< OneD, const NekDouble > &direction, Array< OneD, Array< OneD, NekDouble > > &dfdir)
Definition Expansion.cpp:783
Nektar::LocalRegions::Expansion::v_ComputeTraceNormal
virtual void v_ComputeTraceNormal(const int id)
Definition Expansion.cpp:1061
Nektar::LocalRegions::Expansion::GetMFMag
Array< OneD, NekDouble > GetMFMag(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
Definition Expansion.cpp:889
Nektar::LocalRegions::Expansion::v_GetLocTracePhysVals
virtual void v_GetLocTracePhysVals(const int trace, const StdRegions::StdExpansionSharedPtr &TraceExp, const NekDouble *inarray, Array< OneD, NekDouble > &outarray)
Definition Expansion.cpp:1025
Nektar::LocalRegions::Expansion::GetLeftAdjacentElementTrace
int GetLeftAdjacentElementTrace() const
Definition Expansion.h:544
Nektar::LocalRegions::Expansion::m_traceExp
std::map< int, ExpansionWeakPtr > m_traceExp
Definition Expansion.h:305
Nektar::LocalRegions::Expansion::AddEdgeNormBoundaryInt
void AddEdgeNormBoundaryInt(const int edge, const std::shared_ptr< Expansion > &EdgeExp, const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
Definition Expansion.cpp:118
Nektar::LocalRegions::Expansion::v_GetTraceOrient
virtual StdRegions::Orientation v_GetTraceOrient(int trace)
Definition Expansion.cpp:993
Nektar::LocalRegions::Expansion::GetTraceExp
ExpansionSharedPtr GetTraceExp(const int traceid)
Definition Expansion.h:491
Nektar::LocalRegions::Expansion::GetLocTracePhysVals
void GetLocTracePhysVals(const int trace, const StdRegions::StdExpansionSharedPtr &TraceExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition Expansion.h:222
Nektar::LocalRegions::Expansion::GetTraceOrient
StdRegions::Orientation GetTraceOrient(int trace)
Definition Expansion.h:181
Nektar::LocalRegions::Expansion::GetTraceNormals
const std::map< int, NormalVector > & GetTraceNormals(void)
Definition Expansion.cpp:264
Nektar::LocalRegions::Expansion::GetIndexMap
IndexMapValuesSharedPtr GetIndexMap(const IndexMapKey &ikey)
Definition Expansion.h:159
Nektar::LocalRegions::Expansion::GetLocMatrix
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition Expansion.cpp:88
Nektar::LocalRegions::Expansion::GetMFDiv
Array< OneD, NekDouble > GetMFDiv(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
Definition Expansion.cpp:866
Nektar::LocalRegions::Expansion::TraceNormLen
void TraceNormLen(const int traceid, NekDouble &h, NekDouble &p)
Definition Expansion.h:275
Nektar::LocalRegions::Expansion::GetTraceNormal
const NormalVector & GetTraceNormal(const int id)
Definition Expansion.cpp:250
Nektar::LocalRegions::Expansion::GetMF
Array< OneD, NekDouble > GetMF(const int dir, const int shapedim, const StdRegions::VarCoeffMap &varcoeffs)
Definition Expansion.cpp:813
Nektar::LocalRegions::Expansion::m_geomFactors
SpatialDomains::GeomFactorsUniquePtr m_geomFactors
Definition Expansion.h:307
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition NekMemoryManager.hpp:167
Nektar::SpatialDomains::Geometry2D
2D geometry information
Definition Geometry2D.h:50
Nektar::SpatialDomains::Geometry
Base class for shape geometry information.
Definition Geometry.h:84
Nektar::SpatialDomains::Geometry::GetGlobalID
int GetGlobalID(void) const
Get the ID of this object.
Definition Geometry.h:314
Nektar::SpatialDomains::Geometry::GetVertex
PointGeom * GetVertex(int i) const
Returns vertex i of this object.
Definition Geometry.h:353
Nektar::SpatialDomains::Geometry::GetCoordim
int GetCoordim() const
Return the coordinate dimension of this object (i.e. the dimension of the space in which this object ...
Definition Geometry.h:277
Nektar::SpatialDomains::Geometry::GetDir
int GetDir(const int i, const int j=0) const
Returns the element coordinate direction corresponding to a given face coordinate direction.
Definition Geometry.h:662
Nektar::SpatialDomains::Geometry::GetNumVerts
int GetNumVerts() const
Get the number of vertices of this object.
Definition Geometry.h:395
Nektar::SpatialDomains::Geometry::GetEdge
Geometry1D * GetEdge(int i) const
Returns edge i of this object.
Definition Geometry.h:361
Nektar::SpatialDomains::Geometry::GetEorient
StdRegions::Orientation GetEorient(const int i) const
Returns the orientation of edge i with respect to the ordering of edges in the standard element.
Definition Geometry.h:378
Nektar::SpatialDomains::PointGeom::Sub
void Sub(PointGeom &a, PointGeom &b)
Definition PointGeom.cpp:118
Nektar::SpatialDomains::PointGeom::dot
NekDouble dot(PointGeom &a)
retun the dot product between this and input a
Definition PointGeom.cpp:181
Nektar::SpatialDomains::PointGeom::dist
NekDouble dist(PointGeom &a)
return distance between this and input a
Definition PointGeom.cpp:174
Nektar::StdRegions::StdExpansion2D::v_IProductWRTBaseKernel
virtual void v_IProductWRTBaseKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const Array< OneD, NekDouble > &jac, const bool Deformed, bool CollDir0=false, bool CollDir1=false)=0
Nektar::StdRegions::StdExpansion2D::v_MultiplyByStdQuadratureMetric
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition StdExpansion2D.cpp:306
Nektar::StdRegions::StdExpansion::GetBoundaryMap
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition StdExpansion.h:760
Nektar::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition StdExpansion.h:168
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition StdExpansion.h:178
Nektar::StdRegions::StdExpansion::FillMode
void FillMode(const int mode, Array< OneD, NekDouble > &outarray)
This function fills the array outarray with the mode-th mode of the expansion.
Definition StdExpansion.h:561
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition StdExpansion.h:686
Nektar::StdRegions::StdExpansion::GetPointsKeys
const LibUtilities::PointsKeyVector GetPointsKeys() const
Definition StdExpansion.h:1185
Nektar::StdRegions::StdExpansion::GetLocStaticCondMatrix
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
Definition StdExpansion.h:729
Nektar::StdRegions::StdExpansion::GetVertexMap
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
Definition StdExpansion.h:770
Nektar::StdRegions::StdExpansion::IProductWRTBase
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:598
Nektar::StdRegions::StdExpansion::GetNtraces
int GetNtraces() const
Returns the number of trace elements connected to this element.
Definition StdExpansion.h:398
Nektar::StdRegions::StdExpansion::GetNverts
int GetNverts() const
This function returns the number of vertices of the expansion domain.
Definition StdExpansion.h:290
Nektar::StdRegions::StdExpansion::GetTraceToElementMap
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:775
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition StdExpansion.h:414
Nektar::StdRegions::StdExpansion::GetTraceInteriorToElementMap
void GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eForwards)
Definition StdExpansion.h:799
Nektar::StdRegions::StdExpansion::GetTraceNcoeffs
int GetTraceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th trace.
Definition StdExpansion.h:305
Nektar::StdRegions::StdExpansion::GenMatrix
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition StdExpansion.h:935
Nektar::StdRegions::StdExpansion::FwdTrans
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition StdExpansion.h:501
Nektar::StdRegions::StdExpansion::GetTraceBasisKey
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:345
Nektar::StdRegions::StdExpansion::m_stdMatrixManager
LibUtilities::NekManager< StdMatrixKey, DNekMat, StdMatrixKey::opLess > m_stdMatrixManager
Definition StdExpansion.h:1279
Nektar::StdRegions::StdExpansion::GetBasisNumModes
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
Definition StdExpansion.h:213
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition StdExpansion.h:1268
Nektar::StdRegions::StdExpansion::v_StdPhysDeriv
virtual void v_StdPhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2, Array< OneD, NekDouble > &out_d3)
Definition StdExpansion.cpp:1503
Nektar::StdRegions::StdMatrixKey::HasVarCoeffForMatrixType
bool HasVarCoeffForMatrixType(const StdRegions::MatrixType &mtype) const
Definition StdMatrixKey.h:180
Nektar::StdRegions::StdMatrixKey::GetShapeType
LibUtilities::ShapeType GetShapeType() const
Definition StdMatrixKey.h:89
Nektar::StdRegions::StdMatrixKey::GetVarFactors
const Array< OneD, const NekDouble > & GetVarFactors(const StdRegions::ConstFactorType &type) const
Definition StdMatrixKey.h:208
Nektar::StdRegions::StdMatrixKey::GetVarCoeffs
const VarCoeffMap & GetVarCoeffs() const
Definition StdMatrixKey.h:170
Nektar::StdRegions::StdMatrixKey::GetMatrixType
MatrixType GetMatrixType() const
Definition StdMatrixKey.h:84
Nektar::StdRegions::StdMatrixKey::HasVarCoeff
bool HasVarCoeff(const StdRegions::VarCoeffType &coeff) const
Definition StdMatrixKey.h:175
Nektar::StdRegions::StdMatrixKey::HasVarFactors
bool HasVarFactors(const StdRegions::ConstFactorType &type) const
Definition StdMatrixKey.h:228
Nektar::StdRegions::StdMatrixKey::GetConstFactors
const ConstFactorMap & GetConstFactors() const
Definition StdMatrixKey.h:139
Nektar::StdRegions::StdMatrixKey::GetVarCoeff
const Array< OneD, const NekDouble > & GetVarCoeff(const StdRegions::VarCoeffType &coeff) const
Definition StdMatrixKey.h:149
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition StdMatrixKey.h:125
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition StdMatrixKey.h:134
Nektar::StdRegions::StdMatrixKey::GetVarCoeffAsMap
const VarCoeffMap GetVarCoeffAsMap(const VarCoeffType &coeff) const
Definition StdMatrixKey.h:159
Blas::Daxpy
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:117
Nektar::LibUtilities::Interp1D
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
Nektar::LibUtilities::InterpCoeff1D
void InterpCoeff1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
Definition InterpCoeff.cpp:42
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition Points.h:313
Nektar::LibUtilities::eGauss_Lagrange
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
Definition BasisType.h:57
Nektar::LibUtilities::eGLL_Lagrange
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition BasisType.h:56
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition BasisType.h:48
Nektar::LocalRegions::ExpansionSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition Expansion.h:66
Nektar::LocalRegions::IndexMapValuesSharedPtr
std::shared_ptr< IndexMapValues > IndexMapValuesSharedPtr
Definition IndexMapKey.h:126
Nektar::LocalRegions::Expansion1DSharedPtr
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition Expansion1D.h:50
Nektar::LocalRegions::NormalVector
Array< OneD, Array< OneD, NekDouble > > NormalVector
Definition Expansion.h:53
Nektar::SpatialDomains::eNoGeomType
@ eNoGeomType
No type defined.
Definition SpatialDomains.hpp:51
Nektar::SpatialDomains::eDeformed
@ eDeformed
Geometry is curved or has non-constant factors.
Definition SpatialDomains.hpp:54
Nektar::StdRegions::ConstFactorMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition StdRegions.hpp:489
Nektar::StdRegions::NullConstFactorMap
static ConstFactorMap NullConstFactorMap
Definition StdRegions.hpp:490
Nektar::StdRegions::NullVarCoeffMap
static VarCoeffMap NullVarCoeffMap
Definition StdRegions.hpp:433
Nektar::StdRegions::VarCoeffMap
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
Definition StdRegions.hpp:432
Nektar::DNekScalMatSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition NekMatrixFwd.hpp:66
Nektar::DNekScalBlkMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition NekTypeDefs.hpp:79
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition BasicUtils/SharedArray.hpp:928
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition BlockMatrix.hpp:235
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition NekTypeDefs.hpp:75
Vmath::Vmul
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
Vmath::Svtvp
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
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition Vmath.hpp:292
Vmath::Vvtvp
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
Vmath::Smul
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
Vmath::Sdiv
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
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition Vmath.hpp:273
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition Vmath.hpp:54
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition Vmath.hpp:844
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition Vmath.hpp:825
std
STL namespace.
tinysimd::max
scalarT< T > max(scalarT< T > lhs, scalarT< T > rhs)
Definition scalar.hpp:305
tinysimd::min
scalarT< T > min(scalarT< T > lhs, scalarT< T > rhs)
Definition scalar.hpp:300
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition scalar.hpp:290
Generated on Sat Jul 4 2026 01:02:22 for Nektar++ by doxygen 1.9.8